Advertisement

Between-herd variation in resilience and relations to herd performance

Open AccessPublished:November 30, 2020DOI:https://doi.org/10.3168/jds.2020-18525

      ABSTRACT

      Resilient cows are minimally affected in their functioning by infections and other disturbances, and recover quickly. Herd management is expected to have an effect on disturbances and the resilience of cows, and this effect was investigated in this study. Two resilience indicators were first recorded on individual cows. The effect of herd-year on these resilience indicators was then estimated and corrected for genetic and year-season effects. The 2 resilience indicators were the variance and the lag-1 autocorrelation of daily milk yield deviations from an expected lactation curve. Low variance and autocorrelation indicate that a cow does not fluctuate much around her expected milk yield and is, thus, subject to few disturbances, or little affected by disturbances (resilient). The herd-year estimates of the resilience indicators were estimated for 9,917 herd-year classes based on records of 227,655 primiparous cows from 2,644 herds. The herd-year estimates of the resilience indicators were then related to herd performance variables. Large differences in the herd-year estimates of the 2 resilience indicators (variance and autocorrelation) were observed between herd-years, indicating an effect of management on these traits. Furthermore, herd-year classes with a high variance tended to have a high proportion of cows with a rumen acidosis indication (r = 0.31), high SCS (r = 0.19), low fat content (r = −0.18), long calving interval (r = 0.14), low survival to second lactation (r = −0.13), large herd size (r = 0.12), low lactose content (r = −0.12), and high production (r = 0.10). These correlations support that herds with high variance are not resilient. The correlation between the variance and the proportion of cows with a rumen acidosis indication suggests that feed management may have an important effect on the variance. Herd-year classes with a high autocorrelation tended to have a high proportion of cows with a ketosis indication (r = 0.14) and a high production (r = 0.13), but a low somatic cell score (r = −0.17) and a low proportion of cows with a rumen acidosis indication (r = −0.12). These correlations suggest that high autocorrelation at herd level indicates either good or poor resilience, and is thus a poor resilience indicator. However, the combination of a high variance and a high autocorrelation is expected to indicate many fluctuations with slow recovery. In conclusion, herd management, in particular feed management, seems to affect herd resilience.

      Key words

      INTRODUCTION

      Dairy cows have to cope with a variety of environmental disturbances, such as pathogens or extreme weather. Resilient cows are minimally affected in their functioning by such disturbances, and if they are affected, they quickly recover (
      • Colditz I.G.
      • Hine B.C.
      Resilience in farm animals: biology, management, breeding and implications for animal welfare.
      ;
      • Berghof T.V.L.
      • Poppe M.
      • Mulder H.A.
      Opportunities to improve resilience in animal breeding programs.
      ). Herd management is expected to affect the resilience of cows because management is known to affect several components of resilience, such as resistance, tolerance, and recovery to and from certain infections (
      • Deng Z.
      • Koop G.
      • Lam T.J. G.M.
      • van der Lans I.A.
      • Vernooij J.C.M.
      • Hogeveen H.
      Farm-level risk factors for bovine mastitis in Dutch automatic milking dairy herds.
      ). In addition, herd management is known to affect the number of environmental disturbances on a farm, such as infection pressure (
      • Deng Z.
      • Koop G.
      • Lam T.J. G.M.
      • van der Lans I.A.
      • Vernooij J.C.M.
      • Hogeveen H.
      Farm-level risk factors for bovine mastitis in Dutch automatic milking dairy herds.
      ) and exposure to hot weather (
      • Kendall P.E.
      • Nielsen P.P.
      • Webster J.R.
      • Verkerk G.A.
      • Littlejohn R.P.
      • Matthews L.R.
      The effects of providing shade to lactating dairy cows in a temperate climate.
      ;
      • Fournel S.
      • Ouellet V.
      • Charbonneau É.
      Practices for alleviating heat stress of dairy cows in humid continental climates: A literature review.
      ).
      Recently, 2 traits were developed that indicate the resilience of individual cows. These traits were the natural log-transformed variance and lag-1 autocorrelation of daily milk yield deviations from a lactation curve (
      • Poppe M.
      • Veerkamp R.F.
      • van Pelt M.L.
      • Mulder H.A.
      Exploration of variance, autocorrelation, and skewness of deviations from lactation curves as resilience indicators for breeding.
      ). High variance indicates high variability in milk yield caused by disturbances, and high autocorrelation indicates slow return rate of milk yield to expected levels upon disturbances. Averaged per herd, variance and autocorrelation are expected to indicate the level of herd resilience, which is the ability of herd management to control resilience of the cows, and to control the number and severity of disturbances.
      Several studies used a mixed animal model to study variation in traits between herds (
      • Koivula M.
      • Nousiainen J.I.
      • Nousiainen J.
      • Mäntysaari E.A.
      Use of herd solutions from a random regression test-day model for diagnostic dairy herd management.
      ;
      • Caccamo M.
      • Veerkamp R.F.
      • De Jong G.
      • Pool M.H.
      • Petriglieri R.
      • Licitra G.
      Variance components for test-day milk, fat, and protein yield, and somatic cell score for analyzing management information.
      ;
      • Stoop W.M.
      • Van Arendonk J.A.M.
      • Heck J.M.L.
      • Van Valenberg H.J.F.
      • Bovenhuis H.
      Genetic parameters for major milk fatty acids and milk production traits of Dutch Holstein-Friesians.
      ). Such models correct for genetic effects, which makes it possible to compare differences in traits between herds due to differences in management. Differences in herd estimates for the variance and autocorrelation are expected to exist, but it is still unknown how large variation between herds is, and how differences between herds are related to herd management.
      The first aim of this study was, therefore, to investigate differences in herd resilience between herds, using the herd estimates of variance and autocorrelation of milk yield deviations measured on individual cows. The second aim was to explain the differences in herd resilience between herds by indicators of the herd (size) and its cow management (rumen acidosis, ketosis, SCS, longevity, fertility).

      MATERIALS AND METHODS

      To address the aims of this study, for each herd included in this study the following information was needed: (1) an indication of herd resilience, and (2) information about variables indicating the management and performance of the herd. First, the calculation of the herd resilience indicators will be described. After that, the calculation of herd variables related to herd management and performance based on milk production recording data will be described. For most herds, several years of data were available. Management within herds was expected to vary between years, for example due to changes in feed management or expansion of the herd. Therefore, the data of each herd were split into herd-year classes, and the resilience indicators and herd management variables were computed for each herd-year class instead of for each herd. Data editing was performed by the AWK programming language (
      • Aho A.V.
      • Kernighan B.W.
      • Weinberger P.J.
      The AWK Programming Language..
      ) and R (version 3.4.3; R Project for Statistical Computing, Vienna, Austria).

       Calculation of Resilience Indicators

      Resilience indicators were first calculated for individual cows within herds using daily milk yield records obtained by automatic milking systems (AMS). These indicators were then used to compute resilience indicators at the herd-year level using mixed model equations (see Equation [2]). The data used was the same as in
      • Poppe M.
      • Veerkamp R.F.
      • van Pelt M.L.
      • Mulder H.A.
      Exploration of variance, autocorrelation, and skewness of deviations from lactation curves as resilience indicators for breeding.
      , and originally contained 1,782,373,113 milk yield records on 1,120,550 cows obtained by AMS and conventional milking systems. The data were provided by Cooperation CRV (Arnhem, the Netherlands).
      Data preparation started off with selecting only primiparous cows that were milked by AMS, that were at least 87.5% Holstein Friesian, that were herd-book registered, and that calved after 640 d of age and before June 1, 2017. Milk yield records obtained during single AMS visits were first converted to daily milk yield records per cow. Daily milk yield records after 350 DIM were excluded. The remaining milk yield records were used to model a lactation curve for each cow that reflected her expected milk yield in the absence of disturbances. Lactation curves were fitted using fourth-order polynomial quantile regression with a 0.7 quantile:
      yt=β0+β1×t1+β2×t2+β3×t3+β4×t4+ε,
      [1]


      where yt is the observed milk yield on DIM t, tn are DIM to the power of n, where n is 1, 2, 3, or 4, βn are regression coefficients describing the relationships between tn and yt, and is the error term. Quantile regression was used instead of classical regression because it is better able to model a lactation curve without disturbances, which makes the negative deviations larger and more informative about resilience (
      • Koenker R.
      Quantile regression.
      ; Poppe et al., 2020). The model was fitted using the quantreg package (
      • Koenker R.
      quantreg: Quantile Regression. R package version 5.36.
      ) and the poly function in R. After fitting a lactation curve for each cow, the deviations in milk yield were calculated as yy. The deviations on the first and last 10 d of each lactation were removed because of poor fit of the model in the beginning and end of the lactation. If a cow had at least 50 remaining daily milk yield deviations and not more than 5% was missing, these records were used to calculate the resilience indicators, which were the natural log-transformed variance (LnVar) and the lag-1 autocorrelation (rauto) of the deviations. A low LnVar indicates that the milk yield of a cow does not fluctuate much around the expected milk yield, and thus indicates good resilience. A low rauto indicates that milk yield fluctuates quickly and independent from the previous day, and thus indicates good resilience (
      • Berghof T.V.L.
      • Poppe M.
      • Mulder H.A.
      Opportunities to improve resilience in animal breeding programs.
      ).
      The data preparation until this point was the same as in
      • Poppe M.
      • Veerkamp R.F.
      • van Pelt M.L.
      • Mulder H.A.
      Exploration of variance, autocorrelation, and skewness of deviations from lactation curves as resilience indicators for breeding.
      and resulted in 255,096 first-parity cows with an observation for LnVar and rauto. The remaining 2 data preparation steps differ in the current study compared with
      • Poppe M.
      • Veerkamp R.F.
      • van Pelt M.L.
      • Mulder H.A.
      Exploration of variance, autocorrelation, and skewness of deviations from lactation curves as resilience indicators for breeding.
      with respect to outliers and the number of cows per herd-year class. Here, records for LnVar or rauto that deviated more than 4 standard deviations from the mean were removed, which were 47 LnVar records and 227 rauto records. In
      • Poppe M.
      • Veerkamp R.F.
      • van Pelt M.L.
      • Mulder H.A.
      Exploration of variance, autocorrelation, and skewness of deviations from lactation curves as resilience indicators for breeding.
      , the same approach was used, but they computed each resilience indicator based on 4 different lactation curves. If one resilience indicator was an outlier based on all lactation curves, all other resilience indicators were also set to missing, resulting in removal of more records than in the current study. Finally, herd-year of calving classes were made, and classes with less than 10 cows were removed. After removal, 9,917 herd-year classes remained, including 2,644 herds from the years 2011 to 2017. These herd-year classes contained 227,615 cows with a record for LnVar and 227,453 cows with a record for rauto. For a detailed description of the data preparation steps, including the number of records and cows remaining after each preparation step, see Table 1 in
      • Poppe M.
      • Veerkamp R.F.
      • van Pelt M.L.
      • Mulder H.A.
      Exploration of variance, autocorrelation, and skewness of deviations from lactation curves as resilience indicators for breeding.
      .
      After calculating the LnVar and rauto of individual cows, for each herd-year the level of these traits was estimated, corrected for genetic effects and general year-season effects using a mixed model. By correcting for genetic and year-season effects, differences in LnVar or rauto between herds can be attributed only to herd management. The estimates of LnVar and rauto per herd-year were obtained using the following mixed animal model:
      yijk=HYi+YSj+ak+eijk,
      [2]


      where yijk is the LnVar or rauto of cow k in herd-year class i and year-season class j, HYi is the fixed effect of herd-year of calving i, YSj is the fixed effect of year-season of calving j (4 seasons: January–March, April–June, July–September, October–December), ak is the random genetic effect of animal k, and eijk is a random error term. The following assumptions were made about the vector of random genetic effects, a, and the vector of residuals, e:a~N(0,Aσa2), and e~N(0,Iσe2) where A is the additive genetic relationship matrix, I is the identity matrix, σa2 is the additive genetic variance, and σe2 is the residual variance. The model was applied with ASReml 4.1 (
      • Gilmour A.R.
      • Gogel B.J.
      • Cullis B.R.
      • Welham S.J.
      • Thompson R.
      ASReml user guide release 4.1 functional specification.
      ) and the pedigree included 5 generations of ancestors resulting in 758,921 animals. The estimates of LnVar and rauto for each herd-year computed by the mixed model were used in our further analyses as measures of herd resilience. Herd-year classes with a high estimate for LnVar contain cows that have on average a highly variable daily milk yield, and herd-year classes with a high estimate for rauto contain cows that have on average a slowly fluctuating milk yield.

       Calculation of Herd Performance Variables from Milk Production Registration Data

      The second aim of this study was to investigate the relationships between the resilience indicators and variables indicating herd management. Variables indicating herd management were derived from milk production registration data, birth and calving dates, and information about breed and herd book status. These variables will be referred to as herd performance variables. Herd performance variables were calculated for each herd-year class that also contained cows with resilience indicators. The data set with milk production registration records included 3-, 4-, 5-, or 6-weekly records on milk yield (kg), fat percentage (%), protein percentage (%), lactose percentage (%), SCC, ureum, and a ketosis indication based on fat-protein ratio and Fourier-transform infrared measurements of milk acetone and milk β-hydroxybutyric acid (binary, only available for first 60 d after calving;
      • Vosman J.J.
      • De Jong G.
      • Eding H.
      • Knijn H.
      Genetic evaluation for ketosis in the Netherlands based on FTIR measurements.
      ). The original data set contained 9,272,501 records on 1,065,931 first-parity cows from 4,947 herds. However, only the 2,159,817 records on the 227,655 cows in the 9,917 herd-year classes that also had resilience indicators were used.
      The herd performance variables are listed in Table 1. Part of the herd performance variables were based on the milk production recording data of the cows that also had a resilience indicator (that is, primiparous cows that met the inclusion requirements, such as being at least 87.5% Holstein Friesian and being herd-book registered). The other part of the herd performance variables, such as the herd size and average age, were based on pedigree data and calf dates of all cows in the herd-year, including older cows, crossbred cows, and so on. The herd performance variables that describe a mean of a trait were calculated by averaging the trait per cow, and then taking the mean of all cows in the herd-year. The herd performance variables that describe a proportion are the number of cows with a certain condition divided by the total number of cows that calved in the herd-year.
      Table 1The herd performance variables, units, the type of cows on which each herd performance variable was based, and the number of records for each herd performance variable
      Herd performance variablesUnitCows on which variable is based
      1 = primiparous cows with a resilience indicator in the herd-year class; 2 = all cows in the herd-year class.
      Number of herd-years with record
      Mean milk yieldKilograms19,917
      Mean fat content%19,914
      Mean protein content%19,914
      Mean lactose content%19,914
      Mean ureum content19,566
      Mean SCS19,830
      Mean calving interval from first to second lactationDays18,662
      Mean age at first calvingMonths19,917
      Proportion of cows with at least 1 elevated SCC19,830
      Proportion of cows with at least 1 rumen acidosis indication19,914
      Proportion of cows with at least 1 ketosis indication18,977
      Proportion of cows that survived to second lactation19,917
      Mean ageYears29,917
      Number of cows calved29,917
      Proportion of cows that are not 100% Holstein Friesian29,917
      Proportion of cows that are herd-book registered29,917
      1 1 = primiparous cows with a resilience indicator in the herd-year class; 2 = all cows in the herd-year class.
      The SCS was derived from SCC as (
      • CRV
      E-18 Somatic cell score with test day model.
      ):
      SCS=1,000+100×[*2log(SCC1,000)].
      [3]


      For the proportion of cows with 1 or more cases of elevated SCC, an SCC record was considered elevated if it was greater than 100,000 (
      • Windig J.J.
      • Calus M.P.L.
      • De Jong G.
      • Veerkamp R.F.
      The association between somatic cell count patterns and milk production prior to mastitis.
      ;
      • Hooijer G.A.
      • De Haas Y.
      • Horneman M.
      Uiergezondheid: Analyse van de celgetaldynamiek in het “stoplicht” van PIR-DAP.
      ). For the proportion of cows with 1 or more rumen acidosis indications, a rumen acidosis indication was given based on fat and protein content: when at least once during the lactation of a cow the fat content was lower than the protein content and at the same time lower than 4.00% (), the cow received a score of 1, and otherwise 0. For the proportion of cows that survived to the second lactation, survival to second lactation was based on the second calving date: when a second calving date was known, a cow received a score of 1 and otherwise 0. If a herd-year class contained less than 10 cows with a record for a certain trait, the record for that herd-year class was set to missing. Table 1 shows for each herd performance variable the number of remaining herd-year classes that were used for calculation of correlations (see Analyses section). For multiple regression, only herd-year classes could be used without missing herd performance variables, which were 7,828 classes (see Analyses section).

       Analyses

       Variation in LnVar and rauto Between Herd-Years and Consistency Between Years Within Herd.

      To investigate variation in LnVar and rauto between herd-years, the standard deviation, minimum, and maximum of the herd-year estimates for LnVar and rauto were computed. To investigate the association between LnVar and rauto at herd-year level, Pearson correlations between the herd-year estimates of LnVar and rauto were calculated. In addition, for both traits, Pearson correlations were calculated between herd-year estimates within herd between different years, to investigate how consistent LnVar and rauto were within herds over years.

       Associations Between Herd-Year Estimates and Herd Performance Variables.

      Pearson correlations were calculated between the herd-year estimates of LnVar and rauto and the herd performance variables derived from milk production recording data. In addition, a multiple linear regression with stepwise model selection was performed with the herd-year estimates of LnVar and rauto as dependent variables and the herd performance variables as independent variables. Such a multiple regression yields partial effects of herd performance variables on LnVar and rauto at herd-year level, conditional on the other variables involved. Stepwise model selection was performed using the StepAIC function from the MASS package (
      • Venables W.N.
      • Ripley B.D.
      Modern Applied Statistics with S.
      ) in R. The stepwise model selection could only be performed on the 7,828 herd-years where none of the herd performance variables were missing. Multiple linear regression with stepwise model selection was performed both with mean milk yield included as a fixed independent variable, and without mean milk yield included as an independent variable, to investigate how milk yield level influences the effect of other herd performance variables on resilience. After performing stepwise model selection, the relative importance of each remaining variable was obtained using the calc.relimp function with the lmg metrics from the relaimpo package (
      • Grömping U.
      Relative importance for linear regression in R: The package relaimpo.
      ). The lmg metric is based on sequential R2, but accounts for the effect of different orderings of regressors on the sequential R2 and therefore takes the average across different orderings by permuting the order of the regressors (
      • Grömping U.
      Relative importance for linear regression in R: The package relaimpo.
      ). The relative importance of a variable is thus in other words the average relative contribution of each variable to R2.

      RESULTS

       Variation in LnVar and rauto Between Herd-Years

      The LnVar and rauto differed extensively between herd-year classes. The highest LnVar estimate was more than 6 times larger than the lowest LnVar estimate (Table 2). To illustrate this: the primiparous cows of the herd-year with the smallest LnVar had on average a standard deviation of 1.19 kg in their deviations from expected yield, whereas the primiparous cows of the herd-year with the largest LnVar had on average a standard deviation of 3.66 kg in their deviations from expected yield. The highest rauto estimate was more than 4 times larger than the lowest rauto estimate (Table 2). The first-parity cows of the herd-year class with the smallest rauto had on average a correlation of 0.32 between subsequent deviations from expected yield, whereas the first-parity cows of the herd-year class with the largest rauto had on average a correlation of 0.59 between subsequent deviations from expected yield. The herd-year estimates of LnVar and rauto were both normally distributed upon visual inspection. At herd-year level, LnVar was negatively correlated with rauto (−0.34, P < 0.001), which means that herd-years with high LnVar tended to have low rauto. Thus, if in a certain herd-year one of the indicators indicated good resilience, the other indicator tended to indicate poor resilience. In summary, herds differed extensively in the 2 resilience indicators LnVar and rauto, and the 2 resilience indicators were negatively correlated at herd level.
      Table 2Descriptive statistics of the herd-year estimates of natural log-transformed variance of milk yield deviations (LnVar) and lag-1 autocorrelation of milk yield deviations (rauto)
      Herd-year estimateMeanMinimumMaximumSD
      LnVar1.3390.3812.5570.269
      rauto0.5540.1830.8170.084

       Consistency of LnVar and rauto Within Herds Between Years

      Herd-year estimates of the same herds between years were positively correlated (Table 3): if a herd had a high LnVar or rauto in a certain year, the herd also tended to have a high LnVar or rauto in other years. The correlations between years were stronger for LnVar (average: 0.58) than for rauto (average: 0.53). The correlations between years decreased with the interval between years. The average correlation between herd-year estimates of subsequent years was 0.69 for LnVar, and 0.64 for rauto (Table 3). In summary, the resilience indicators showed consistency within herds between years.
      Table 3Pearson correlations of herd-year estimates between years within herds
      Correlations for herd-year estimates of natural log-transformed variance of milk yield deviations are above the diagonal, and correlations for herd-year estimates of lag-1 autocorrelation of milk yield deviations are below the diagonal.
      YearYear
      2011201220132014201520162017
      20110.720.560.540.470.480.38
      20120.680.680.580.560.530.48
      20130.550.600.650.610.570.49
      20140.560.520.630.660.600.52
      20150.550.470.550.610.710.60
      20160.500.470.550.550.670.69
      20170.310.320.440.430.570.63
      1 Correlations for herd-year estimates of natural log-transformed variance of milk yield deviations are above the diagonal, and correlations for herd-year estimates of lag-1 autocorrelation of milk yield deviations are below the diagonal.

       Association Between Herd-Year Estimates of Resilience Indicators and Herd Performance Variables

      The herd-year classes had a mean milk yield of 25.23 kg per day with 4.42% fat and 3.59% protein, and an SCS between 1,409 and 1,726 (Table 4). The correlations between LnVar and the herd performance variables ranged from −0.18 to 0.31 (Table 5). Herd-years with high LnVar (low resilience) tended to have high production (correlation 0.10), low fat content (−0.18), low lactose content (−0.12), high proportion of cows with a rumen acidosis indication (0.31), high SCS (0.19), high proportion of cows with an elevated SCC (0.20), long calving interval (0.14), low survival to second lactation (−0.13), large herd size (0.12), and low participation in the herd book (−0.10; P < 0.001). The correlations between rauto and the herd performance variables ranged from −0.17 to 0.14 (Table 5). Herd-years with high rauto (low resilience) tended to have a high production (correlation 0.13) and a high proportion of cows with a ketosis indication (0.14), but a low proportion of cows with a rumen acidosis indication (−0.12), low SCS (−0.17), and low proportion of cows with elevated SCC (−0.15; P < 0.001). In summary, LnVar and rauto at herd-year level were correlated with several herd performance variables, and especially the correlation between the herd-year estimate of LnVar and the proportion of cows with a rumen acidosis indication was considerable (0.31).
      Table 4Descriptive statistics of the herd performance variables based on milk production recording data
      Herd performance variable
      Milk = mean milk yield; fat = mean fat content; protein = mean protein content; lactose = mean lactose content; ureum = mean ureum content; SCS = mean somatic cell score; CIN = mean calving interval from first to second lactation; CA = mean age at first calving; PropSCC = proportion of cows with at least 1 elevated somatic cell count; PropACI = proportion of cows with at least 1 rumen acidosis indication; PropKET = proportion of cows with at least 1 ketosis indication; PropSURV = proportion of cows that survived to second lactation; age = mean age; herd size = number of cows calved; PropNonHF = proportion of cows that are not 100% Holstein Friesian; PropREG = proportion of cows that are herd-book registered.
      MeanSDMinimumMaximum
      Milk (kg)25.232.9013.7037.07
      Fat (%)4.420.193.735.24
      Protein (%)3.590.113.174.05
      Lactose (%)4.640.054.414.81
      Ureum23.002.4511.3733.70
      SCS1,56445.431,4091,726
      CIN (d)402.1028.23336.50676.70
      CA (mo)25.401.3421.4534.17
      PropSCC0.670.170.081.00
      PropACI0.200.140.001.00
      PropKET0.080.090.000.90
      PropSURV0.850.130.001.00
      Age (yr)4.000.341.975.61
      Herd size108.2047.6313485
      PropNonHF0.220.180.001.00
      PropReg0.950.070.221.00
      1 Milk = mean milk yield; fat = mean fat content; protein = mean protein content; lactose = mean lactose content; ureum = mean ureum content; SCS = mean somatic cell score; CIN = mean calving interval from first to second lactation; CA = mean age at first calving; PropSCC = proportion of cows with at least 1 elevated somatic cell count; PropACI = proportion of cows with at least 1 rumen acidosis indication; PropKET = proportion of cows with at least 1 ketosis indication; PropSURV = proportion of cows that survived to second lactation; age = mean age; herd size = number of cows calved; PropNonHF = proportion of cows that are not 100% Holstein Friesian; PropREG = proportion of cows that are herd-book registered.
      Table 5Pearson correlations of herd-year estimates of the resilience indicators natural log-transformed variance of milk yield deviations (LnVar) and lag-1 autocorrelation of milk yield deviations (rauto) with the herd performance variables
      Herd performance variable
      Milk = mean milk yield; fat = mean fat content; protein = mean protein content; lactose = mean lactose content; ureum = mean ureum content; SCS = mean somatic cell score; CIN = mean calving interval from first to second lactation; CA = mean age at first calving; PropSCC = proportion of cows with at least 1 elevated somatic cell count; PropACI = proportion of cows with at least 1 rumen acidosis indication; PropKET = proportion of cows with at least 1 ketosis indication; PropSURV = proportion of cows that survived to second lactation; age = mean age; herd size = number of cows calved; PropNonHF = proportion of cows that are not 100% Holstein Friesian; PropREG = proportion of cows that are herd-book registered.
      Resilience indicator
      LnVarrauto
      Correlation with herd performance variableP-valueCorrelation with herd performance variableP-value
      Milk (kg)0.10<0.0010.13<0.001
      Fat (%)−0.18<0.0010.030.011
      Protein (%)−0.020.115−0.05<0.001
      Lactose (%)−0.12<0.0010.04<0.001
      Ureum0.010.3720.020.024
      SCS0.19<0.001−0.17<0.001
      CIN (d)0.14<0.0010.000.772
      CA (mo)0.08<0.0010.000.940
      PropSCC0.20<0.001−0.15<0.001
      PropACI0.31<0.001−0.12<0.001
      PropKET0.030.0090.14<0.001
      PropSURV−0.13<0.0010.030.003
      Age (yr)−0.05<0.0010.020.014
      Herd size0.12<0.001−0.08<0.001
      PropNonHF0.010.225−0.030.009
      PropREG−0.10<0.0010.05<0.001
      1 Milk = mean milk yield; fat = mean fat content; protein = mean protein content; lactose = mean lactose content; ureum = mean ureum content; SCS = mean somatic cell score; CIN = mean calving interval from first to second lactation; CA = mean age at first calving; PropSCC = proportion of cows with at least 1 elevated somatic cell count; PropACI = proportion of cows with at least 1 rumen acidosis indication; PropKET = proportion of cows with at least 1 ketosis indication; PropSURV = proportion of cows that survived to second lactation; age = mean age; herd size = number of cows calved; PropNonHF = proportion of cows that are not 100% Holstein Friesian; PropREG = proportion of cows that are herd-book registered.

       Partial Effects of the Herd Performance Variables on LnVar and rauto

      The model for LnVar remaining after stepwise model selection explained 20% of the variation, and when including mean milk yield this increased to 21% (Table 6). All regression coefficients were in the same direction as the correlations in Table 5, regardless if mean milk yield was included or not. Moreover, the relative importances of the explanatory variables were similar in the models with and without mean milk yield. The proportion of cows with a rumen acidosis indication had the highest relative importance in both the model without mean milk yield and the model with mean milk yield (0.37 and 0.40, respectively), followed by the herd size (0.12 and 0.11, respectively). Mean milk yield itself had a relative importance of 0.08. In summary, the proportion of cows with a rumen acidosis indication explained most variation in LnVar across herd-years, and the mean milk yield level had little effect on this association.
      Table 6Results from multiple linear regression with stepwise model selection of the herd-year estimates of natural log-transformed variance of milk yield deviations on the herd performance variables
      Regression coefficients and relative importances are shown for the herd performance variables remaining in the model after stepwise model selection, both when mean milk yield based on milk production recording (Milk) was excluded as an independent variable, and when Milk was included as an independent variable. Adjusted coefficients of determination (R2) of the models are shown in the bottom row.
      Herd performance variable
      Milk = mean milk yield; fat = mean fat content; protein = mean protein content; lactose = mean lactose content; ureum = mean ureum content; SCS = mean somatic cell score; CIN = mean calving interval from first to second lactation; CA = mean age at first calving; PropSCC = proportion of cows with at least 1 elevated somatic cell count; PropACI = proportion of cows with at least 1 rumen acidosis indication; PropKET = proportion of cows with at least 1 ketosis indication; PropSURV = proportion of cows that survived to second lactation; age = mean age; herd size = number of cows calved; PropNonHF = proportion of cows that are not 100% Holstein Friesian; PropREG = proportion of cows that are herd-book registered.
      Milk excludedMilk fixed
      Regression coefficientRelative importanceRegression coefficientRelative importance
      Milk (kg)1.16E-040.08
      Fat (%)−1.20E-050.11
      Protein (%)−8.13E-060.00
      Lactose (%)−7.24E-050.07−7.98E-050.07
      Ureum
      SCS4.58E-040.095.29E-040.09
      CIN (d)6.55E-040.046.42E-040.04
      CA (mo)7.36E-030.018.56E-030.01
      PropACI4.93E-010.375.14E-010.40
      PropKET1.33E-010.015.45E-020.00
      PropSURV−2.06E-010.04−2.18E-010.04
      PropSCC1.99E-010.102.20E-010.11
      Age (yr)−5.31E-020.02−3.79E-020.01
      Herd size8.98E-040.128.52E-040.11
      PropNonHF
      PropREG−2.62E-010.03−3.11E-010.03
      Adjusted R20.200.21
      1 Regression coefficients and relative importances are shown for the herd performance variables remaining in the model after stepwise model selection, both when mean milk yield based on milk production recording (Milk) was excluded as an independent variable, and when Milk was included as an independent variable. Adjusted coefficients of determination (R2) of the models are shown in the bottom row.
      2 Milk = mean milk yield; fat = mean fat content; protein = mean protein content; lactose = mean lactose content; ureum = mean ureum content; SCS = mean somatic cell score; CIN = mean calving interval from first to second lactation; CA = mean age at first calving; PropSCC = proportion of cows with at least 1 elevated somatic cell count; PropACI = proportion of cows with at least 1 rumen acidosis indication; PropKET = proportion of cows with at least 1 ketosis indication; PropSURV = proportion of cows that survived to second lactation; age = mean age; herd size = number of cows calved; PropNonHF = proportion of cows that are not 100% Holstein Friesian; PropREG = proportion of cows that are herd-book registered.
      The model for rauto remaining after stepwise model selection explained 8% of the variation, and when including mean milk yield this increased to 10% (Table 7). All regression coefficients, except the one for the proportion of non-Holstein Friesian cows, were in the same direction as the correlations in Table 5, regardless if mean milk yield was included or not. When including mean milk yield in the model, the relative importance of the proportion of cows with a ketosis indication decreased from 0.22 to 0.16 compared with excluding mean milk yield from the model. The relative importance of mean SCS decreased from 0.19 to 0.14. In the model without mean milk yield, the proportion of cows with a ketosis indication had the highest relative importance (0.22), whereas in the model with mean milk yield the proportion of first-parity cows with a rumen acidosis indication had the highest relative importance (0.19). Mean milk yield itself had a relative importance of 0.16. In summary, differences in rauto between herd-years could be partly explained by herd performance variables obtained from milk production registration recording, and mean milk yield affected the associations between some herd performance variables and rauto.
      Table 7Results from multiple linear regression with stepwise model selection of the herd-year estimates of lag-1 autocorrelation of milk yield deviations on the herd performance variables
      Regression coefficients and relative importances are shown for the herd performance variables remaining in the model after stepwise model selection, both when the herd performance variable mean milk yield based on milk production recording (milk) was excluded as an independent variable, and when milk was included as an independent variable. Adjusted coefficients of determination (R2) of the models are shown in the bottom row.
      Herd performance variable
      Milk = mean milk yield; fat = mean fat content; protein = mean protein content; lactose = mean lactose content; ureum = mean ureum content; SCS = mean somatic cell score; CIN = mean calving interval from first to second lactation; CA = mean age at first calving; PropSCC = proportion of cows with at least 1 elevated somatic cell count; PropACI = proportion of cows with at least 1 rumen acidosis indication; PropKET = proportion of cows with at least 1 ketosis indication; PropSURV = proportion of cows that survived to second lactation; age = mean age; herd size = number of cows calved; PropNonHF = proportion of cows that are not 100% Holstein Friesian; PropREG = proportion of cows that are herd-book registered.
      Milk excludedMilk fixed
      Regression coefficientRelative importanceRegression coefficientRelative importance
      Milk (kg)4.68E-050.16
      Fat (%)−2.83E-060.02
      Protein (%)−1.40E-060.02
      Lactose (%)5.61E-060.01
      Ureum6.54E-040.008.76E-040.00
      SCS−2.37E-040.19−2.14E-040.14
      CIN (d)1.47E-040.021.45E-040.02
      CA (mo)1.58E-030.011.64E-030.01
      PropACI−8.19E-020.18−8.53E-020.19
      PropKET1.26E-010.229.90E-020.16
      PropSURV3.46E-020.032.68E-020.02
      PropSCC−3.16E-020.14−2.21E-020.10
      Age (yr)8.57E-030.01
      Herd size−1.76E-040.14−1.98E-040.14
      PropNonHF1.58E-020.002.71E-020.01
      PropREG6.39E-020.026.41E-020.02
      Adjusted R20.080.10
      1 Regression coefficients and relative importances are shown for the herd performance variables remaining in the model after stepwise model selection, both when the herd performance variable mean milk yield based on milk production recording (milk) was excluded as an independent variable, and when milk was included as an independent variable. Adjusted coefficients of determination (R2) of the models are shown in the bottom row.
      2 Milk = mean milk yield; fat = mean fat content; protein = mean protein content; lactose = mean lactose content; ureum = mean ureum content; SCS = mean somatic cell score; CIN = mean calving interval from first to second lactation; CA = mean age at first calving; PropSCC = proportion of cows with at least 1 elevated somatic cell count; PropACI = proportion of cows with at least 1 rumen acidosis indication; PropKET = proportion of cows with at least 1 ketosis indication; PropSURV = proportion of cows that survived to second lactation; age = mean age; herd size = number of cows calved; PropNonHF = proportion of cows that are not 100% Holstein Friesian; PropREG = proportion of cows that are herd-book registered.

      DISCUSSION

      This study investigated herd differences between 2 resilience indicators, LnVar and rauto, measured on individual cows. Low LnVar of a cow indicates that milk yield does not fluctuate a lot from day to day, and low rauto indicates that milk yield of a cow recovers upon a disturbance quickly rather than slowly. Therefore, low LnVar and rauto were expected to indicate cows with good resilience (
      • Berghof T.V.L.
      • Poppe M.
      • Mulder H.A.
      Opportunities to improve resilience in animal breeding programs.
      ). Because herd-years with low herd-year estimates for the resilience indicators had low LnVar and rauto among cows, low herd-year estimates were expected to indicate good herd resilience. This study showed that herd-year estimates of LnVar and rauto differed between herd-years. Because the variation between herds was corrected for genetic effects and general year-season effects, this indicates that herd management affects LnVar and rauto. Other studies also found variation in herd estimates resulting from genetic analyses, especially for production traits and SCS (
      • Koivula M.
      • Nousiainen J.I.
      • Nousiainen J.
      • Mäntysaari E.A.
      Use of herd solutions from a random regression test-day model for diagnostic dairy herd management.
      ;
      • Caccamo M.
      • Veerkamp R.F.
      • De Jong G.
      • Pool M.H.
      • Petriglieri R.
      • Licitra G.
      Variance components for test-day milk, fat, and protein yield, and somatic cell score for analyzing management information.
      ;
      • Stoop W.M.
      • Van Arendonk J.A.M.
      • Heck J.M.L.
      • Van Valenberg H.J.F.
      • Bovenhuis H.
      Genetic parameters for major milk fatty acids and milk production traits of Dutch Holstein-Friesians.
      ). This study also showed that herds with high LnVar tended to have worse health and a lower survival than herds with low LnVar. These associations confirm the importance of herd LnVar as an indicator of herd resilience. However, the associations between rauto and herd performance variables were ambiguous, which may indicate that rauto at herd level is less directly associated with health-related traits and therefore might be perceived as not being a good indicator of resilience. In the following paragraphs, we will first discuss the meaning of herd resilience, we will then discuss the associations of LnVar and rauto with the herd performance indicators in more detail, and finally we will discuss potential application for improvement of herd management.
      Throughout the paper, we have considered LnVar and rauto at herd level as indicators of herd resilience. It is, however, important to discuss that herd resilience consists of 2 aspects, which could not be disentangled in our study. The first aspect of herd resilience is the control of number and severity of disturbances in the herd. With respect to this aspect, herds are resilient if management reduces exposure of cows to disturbances. For example, hygiene practices can reduce exposure to pathogens (
      • Deng Z.
      • Koop G.
      • Lam T.J. G.M.
      • van der Lans I.A.
      • Vernooij J.C.M.
      • Hogeveen H.
      Farm-level risk factors for bovine mastitis in Dutch automatic milking dairy herds.
      ), and roof insulation can reduce exposure to extreme weather (
      • Fournel S.
      • Ouellet V.
      • Charbonneau É.
      Practices for alleviating heat stress of dairy cows in humid continental climates: A literature review.
      ). The second aspect of herd resilience is the ability of cows in a herd to cope with disturbances. With respect to this aspect, herds are resilient if management reduces vulnerability of cows to disturbances. For example, feeding an adequate amount of vitamins and minerals can reduce vulnerability to mastitis pathogens (
      • Heinrichs A.J.
      • Costello S.S.
      • Jones C.M.
      Control of heifer mastitis by nutrition.
      ). The resilience indicators in this study capture both aspects of herd resilience simultaneously, and in the case of high LnVar or rauto it is unknown whether it is the result of a high number of disturbances or the lack of ability of the cows to respond to disturbances, or a combination of the 2 aspects. Although it is unknown from our resilience indicators which aspect leads to poor herd resilience, high LnVar and rauto can still provide a sign that management may need improvement.
      Herd resilience is expected to be affected by multiple factors and is not a direct measure of health, survival, or fertility. Nevertheless, herds with a high LnVar among cows tended to have a high proportion of cows with a rumen acidosis indication, high SCS, a high proportion of cows with elevated SCC, low survival to second lactation, long calving interval, a high proportion of cows with a ketosis indication and high age at first calving. These relations suggest that high LnVar is also indicative of the health, survival, and fertility status of the herd. It is not surprising that LnVar is indicative of health traits, because rumen acidosis (
      • Krause K.M.
      • Oetzel G.R.
      Understanding and preventing subacute ruminal acidosis in dairy herds: a review.
      ;
      • Enemark J.M.D.
      The monitoring, prevention and treatment of sub-acute ruminal acidosis (SARA): A review.
      ;
      • Humer E.
      • Aschenbach J.R.
      • Neubauer V.
      • Kröger I.
      • Khiaosa-ard R.
      • Baumgartner W.
      • Zebeli Q.
      Signals for identifying cows at risk of subacute ruminal acidosis in dairy veterinary practice.
      ), mastitis (
      • Rajala-Schultz P.J.
      • Gröhn Y.T.
      • McCulloch C.E.
      • Guard C.L.
      Effects of clinical mastitis on milk yield in dairy cows.
      ;
      • Gröhn Y.T.
      • Wilson D.J.
      • González R.N.
      • Hertl J.A.
      • Schulte H.
      • Bennett G.
      • Schukken Y.H.
      Effect of pathogen-specific clinical mastitis on milk yield in dairy cows.
      ;
      • Halasa T.
      • Nielen M.
      • De Roos A.P.W.
      • Van Hoorne R.
      • De Jong G.
      • Lam T.J. G.M.
      • Van Werven T.
      • Hogeveen H.
      Production loss due to new subclinical mastitis in Dutch dairy cows estimated with a test-day model.
      ), and ketosis (
      • Rajala-Schultz P.J.
      • Gröhn Y.T.
      • McCulloch C.E.
      Effects of milk fever, ketosis, and lameness on milk yield in dairy cows.
      ) can all lead to drops in milk yield in individual cows. A high number of cows with drops in yield in a herd will lead to increased mean LnVar at the herd level as well. In addition, the effect of health problems on LnVar in individual cows will likely be reinforced by the increased vulnerability to other disturbances once a cow has a health problem, such as subclinical ketosis (
      • Raboisson D.
      • Mounié M.
      • Maigné E.
      Diseases, reproductive performance, and changes in milk production associated with subclinical ketosis in dairy cows: A meta-analysis and review.
      ). Furthermore, if high LnVar indicates reduced health, and we assume that reduced health is related to low survival (
      • Beaudeau F.
      • Ducrocq V.
      • Fourichon C.
      • Seegers H.
      Effect of disease on length of productive life of French Holstein dairy cows assessed by survival analysis.
      ;
      • Neerhof H.J.
      • Madsen P.
      • Ducrocq V.P.
      • Vollema A.R.
      • Jensen J.
      • Korsgaard I.R.
      Relationships between mastitis and functional longevity in Danish Black and White dairy cattle estimated using survival analysis.
      ) and poor fertility (
      • Fourichon C.
      • Seegers H.
      • Malher X.
      Effect of disease on reproduction in the dairy cow: A meta-analysis.
      ;
      • Wolfenson D.
      • Leitner G.
      • Lavon Y.
      The disruptive effects of mastitis on reproduction and fertility in dairy cows.
      ), then high LnVar is also expected to be related to low survival and fertility as we observed. In addition, there may be management practices, such as feeding frequency (
      • DeVries T.J.
      Feeding behavior, feed space, and bunk design and management for adult dairy cattle.
      ), that underlie both LnVar and health and fertility traits. Although correlations between the herd performance variables and LnVar were in the expected direction, most of them were weak. However, strong correlations were not expected because LnVar indicates general resilience that is affected by many factors (
      • Putz A.M.
      • Harding J.C.S.
      • Dyck M.K.
      • Fortin F.
      • Plastow G.S.
      • Dekkers J.C.M.
      • PigGen Canada
      Novel resilience phenotypes using feed intake data from a natural disease challenge model in wean-to-finish pigs.
      ;
      • Poppe M.
      • Veerkamp R.F.
      • van Pelt M.L.
      • Mulder H.A.
      Exploration of variance, autocorrelation, and skewness of deviations from lactation curves as resilience indicators for breeding.
      ). If the correlation with, for example, SCC, would be strong, LnVar would be a mastitis indicator rather than a general resilience indicator. The fact that all correlations with health and survival traits were favorable, without any exception, is an important confirmation that LnVar indicates general herd resilience.
      The positive association between LnVar and the proportion of cows with a rumen acidosis indication was by far the most important one. This association may suggest that management practices that lead to rumen acidosis, such as feeding too many rapidly fermentable carbohydrates (
      • Krause K.M.
      • Oetzel G.R.
      Understanding and preventing subacute ruminal acidosis in dairy herds: a review.
      ;
      • Enemark J.M.D.
      The monitoring, prevention and treatment of sub-acute ruminal acidosis (SARA): A review.
      ;
      • Humer E.
      • Aschenbach J.R.
      • Neubauer V.
      • Kröger I.
      • Khiaosa-ard R.
      • Baumgartner W.
      • Zebeli Q.
      Signals for identifying cows at risk of subacute ruminal acidosis in dairy veterinary practice.
      ), are important factors contributing to an increased LnVar. This association is as expected because a decreased rumen pH can result in a reduction in feed intake, followed by an increase in feed intake (
      • Humer E.
      • Aschenbach J.R.
      • Neubauer V.
      • Kröger I.
      • Khiaosa-ard R.
      • Baumgartner W.
      • Zebeli Q.
      Signals for identifying cows at risk of subacute ruminal acidosis in dairy veterinary practice.
      ). This fluctuating feed intake pattern can lead to a variable milk yield (increased LnVar). However, in this study, we used an indicator of rumen acidosis based on an inverted fat-protein ratio, which is not the same as true rumen acidosis. Not only intake of large amounts of rapidly fermentable carbohydrates, but also grazing can increase the probability of an inverted fat-protein ratio (
      • Elgersma A.
      • Ellen G.
      • Van Der Horst H.
      • Boer H.
      • Dekker P.R.
      • Tamminga S.
      Quick changes in milk fat composition from cows after transition from fresh grass to a silage diet.
      ;
      • Couvreur S.
      • Hurtaud C.
      • Lopez C.
      • Delaby L.
      • Peyraud J.L.
      The linear relationship between the proportion of fresh grass in the cow diet, milk fatty acid composition, and butter properties.
      ). Yet, grazing does not necessarily deprive health (
      • Dijkstra J.
      • van Gastelen S.
      • Dieho K.
      • Nichols K.
      • Bannink A.
      Review: Rumen sensors: Data and interpretation for key rumen metabolic processes.
      ). Grazing is also likely to lead to a higher LnVar than TMR due to variations in grass quality, which may explain the positive association between LnVar and the proportion of cows with a rumen acidosis indication. Because of the expected association between LnVar and the application of grazing, as well as application of different grazing strategies within and between herds, it is important to properly account for grazing in future research, especially in the relationship between LnVar and the rumen acidosis indication.
      A positive correlation was observed between LnVar and mean milk yield, which was in accordance to some hypotheses, but in contrast to others. To start with the contrasting one, we may hypothesize that a high mean milk yield is related to low LnVar (negative association) because environments that lead to high milk yield have been shown to be related to good udder health (
      • Haile-Mariam M.
      • Bowman P.J.
      • Goddard M.E.
      Genetic and environmental relationship among calving interval, survival, persistency of milk yield and somatic cell count in dairy cattle.
      ;
      • Montaldo H.H.
      • Castillo-Juárez H.
      • Valencia-Posadas M.
      • Cienfuegos-Rivas E.G.
      • Ruiz-López F.J.
      Genetic and environmental parameters for milk production, udder health, and fertility traits in Mexican Holstein cows.
      ) and few drops in milk yield (
      • Windig J.J.
      • Calus M.P.L.
      • Veerkamp R.F.
      Influence of herd environment on health and fertility and their relationship with milk production.
      ). However, a hypothesis that is in accordance with our results is that LnVar and mean milk yield are positively associated because of scaling. From a statistical perspective, an increase in mean of a trait leads to a proportional increase in variance (
      • Falconer D.S.
      • Mackay T.F.C.
      Introduction to Quantitative Genetics.
      ). However, scaling cannot be the only reason for the positive association between LnVar and milk yield because herds with high LnVar also tended to have a higher coefficient of variation (SD/mean) than herds with a low LnVar (data not shown). From a feeding perspective, we may also hypothesize a positive association between LnVar and mean milk yield. Feeding large proportions of rapidly fermentable carbohydrates results in a high milk yield, but also increases the risk of rumen acidosis (
      • Krause K.M.
      • Oetzel G.R.
      Understanding and preventing subacute ruminal acidosis in dairy herds: a review.
      ;
      • Enemark J.M.D.
      The monitoring, prevention and treatment of sub-acute ruminal acidosis (SARA): A review.
      ;
      • Humer E.
      • Aschenbach J.R.
      • Neubauer V.
      • Kröger I.
      • Khiaosa-ard R.
      • Baumgartner W.
      • Zebeli Q.
      Signals for identifying cows at risk of subacute ruminal acidosis in dairy veterinary practice.
      ), which was indicated to be related to high LnVar. The observed positive correlation between LnVar and mean milk yield is probably a result of all 3 hypotheses and therefore also not very strong.
      Previous studies have suggested that high rauto indicates reduced resilience (
      • Scheffer M.
      Critical Transitions in Nature and Society..
      ;
      • Berghof T.V.L.
      • Poppe M.
      • Mulder H.A.
      Opportunities to improve resilience in animal breeding programs.
      ;
      • Poppe M.
      • Veerkamp R.F.
      • van Pelt M.L.
      • Mulder H.A.
      Exploration of variance, autocorrelation, and skewness of deviations from lactation curves as resilience indicators for breeding.
      ). The positive association between rauto and the proportion of cows with a ketosis indication suggests that high rauto among cows indicates reduced herd resilience: ketosis reduces milk yield (
      • Rajala-Schultz P.J.
      • Gröhn Y.T.
      • McCulloch C.E.
      Effects of milk fever, ketosis, and lameness on milk yield in dairy cows.
      ;
      • Raboisson D.
      • Mounié M.
      • Maigné E.
      Diseases, reproductive performance, and changes in milk production associated with subclinical ketosis in dairy cows: A meta-analysis and review.
      ), and reduces resilience to further disturbances (
      • Raboisson D.
      • Mounié M.
      • Maigné E.
      Diseases, reproductive performance, and changes in milk production associated with subclinical ketosis in dairy cows: A meta-analysis and review.
      ). However, herds with high rauto among their cows also tended to have low SCS, low incidence of cows with an elevated SCC, low proportion of cows with a rumen acidosis indication, and high survival. These associations indicate that a low instead of a high rauto is a sign of decreased resilience. Moreover, rauto and LnVar were negatively correlated (−0.34) at herd level, whereas a positive association was expected if high rauto indicated poor resilience. An explanation for the finding that poor resilience could be indicated by both high and low rauto could be that drops with slow recovery lead to increased rauto, whereas deep drops that immediately recover lead to low rauto. Both types of drops are signs of poor resilience, but lead to opposing rauto.
      Because of the ambiguous correlations of rauto with the herd performance variables, rauto by itself is not a good indicator of resilience at the herd level. However, in combination with increased LnVar, the herd-year estimates of rauto may provide information about recovery. If a herd has both a high LnVar and a high rauto, drops in milk yield occur and the recovery is slow. If a herd has a high LnVar but a low rauto, drops occur but recovery is fast. Therefore, resilience may be more severely affected in herds with high LnVar and high rauto than in herds with high LnVar but low rauto. In summary, rauto is by itself not a good resilience indicator, but an index of LnVar and rauto may be a better indicator of resilience at herd level than LnVar alone.
      The question still remains how the index of LnVar and rauto should be used as a management tool, given that AMS systems provide warning tools related to drops in milk yield as well. The most important application is the tactical and strategic benchmarking of the resilience of a herd relative to other herds and relative to the previous year. This benchmarking may assist farmers in making tactical or strategic decisions, such as changing the feeding regimen. Warnings provided by AMS have a different purpose than the resilience indicators because they act on the operational level rather than the tactical or strategic level. A warning is given when milk yield of individual cows drop or if the average herd yield drops, and action is expected immediately upon a warning. However, these warnings do not inform farmers if they have a lot of individual or herd drops in milk yield compared with other farms or compared with previous years, which the resilience indicators do. In summary, warnings provided by AMS are useful for operational management, whereas the herd resilience indicators are useful for tactical or strategic management.
      Herd-year estimates of LnVar, perhaps in combination with rauto, can inform about resilience level, but they do not indicate how herd resilience can be improved. The information from milk production recording, such as SCC, can provide clues, but many more factors are related to herd resilience (e.g., protection from weather influences, nutrition, cow density) that are not covered by milk production recording, as the low R2 of our regression models showed. Therefore, there is probably not a single solution to improve resilience, and the optimal solutions will likely differ between farms. Using their own expert knowledge of their farm, farmers can deduce which of their management practices likely contribute to decreased resilience. In some cases data from milk production recording and the AMS may help. Farmers can then decide if and how they will adjust management to improve resilience. In this sense, improving herd resilience is no different from improving other multifactorial traits, such as herd milk yield. There are also multiple methods to increase milk yield level; the optimal method will differ between farms, and some farmers accept a low milk yield level whereas others want to increase it. In summary, resilience of a farm could be a measure of the management quality and the general level of control at the farm, and provide information additional to the existing data on health and productivity level.
      This study investigated herd estimates of resilience indicators based on AMS data on a yearly basis. However, herd estimates of the resilience indicators could potentially also be calculated for farms with a conventional milking system with electronic milk meters, although this should be studied first. In addition, the herd resilience indicators could be even more useful when calculated for shorter periods, because this would show effects of herd management adjustments on herd resilience more quickly. Such an approach has been taken for monthly herd estimates of production traits and SCS by
      • Koivula M.
      • Nousiainen J.I.
      • Nousiainen J.
      • Mäntysaari E.A.
      Use of herd solutions from a random regression test-day model for diagnostic dairy herd management.
      and
      • Caccamo M.
      • Veerkamp R.F.
      • De Jong G.
      • Pool M.H.
      • Petriglieri R.
      • Licitra G.
      Variance components for test-day milk, fat, and protein yield, and somatic cell score for analyzing management information.
      . The resilience indicators could for example be provided on a monthly basis together with the milk production recording data. In some cases the data of the resilience indicators and the milk production recording data may support each other, for example if herd resilience decreased and SCC increased compared with the previous month. In other cases, the resilience indicators may warn of resilience problems, while no problems are suggested (yet) by the milk production recording data. This would be an important benefit of the resilience indicators compared with already existing tools because they allow farmers to detect decreasing resilience and intervene before any real health problems arise.

      CONCLUSIONS

      This study showed that LnVar and rauto, which were supposed to indicate resilience, varied widely between herd-year classes. Moreover, the differences between herd-year classes were related to several herd performance variables obtained from milk production recording. The associations with the herd performance variables demonstrated that low LnVar indicates good resilience at herd level, and suggest that resilience can be improved through management. The ambiguous findings for rauto suggest that this trait is less suitable as a resilience indicator at herd level than LnVar. However, for herds with high LnVar, rauto may indicate rate of recovery. In conclusion, differences in resilience indicators between herds exist, and these differences can be partly explained by herd management.

      ACKNOWLEDGMENTS

      We acknowledge the Dutch Ministry of Economic Affairs (TKI Agri & Food project 16022) and the Breed4Food partners Cobb Europe, CRV, Hendrix Genetics, and Topigs Norsvin for their financial support. In addition, we acknowledge Cooperation CRV and CRV BV for providing the data. Furthermore, we acknowledge European Union's Horizon 2020 research and innovation program (GenTORE) under grant agreement No. 727213 for their financial support. In addition, we acknowledge Yvette de Haas (Wageningen University & Research, Wageningen, the Netherlands) for her general input and Jan Dijkstra (Wageningen University & Research) for his input about rumen acidosis. The authors have not stated any conflicts of interest.

      REFERENCES

        • Aho A.V.
        • Kernighan B.W.
        • Weinberger P.J.
        The AWK Programming Language..
        Addison-Wesley Publishing Company, Reading, MA1988
        • Beaudeau F.
        • Ducrocq V.
        • Fourichon C.
        • Seegers H.
        Effect of disease on length of productive life of French Holstein dairy cows assessed by survival analysis.
        J. Dairy Sci. 1995; 78 (7738248): 103-117
        • Berghof T.V.L.
        • Poppe M.
        • Mulder H.A.
        Opportunities to improve resilience in animal breeding programs.
        Front. Genet. 2019; 9 (30693014): 692
        • Caccamo M.
        • Veerkamp R.F.
        • De Jong G.
        • Pool M.H.
        • Petriglieri R.
        • Licitra G.
        Variance components for test-day milk, fat, and protein yield, and somatic cell score for analyzing management information.
        J. Dairy Sci. 2008; 91 (18650304): 3268-3276
        • Colditz I.G.
        • Hine B.C.
        Resilience in farm animals: biology, management, breeding and implications for animal welfare.
        Anim. Prod. Sci. 2016; 56: 1961-1983
        • Couvreur S.
        • Hurtaud C.
        • Lopez C.
        • Delaby L.
        • Peyraud J.L.
        The linear relationship between the proportion of fresh grass in the cow diet, milk fatty acid composition, and butter properties.
        J. Dairy Sci. 2006; 89 (16702259): 1956-1969
        • CRV
        Handleiding MPR Voeding.
        • CRV
        E-18 Somatic cell score with test day model.
        • Deng Z.
        • Koop G.
        • Lam T.J. G.M.
        • van der Lans I.A.
        • Vernooij J.C.M.
        • Hogeveen H.
        Farm-level risk factors for bovine mastitis in Dutch automatic milking dairy herds.
        J. Dairy Sci. 2019; 102 (30852004): 4522-4535
        • DeVries T.J.
        Feeding behavior, feed space, and bunk design and management for adult dairy cattle.
        Vet. Clin. North Am. Food Anim. Pract. 2019; 35 (30686466): 61-76
        • Dijkstra J.
        • van Gastelen S.
        • Dieho K.
        • Nichols K.
        • Bannink A.
        Review: Rumen sensors: Data and interpretation for key rumen metabolic processes.
        Animal. 2020; 14 (32024561): s176-s186
        • Elgersma A.
        • Ellen G.
        • Van Der Horst H.
        • Boer H.
        • Dekker P.R.
        • Tamminga S.
        Quick changes in milk fat composition from cows after transition from fresh grass to a silage diet.
        Anim. Feed Sci. Technol. 2004; 117: 13-27
        • Enemark J.M.D.
        The monitoring, prevention and treatment of sub-acute ruminal acidosis (SARA): A review.
        Vet. J. 2008; 176 (18343172): 32-43
        • Falconer D.S.
        • Mackay T.F.C.
        Introduction to Quantitative Genetics.
        4th ed. P.E. Limited, Essex, UK1996
        • Fourichon C.
        • Seegers H.
        • Malher X.
        Effect of disease on reproduction in the dairy cow: A meta-analysis.
        Theriogenology. 2000; 53 (10968418): 1729-1759
        • Fournel S.
        • Ouellet V.
        • Charbonneau É.
        Practices for alleviating heat stress of dairy cows in humid continental climates: A literature review.
        Animals (Basel). 2017; 7 (28468329): 37
        • Gilmour A.R.
        • Gogel B.J.
        • Cullis B.R.
        • Welham S.J.
        • Thompson R.
        ASReml user guide release 4.1 functional specification.
        VSN International Ltd., Hemel Hempstead, UK2015
        • Gröhn Y.T.
        • Wilson D.J.
        • González R.N.
        • Hertl J.A.
        • Schulte H.
        • Bennett G.
        • Schukken Y.H.
        Effect of pathogen-specific clinical mastitis on milk yield in dairy cows.
        J. Dairy Sci. 2004; 87 (15377615): 3358-3374
        • Grömping U.
        Relative importance for linear regression in R: The package relaimpo.
        J. Stat. Softw. 2006; 17: 1-27
        • Haile-Mariam M.
        • Bowman P.J.
        • Goddard M.E.
        Genetic and environmental relationship among calving interval, survival, persistency of milk yield and somatic cell count in dairy cattle.
        Livest. Prod. Sci. 2003; 80: 189-200
        • Halasa T.
        • Nielen M.
        • De Roos A.P.W.
        • Van Hoorne R.
        • De Jong G.
        • Lam T.J. G.M.
        • Van Werven T.
        • Hogeveen H.
        Production loss due to new subclinical mastitis in Dutch dairy cows estimated with a test-day model.
        J. Dairy Sci. 2009; 92 (19164670): 599-606
        • Heinrichs A.J.
        • Costello S.S.
        • Jones C.M.
        Control of heifer mastitis by nutrition.
        Vet. Microbiol. 2009; 134 (18947943): 172-176
        • Hooijer G.A.
        • De Haas Y.
        • Horneman M.
        Uiergezondheid: Analyse van de celgetaldynamiek in het “stoplicht” van PIR-DAP.
        Tijdschr. Diergeneeskd. 2008; 133 (18505231): 340-341
        • Humer E.
        • Aschenbach J.R.
        • Neubauer V.
        • Kröger I.
        • Khiaosa-ard R.
        • Baumgartner W.
        • Zebeli Q.
        Signals for identifying cows at risk of subacute ruminal acidosis in dairy veterinary practice.
        J. Anim. Physiol. Anim. Nutr. (Berl.). 2018; 102 (29218772): 380-392
        • Kendall P.E.
        • Nielsen P.P.
        • Webster J.R.
        • Verkerk G.A.
        • Littlejohn R.P.
        • Matthews L.R.
        The effects of providing shade to lactating dairy cows in a temperate climate.
        Livest. Sci. 2006; 103: 148-157
        • Koenker R.
        Quantile regression.
        in: Prat A. Bonhomme S. Econometric Society Monographs. Cambridge University Press, New York, NY2005
        • Koenker R.
        quantreg: Quantile Regression. R package version 5.36.
        • Koivula M.
        • Nousiainen J.I.
        • Nousiainen J.
        • Mäntysaari E.A.
        Use of herd solutions from a random regression test-day model for diagnostic dairy herd management.
        J. Dairy Sci. 2007; 90 (17430961): 2563-2568
        • Krause K.M.
        • Oetzel G.R.
        Understanding and preventing subacute ruminal acidosis in dairy herds: a review.
        Anim. Feed Sci. Technol. 2006; 126: 215-236
        • Montaldo H.H.
        • Castillo-Juárez H.
        • Valencia-Posadas M.
        • Cienfuegos-Rivas E.G.
        • Ruiz-López F.J.
        Genetic and environmental parameters for milk production, udder health, and fertility traits in Mexican Holstein cows.
        J. Dairy Sci. 2010; 93 (20412932): 2168-2175
        • Neerhof H.J.
        • Madsen P.
        • Ducrocq V.P.
        • Vollema A.R.
        • Jensen J.
        • Korsgaard I.R.
        Relationships between mastitis and functional longevity in Danish Black and White dairy cattle estimated using survival analysis.
        J. Dairy Sci. 2000; 83 (10821581): 1064-1071
        • Poppe M.
        • Veerkamp R.F.
        • van Pelt M.L.
        • Mulder H.A.
        Exploration of variance, autocorrelation, and skewness of deviations from lactation curves as resilience indicators for breeding.
        J. Dairy Sci. 2020; 103 (31759590): 1667-1684
        • Putz A.M.
        • Harding J.C.S.
        • Dyck M.K.
        • Fortin F.
        • Plastow G.S.
        • Dekkers J.C.M.
        • PigGen Canada
        Novel resilience phenotypes using feed intake data from a natural disease challenge model in wean-to-finish pigs.
        Front. Genet. 2019; 9: 660
        • Raboisson D.
        • Mounié M.
        • Maigné E.
        Diseases, reproductive performance, and changes in milk production associated with subclinical ketosis in dairy cows: A meta-analysis and review.
        J. Dairy Sci. 2014; 97 (25306269): 7547-7563
        • Rajala-Schultz P.J.
        • Gröhn Y.T.
        • McCulloch C.E.
        Effects of milk fever, ketosis, and lameness on milk yield in dairy cows.
        J. Dairy Sci. 1999; 82 (10068950): 288-294
        • Rajala-Schultz P.J.
        • Gröhn Y.T.
        • McCulloch C.E.
        • Guard C.L.
        Effects of clinical mastitis on milk yield in dairy cows.
        J. Dairy Sci. 1999; 82 (10386307): 1213-1220
        • Scheffer M.
        Critical Transitions in Nature and Society..
        Princeton University Press, Princeton, NJ2009
        • Stoop W.M.
        • Van Arendonk J.A.M.
        • Heck J.M.L.
        • Van Valenberg H.J.F.
        • Bovenhuis H.
        Genetic parameters for major milk fatty acids and milk production traits of Dutch Holstein-Friesians.
        J. Dairy Sci. 2008; 91 (18096963): 385-394
        • Venables W.N.
        • Ripley B.D.
        Modern Applied Statistics with S.
        Fourth Edition. Springer, New York, NY2002
        • Vosman J.J.
        • De Jong G.
        • Eding H.
        • Knijn H.
        Genetic evaluation for ketosis in the Netherlands based on FTIR measurements.
        Interbull Bull. 2015; 49: 1-5
        • Windig J.J.
        • Calus M.P.L.
        • De Jong G.
        • Veerkamp R.F.
        The association between somatic cell count patterns and milk production prior to mastitis.
        Livest. Prod. Sci. 2005; 96: 291-299
        • Windig J.J.
        • Calus M.P.L.
        • Veerkamp R.F.
        Influence of herd environment on health and fertility and their relationship with milk production.
        J. Dairy Sci. 2005; 88 (15591398): 335-347
        • Wolfenson D.
        • Leitner G.
        • Lavon Y.
        The disruptive effects of mastitis on reproduction and fertility in dairy cows.
        Ital. J. Anim. Sci. 2015; 14: 650-654