If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
Department of Reproduction, Obstetrics and Herd Health, M-team and Mastitis and Milk Quality Research Unit, Ghent University, Salisburylaan 133, 9820 Merelbeke, Belgium
Department of Reproduction, Obstetrics and Herd Health, M-team and Mastitis and Milk Quality Research Unit, Ghent University, Salisburylaan 133, 9820 Merelbeke, Belgium
Department of Reproduction, Obstetrics and Herd Health, M-team and Mastitis and Milk Quality Research Unit, Ghent University, Salisburylaan 133, 9820 Merelbeke, Belgium
Milk yield dynamics during perturbations reflect how cows respond to challenges. This study investigated the characteristics of 62,406 perturbations from 16,604 lactation curves of dairy cows milked with an automated milking system at 50 Belgian, Dutch, and English farms. The unperturbed lactation curve representing the theoretical milk yield dynamics was estimated with an iterative procedure fitting a model on the daily milk yield data that was not part of a perturbation. Perturbations were defined as periods of at least 5 d of negative residuals having at least 1 day that the total daily milk production was below 80% of the estimated unperturbed lactation curve. Every perturbation was characterized and split in a development and a recovery phase. Based hereon, we calculated both the characteristics of the perturbation as a whole, and the duration, slopes, and milk losses in the phases separately. A 2-way ANOVA followed by a pairwise comparison of group means was carried out to detect differences between these characteristics in different lactation stages (early, mid-early, mid-late, and late) and parities (first, second, and third or higher). On average, 3.8 ± 1.9 (mean ± standard deviation) perturbations were detected per lactation in the first 305 d after calving, corresponding to an estimated 92.1 ± 135.8 kg of milk loss. Only 1% of the lactations had no perturbations. On average, 2.3 kg of milk was lost per day in the development phase, while the recovery phase corresponded to an average increase in milk production of 1.5 kg/d, and these phases lasted an average of 10.1 and 11.6 d, respectively. Perturbation characteristics were significantly different across parity and lactation stage groups, and early and mid-early perturbations in higher parities were found to be more severe with faster development rates, slower recovery rates, and higher milk losses. The method to characterize perturbations can be used for precision phenotyping purposes that look into the response of cows to challenges or that monitor applications (e.g., to evaluate the development and recovery of diseases and how these are affected by preventive actions or treatments).
Milk yield and production performance of dairy cows are the result of both genetic and environmental factors. In the past decades, substantial progress has been made in terms of genetic predisposition for milk yield through dedicated selection and via the increased use of advanced breeding tools such as AI and sexed semen (
). The improvement of management and environment have also contributed to a better production performance. These changes include optimization of feed and nutrition; monitoring of health, reproduction, and welfare; better veterinary practices, housing, and climate control; and improved milking routines (
Milk production is crucial for a dairy farm's profitability, and a substantial part of the economic challenge is linked to losses of milk yield and quality due to health issues (
). Both systemic and local effects of disease and inflammation contribute to these perturbations. The systemic effects are mainly caused by a decreased feed intake or a reallocation of energy to the immune system fighting an infection (
Modeling homeorhetic trajectories of milk component yields, body composition and dry-matter intake in dairy cows: Influence of parity, milk production potential and breed.
). Local effects, especially in the case of intramammary infections, can be due to damage of the blood-milk barrier in the udder caused by somatic cell inflow or pathogen toxins (
In the past, milk yield performance and dynamics were mainly analyzed in the context of genetic studies using low-frequency test-day data (e.g., for national herd health and genetic improvement programs;
). Although proven useful, these low-frequency measurements do not allow for characterizing specific perturbation dynamics such as development and recovery rates. These perturbation dynamics can, for example, give information to what extent the cow's production levels recover postchallenge. This kind of information, however, is useful to determine a cow's robustness and resilience and to compare her response to that of herd mates subjected to similar challenges (e.g., during heat stress episodes).
With high-frequency milk meter data, advanced computation, and storage capacities becoming widely available, the milk yield dynamics can be studied in more detail than before. This opens opportunities for precision phenotyping and targeted monitoring that start from insight in lactation dynamics and include the aforementioned recovery and development rates and production losses. To study perturbations, first the theoretical unperturbed lactation curve has to be estimated, for example, as proposed by
. The deviations from this unperturbed curve can then serve as the basis for detecting and characterizing the deviations in milk yield and identification of challenges, especially in absence of detailed and reliable health registers.
In this study, we hypothesized that studying the milk losses and dynamics of perturbations can help to gather insights in (1) how a dairy cow reacts to physiological and environmental challenges by characterizing specific development and recovery phases of the perturbations and (2) how milk yield dynamics during challenges vary with lactation stage and parity. To this end, we present a methodology to calculate and characterize specific milk yield dynamics in perturbations based on individual daily milk yield measurements (DMY). Using this approach, the frequency, timings, milk losses, and development and recovery rates of perturbations in the lactation curve of dairy cows were studied. Additionally, we tested the differences between the mean values of the characteristics across lactation stage and parity groups.
MATERIALS AND METHODS
Data Collection, Selection, and Preprocessing
Back-ups of the data management software of automatic milking systems (AMS) were collected at 30 dairy farms with a Lely AMS (Lely Industries N.V., Maassluis, the Netherlands) and 28 farms with a DeLaval AMS (DeLaval, Tumba, Sweden) throughout Belgium, the south of the Netherlands, and England. Using Microsoft SQL Server Management Studio 18 (SMSS, Microsoft, Redmond, WA), these back-up files were restored, and the data tables containing the historical DMY data and cow and lactation identification were extracted. All further data processing was done in Matlab R2019b (The MathWorks Inc., Natick, MA). An individual table was constructed for each farm containing unique cow and lactation identifiers including birth and calving dates and the DMY together with their corresponding measurement dates and DIM.
From these tables, lactations were selected based on the following 2 criteria: (1) the DMY data were available from before DIM 5 to at least DIM 305, and (2) no more than 2 gaps of at most 5 d each were present in the lactation curves. These gaps typically resulted from data storage issues or hardware problems and they do not represent separated or discarded milkings because these will still be registered by the robots. Only the first 305 d of each lactation were included in the analysis for standardization purposes. This avoided bias in the analysis caused by last part of the lactation curves, which can be substantially influenced by the gestation stage and feed changes toward dry-off. (
The DMY data from the DeLaval software back-ups differed from the DMY data from the Lely software back-ups because DeLaval software takes the sum of the milkings of a day as the daily value, while the Lely software corrects for the varying number of milkings between days typical for AMS data. All other differences were considered negligible. To obtain a similar variability in the DMY across all farms, the DMY data from the DeLaval farms were corrected by replacing each value with the average of the current and previous day.
Unperturbed Lactation Curve
Because health and treatment records are often incomplete and unreliable (
Quantification of antimicrobial consumption in adult cattle on dairy herds in Flanders, Belgium, and associations with udder health, milk quality, and production performance.
), and because not all challenges are noticed or treated, health and treatment registers were not considered in this study to mark perturbations. Moreover, we did not want to limit the study to perturbations caused by clinical events. Alternatively, high-frequency and consistent milk yield measurements were used to obtain a general and robust methodology for the identification of challenges and perturbations. This allowed for inclusion of all farms for which high-granularity time-series measurements were available. To this end, the deviations from unperturbed lactation curves (ULC), corresponding to the estimated theoretical production potential in the absence of perturbations, were calculated.
To calculate each lactation's ULC, we assumed that the DMY during the first 305 d after calving follow the theoretical shape of the incomplete gamma function, referred to as the Wood model (
Wood's model consists of 3 parameters; the a parameter mainly determines the scaling of the curve, whereas b and c specify the moment of peak production and slope, respectively. We preferred the Wood model over other lactation models because of its simplicity, its stability in the presence of missing data, and the computational ease which made this model the most suitable for the purpose of this study (
). To estimate the ULC that represented the milk production potential in absence of perturbations, an iterative fitting procedure was implemented to gradually remove milk yield data during perturbations. Various thresholds for the number of points excluded were tested, and upon visual evaluation of the lactation curves and their model fit (i.e., whether or not perturbations were excluded in the final model), a threshold of 1.6 × standard deviation (SD) was chosen. The methodology is visualized in Figure 1 and is as follows:
(1)
In this first iteration (i = 1), fit the nonlinear lactation curve model (ULC1; red dashed line in Figure 1) on all DMY data of the lactation (dark blue circles);
(2)
Calculate the residuals from this model, their SD (SD1) and the root mean squared error (RMSE1);
(3)
Remove all DMY data below ULC1 −1.6 × SD1 (red circles) to obtain the filtered DMY data resulting from the first iteration (i = 1). We chose 1.6 × SD as a threshold to exclude all data points outside the 95% CI;
(4)
In this next iteration (i), fit the nonlinear lactation curve model (ULCi, light blue dashed line) on the filtered DMY data resulting from the previous iteration (i-1);
(5)
Calculate the residuals from this model, their SD (SDi) and the root mean squared error (RMSEi);
(6)
Remove all DMY data below ULCi − 1.6 × SDi (light blue circles) to obtain the filtered DMY data resulting from the current iteration (i);
(7)
Repeat steps 4 to 6 until the improvement in RMSEi of the current iteration (i) compared with the previous iteration (RMSEi−1) is smaller than 0.1 kg, or after 20 iterations (final model, ULCF; depicted as the thick cyan line).
For each iteration, the parameters (a, b, c) of the previous iteration were used as the initial values for the next iteration. The nonlinear fitting procedure used the trust-region-reflective algorithm with at most 100 iterations and a lower boundary of 0 on the parameters to find the optimal solution. The advantage of using the trust-region method for optimization is that it allows for parameter constraints (
). The resulting ULCF were taken as the reference representing the total DMY that would have been achieved in absence of perturbations. Persistency of these ULCF was calculated as the decrease in predicted DMY between DIM 205 and 305 expressed in kilograms.
Figure 1Example of the initial and final unperturbed lactation curves (ULC) obtained using an iterative fitting procedure that excludes perturbations. The red dotted line represents the initial ULC1, fitted on all daily milk yields (DMY). The red dots represent the milkings below this initial fit minus 1.6 times the SD of the residuals, which are excluded during the first iteration. After 4 iterations (represented by the ULCi curves in cyan), the magenta curve (ULCF) was obtained, fitted on all the blue data and thus excluding the DMY in perturbations indicated in red and cyan.
To detect perturbations, the ULCF of each lactation was subtracted from the DMY to obtain the residuals. These residuals were expected to vary around 0 in the absence of a perturbation, while during perturbations they would be consistently negative. A perturbation was defined as a period of at least 5 successive days of negative residuals for which the DMY dropped at least once below 80% of the expected yield (ULCF), as also discussed in
). This unweighted least squares method filtered out the effect of autocorrelation between days (e.g., when the milking process was not entirely completed, resulting in a low DMY on a day and a high DMY the next day). An example of the raw and the smoothed data for 1 perturbation is shown in Figure 2. Based on the minimum of the resulting smoothed values, the timing of the maximal loss (TML; i.e., the center = d 0 of the perturbation, shown as a dashed black vertical line in Figure 2) of each perturbation was determined. In the rest of the analysis, the days before the TML were considered as the development phase, while the days after the TML represent the recovery phase.
Figure 2Example of a perturbation, its maximum milk loss, and the resulting development and recovery phase. The blue circles indicate the residuals calculated using the unperturbed lactation curve. The thick red line represents the smoothed residuals from which the minimum serves as the reference to determine the development and recovery phase. The milk loss is calculated as the sum of the residuals.
For each perturbation, several characteristics were calculated including (1) the total duration and the duration of the development and recovery phase, (2) the recovery and the development rate as the average of the differences between the smoothed DMY of all successive days in each phase (in kg/d), and (3) the milk losses. These milk losses were calculated as the sum of the nonsmoothed residuals of each perturbation, in total as well as for the development and recovery phases separately, and both expressed in absolute values (kg of milk) and relative losses (% from the predicted milk yield). For a total of 14 characteristics, the mean values were compared in a 2-way ANOVA with 2 groups (parity and lactation stage) and their interaction. The null hypothesis for this analysis was “the mean of each characteristic in each parity and lactation stage group are the same.” The parity factor had 3 levels (first, second, or third and higher), and the lactation stage had 4 levels (early = DIM 0– 63, mid-early = DIM 64–138, mid-late = DIM 139–216, late = DIM 217–305). The split for lactation stage was based on the number of perturbations detected to obtain an approximately balanced amount of perturbations per lactations stage and parity group. Most of the characteristics violated the normality assumption of the ANOVA, but because of (1) the large sample size and (2) the fact that similar tests on the log-transformed characteristics gave normally distributed residuals and the same results, this violation could be neglected. As the log-transformed results were less interpretable, they are not shown in the current manuscript. If the ANOVA pointed out significant effects of the factors, a Bonferroni mean comparison was carried out to determine significant differences between the means of the groups. Bonferroni was preferred over Tukey HSD to avoid capitalization of chance on Type I statistical errors when a high amount of means is compared. In this study, this potentially summed up to 66 comparisons, testing the differences in means for 12 combinations of lactation stage (number of groups = 4) and parity (number of groups).
RESULTS
Data Set Description
The initial data set of 58 farms (30 DeLaval, 28 Lely) contained 13,733,366 records from 20,696 cows and 50,553 lactations. From this data set, no lactations were selected for 8 farms because they either had a high number of gaps (missing data; e.g., because of technical or data storage errors) or did they did not cover a long enough time period (e.g., too recent installation of the AMS system or a reset of the historical data after upgrade to a new software version of the robot). The other 50 farms had at least 1 cow for which 305 d of uninterrupted data were available. Thirty of these farms milked with Lely, and 20 with DeLaval AMS. The resulting data set consisted of 9,681 unique cows and 16,640 unique lactations from which, respectively, 5,590; 4,543; and 6,507 from the first, second, and third or higher parity. Only the data of the first 305 d (DIM 0–305) were considered, resulting in a total data set size of 5,082,913 DMY records corresponding to 37.0% of the initial data set size. The summary statistics of the selected data set are given in Table 1, and Figure 3 shows the average and average ± 3 × SD for the lactation curve data together with 100 randomly selected individual lactation examples per parity (first, second, and third or higher). The random selection was done by sampling 100 numbers from a uniform distribution.
Table 1Summary statistics of the selected data set
Figure 3Overview of the raw daily milk yield data for each parity (first, second, and third or higher) over time. The gray background shows the mean ± 3 × SD and the orange dotted lines show the curves with the highest and lowest milk yield summed over the first 305 d, respectively. The thin blue lines show 100 randomly chosen examples per parity, while the thick red line gives the average lactation curve for that parity.
An overview of the summary statistics of the 16,640 lactations is given in Table 2. The procedure to estimate the unperturbed curve converged on average in 3.7 ± 1.5 (mean ± SD) iterations, and resulted in an average RMSE of 3.6 ± 1.4 kg. The RMSE, which is a measure for goodness-of-fit and the average deviation of all DMY data from the estimated unperturbed curve, increased with parity. The scaling parameter of the Wood models representing the ULCF increased with parity from 15.7 ± 4.6 for first parity lactations to 23.7 ± 6.5 for third and higher lactations. This translates to an average peak yield of 33.9 ± 5.7, 42.8 ± 7.2, and 46.1 ± 7.3 kg for first, second, and third and higher lactations, respectively. For first parity lactations, this peak was estimated at DIM 79.5 ± 36.5, while it was, respectively, at DIM 53.3 ± 20.2 and 50.9 ± 18.3 for second and third and higher parity lactations. The SD for peak DIM was clearly higher for primiparous than for multiparous cows. This was due to the generally flatter (i.e., more persistent) shape of the lactation curve in the first parity, making the estimation of the exact peak less robust. This was confirmed by the average persistency, expressed as the decrease in kilogram of milk from DIM 205 to DIM 305. In this period, the DMY for first parity cows decreased 4.1 ± 2.2 kg on average, corresponding to a persistency of −0.041 kg/d, while second and third and higher parity cows decreased 7.8 ± 2.7 and 9.1 ± 2.9 kg, respectively, corresponding to a persistency of −0.078 and −0.091 kg/d, respectively. Similar conclusions can be drawn from Figure 4, in which each parity had 50 randomly chosen curves plotted together with the average, minimum, and maximum estimated ULCF. The upper orange line is close to the upper limit of the gray shading, as can be expected because only data belonging to perturbations (i.e., strongly negative deviations) were deleted in the iterative procedure to obtain the ULCF. On the other hand, the lower limit of the average UCLF is clearly above the lower limit of the raw data, as the latter includes the perturbations.
Table 2Summary statistics of the unperturbed lactation curves
Figure 4Fifty randomly chosen examples of the final unperturbed lactation curves (ULCF) for each parity. The red dots represent the peak yield, while the orange dotted line gives the mean ± 3 × SD of these UCLF. The gray shading represents mean ± 3 × SD of the raw daily milk yield data.
The data set contained a total of 62,406 perturbations, defined as periods of at least 5 successive days of milk production below the ULCF, from which at least 1 d was less than 80% of the expected DMY. From these perturbations, 21,442; 16,499; and 24,465 belonged to first, second, and third and higher parity lactations, respectively. This corresponded to 3.7 ± 1.8 perturbations per lactation, on average. Only 172 lactations (1.0%) had no perturbations, while 7.3, 19.1, 22.7, 19.2, 13.8, and 16.8% of the lactations had 1, 2, 3, 4, 5, or 6 or more perturbations, respectively. The average duration of each perturbation was 19.7 ± 20.6 d. For 44.6% of the perturbations, the duration was 10 d or shorter, while 37.7% of the perturbations lasted between 10 and 30 d. The remainder of the perturbations lasted more than 30 d.
The timing of the perturbations is graphically represented in Figure 5 for each 10-d time period of the lactation showing the proportion of perturbations that start in that period. In total, 13.3% of all perturbations started in the first 10 d of lactation, and this percentage was even higher for first parity cows (15.8%). Compared with the first 10 DIM, fewer perturbations were detected between DIM 10 and 60, as new perturbations can only be detected when a previous perturbation has ended. Also, between d 60 and 120, and after d 250, a clear increase in number of perturbations was detected.
Figure 5Proportion of the start of the perturbations across lactation stage for each parity. Based on the number of perturbations, lactation stage groups were determined as follows: early = 0 to 63 DIM; mid-early = 64 to 138 DIM; mid-late = 139 to 216 DIM; late = 217 to 305 DIM. The thick gray line shows the average unperturbed lactation model for all parities.
A further overview of summary statistics for number, duration, and timing of the perturbations is given in Table 3. The lactations without recovery are considered the lactations from which the last perturbation did not get back to the ULCF, thus the residuals remained negative until at least DIM 305.
Table 3Overview of the perturbations' number, timing, and milk losses per lactation (305 d)
Item
All mean ± SD [minimum; maximum]
First parity mean ± SD [minimum; maximum]
Second parity mean ± SD [minimum; maximum]
Third parity and higher mean ± SD [minimum; maximum]
Table 4 shows the perturbation characteristics overall and per parity. On average, 92.1 kg of milk was lost per perturbation, with a very high SD of 135 kg. The distributions of the perturbation characteristics were not normally distributed, but skewed. Perturbations in higher parity cows have on average higher losses than first parity perturbations. Also, when expressed in kilograms per day, milk losses increased with parity.
Table 4Summary statistics of perturbations overall and per parity
Item
All mean ± SD [minimum; maximum]
First parity mean ± SD [minimum; maximum]
Second parity mean ± SD [minimum; maximum]
Third parity and higher mean ± SD [minimum; maximum]
In 57.2% of the perturbations, the development slope was steeper than the slope in the recovery phase, with a decrease and increase of −2.3 and +1.5 kg/d, respectively. On average, the development phase was also 1.5 d shorter than the recovery phase, with the latter lasting for 11.6 d on average. Milk losses in both phases were comparable. In absolute numbers, perturbations of the first parity are less severe (fewer milk losses, slower development and recovery rates) than perturbations in higher parity cows. Because we decided to cut off the lactations at 305 DIM, no recovery was reached for some perturbations at the end of the considered period. When leaving out these nonrecovery perturbations, the duration of the development phase decreased to 8.3 ± 10.3 d. The duration of these perturbations was on average 29 ± 19 d, which is 10 d longer than the average duration of all perturbations together (19.8 ± 20.7 d). Of all perturbations, 44.6% lasted 10 d or fewer, while, respectively, 61.2, 71.2, 77.6, and 82.2% were fewer than 15, 20, 25, or 30 d. The other 17.8% had a duration of more than 30 d.
Statistical Analysis
Grouping of the data for parity and lactation stage resulted in a data split (Table 5), and the differences between parity and lactation stage were tested for 14 different characteristics. Seven of these characteristics were linked to milk losses, while the other 7 were related to the shape and dimension of the perturbations. Their names, units, definitions, the range of the means of the groups, and the P-values of the ANOVA are given in Table 6. The degrees of freedom of each test were 2, 3, and 6 for the parity, lactation stage, and interaction between both factors, respectively. The tests had 57,333 and 62,394 error degrees of freedom, respectively, for the characteristics involving specifically the recovery phase and the other characteristics.
Table 5Number of perturbations in each combination of lactation stage and parity group
Table 6Names, units, and description of the characteristics of the individual perturbations for which the effect of parity and lactation stage were tested in a 2-way ANOVA; the range of the absolute value of the means for each group (parity × lactation stage) and respective P-values are also given
The perturbation characteristics for which fewer observations were available because no recovery was reached. The ANOVA error degrees of freedom of these characteristics were 57,333 instead of 62,394.
(kg)
Sum of the residuals in the recovery phase per perturbation
[21.2; 80.2]
<0.0001
<0.0001
<0.0001
RELLOSS (%/d)
Average loss per day in %
[11.1; 15.0]
0.0001
<0.0001
<0.0001
RELLOSS-DEV (%/d)
Average loss per day in % in the development phase
The perturbation characteristics for which fewer observations were available because no recovery was reached. The ANOVA error degrees of freedom of these characteristics were 57,333 instead of 62,394.
(d)
Duration of recovery phase
[6.9; 16.3]
<0.0001
<0.0001
<0.0001
MML (kg)
Minimum of smoothed perturbation
[7.2; 11.2]
<0.0001
<0.0001
<0.0001
RELMML (%)
% loss during minimum of smoothed perturbation
[24.3; 31.2]
<0.0001
<0.0001
<0.0001
SLOPE-DEV (kg/d)
Average slope during development phase (absolute value)
The perturbation characteristics for which fewer observations were available because no recovery was reached. The ANOVA error degrees of freedom of these characteristics were 57,333 instead of 62,394.
(kg/d)
Average slope during recovery phase
[1.3; 1.7]
<0.0001
0.0346
<0.0001
1 LS = lactation stage.
2 The perturbation characteristics for which fewer observations were available because no recovery was reached. The ANOVA error degrees of freedom of these characteristics were 57,333 instead of 62,394.
As shown in Table 6, the effects of parity, lactation stage, and their interactions were significant to highly significant (P = 0.0346 to P < 0.0001) for all 14 characteristics. The Bonferroni mean comparison tests pointed out which mean values and which trends were detected across parities and lactation stages. Moreover, the means of each group are graphically represented in Figure 6 for each characteristic. The exact ranges for the group means can be found in Table 6.
Figure 6Graphical representation of the mean values for each perturbation characteristic. Higher values have a darker color and a letter code further in the alphabet, and the grayscale is proportional to the value such that more similar hues represent means closer to each other. Groups with the same letter code were not statistically different from each other. P indicates the parity (P1 = first parity, P2 = second parity, P3 = third or higher parity), and E, ME, ML, and L represent the lactation stage classes early (0–63 DIM), mid-early (64–138 DIM), mid-late (139–216 DIM), and late (217–305 DIM), respectively. max = maximum; min = minimum.
The milk losses per perturbation expressed in kilograms increase with parity, and are higher during peak lactation (lactation stage = mid-early). For the recovery phase, the milk losses per perturbation were lower at the end of the lactation, but this should be interpreted with caution because the lactations were cut off at 305 DIM and because perturbations that did not recover were not included for the calculation of this characteristic. When considering relative losses expressed in percent loss per day, early- and late-lactation perturbations clearly showed a higher loss for each parity. Moreover, the high relative milk losses in early lactation were mainly related to the development phase, while those in late lactation were mainly linked to the recovery phase. Mid-early (during peak production) perturbations lasted longest. The recovery phases were longer in mid-early lactation for each parity, while the development phases had the longest duration for the late-lactation perturbations. In late lactation, development slopes expressed in kilograms per day were less steep, and were found to decrease with parity. In early lactation, the slope of the development phase was clearly steeper for all parities. For the recovery slope, we observed that there was a clear parity effect, and first parity animals had a slower recovery than higher parity animals. Also, the slope of the recovery phase for first parity cows increased with lactation stage, while for second and third or higher parity animals, the slope was higher in mid-early lactation. Additionally, when considering the maximal milk loss of each perturbation expressed in percent of the expected milk production, early lactation perturbations tended to have higher milk losses compared with perturbations in the later lactation stages.
DISCUSSION
Lactation Modeling and Milk Losses
In this study, we used the theoretical shape of the lactation curve of AMS-milked dairy cows in the first 305 d to identify and characterize perturbations in the DMY dynamics. To this end, we implemented an iterative, but computationally highly efficient, method to estimate the presumed ULC using the gamma function with 3 parameters (Wood model), as also proposed in
). These methods showed very similar results, and thus are not presented in this manuscript. The Wood model is shown to describe dairy cows' lactation curves well in the first 305 d, and only after the first 305 d and depending on the parity, more complex models perform better and are preferred (
). Another asset for our choice for the Wood model, in addition to the model fit and the computational advantage, was its shape stability in the presence of missing data. This was crucial for the iterative procedure we used.
Negative deviations in milk yield for at least 5 d with at least 1 d below 80% of the expected yield were considered as perturbations and included in this study. In this way, short deviations caused by, for example, hardware or measurement errors (
Milk yield and estrous behavior during eight consecutive estruses in Holstein cows fed standardized or high energy diets and grouped according to live weight changes in early lactation.
) were excluded from the analysis. The idea was to focus on perturbations caused principally by clinical and subclinical health events and challenges, but not exclusively because the exact causes were not known. The advantage of this threshold-based procedure is its simplicity and ease of implementation (
), while the disadvantage is that the value of the chosen thresholds might influence the number of perturbations detected. By testing different possibilities and assessing them visually, we identified the threshold that succeeded best in pointing out deviating (perturbed) milk yield dynamics. However, in future research, other thresholds can also be considered. Studying these specific milk yield dynamics, both at the end of the lactation or during specific health events (e.g., analyzed using research farm data for which reliable health and treatment registers are available) can be an interesting extension of the current study, but were outside the scope of the present analysis.
The effect of individual diseases on milk production in dairy cows is well described in literature, but often (1) using lower frequency data such as monthly test-day records or weekly yield aggregates (
). These studies confirm the present results in terms of order of magnitude of the daily milk losses, but did not describe specific dynamics or discriminate between development and recovery phases. In the present study, we also found that higher parities showed a larger number of perturbations, which can be due to higher incidence rates of clinical diseases (
, milk losses were higher for higher parity animals. When the milk losses were expressed in a relative way as kilogram of milk lost per expected kilogram of milk produced, these differences evaded (average milk loss per lactation was 3.2% for first and second parity and 3.6% for third and higher).
Perturbation Characteristics
Only recently did large data sets of high-frequency milk yield time series become available, thanks to the installation of sensor technology and improved data storage capacities on-farm. These data sets allow research on the individual lactation dynamics in more detail beyond averages or projected yield based on test days. Accordingly, specific characteristics can be researched from these time series (
). Distinguishing between development and recovery phases will facilitate studying the dynamics during specific types of challenges to which cows are submitted, which can provide information on how an individual cow reacts to challenges at a certain point in her lactation. The main advantage of this methodology is that it is not restricted to a single specific type of perturbation. For example, when this characterization is applied on a subset of perturbations (e.g., caused by clinical mastitis), it can provide insights into the correlations between the perturbation characteristics and the infection, immune response, cure, and treatment success. This allows not only for detection of challenges, but also to look into the reaction in milk production on these challenges over time. When applying this on, for example, known simultaneous episodes of heat stress, animals with a better or worse resistance to that challenge can be identified (e.g., for phenotyping purposes). Analyzing all of the perturbations together, we noticed that the variability across perturbation characteristics is very high and that a standard perturbation shape does not exist. As no clear correlations between the different characteristics were present, we did not report them for this study. Also, this high variability might suggest that the perturbation characteristics can allow for the identification of differences between challenges and cows to be used for monitoring or phenotyping applications (
Characterization of best linear unbiased estimates generated from national genetic evaluations of reproductive performance, survival, and milk yield in dairy cows.
). Ultimately, this would also allow for, as an example, optimizing the parity composition on-farm, considering trade-offs between production level, resiliency, and robustness across parities.
We found an average of 3.8 perturbations per parity in the first 305 d, corresponding to an average relative milk loss of 3.3% per lactation (range 0–29%, with 0 for cows without any perturbations detected). In dairy goat data,
was that they did not exclude short perturbations below 5 d, and that they also considered perturbations within perturbations, while we regarded each episode of negative residuals as a single perturbation. Accordingly, our results confirmed what was found by these authors in goat data.
The mean of each perturbation characteristic varied significantly with lactation stage and parity. In general, perturbations in early and mid-early lactation seemed more severe with higher relative losses, faster development rates, and slower recovery compared with perturbations in later lactation stages. This can partly be explained by a higher vulnerability to clinical health events (e.g., transition diseases) in the beginning of lactation, while long-lasting subclinical health events might manifest clearer in later lactation stages (
). Other explanations for slowly developing perturbations in later lactation stages are the increased occurrence of contagious intramammary infections with Staphylococcus aureus, as shown by
). Some of the perturbations will also result from technical causes such as a suboptimal model fit of the ULC, which cannot be avoided. Different lactation models and shapes exist, and not a single one fits best for all curves. However, because of the high amount of data we have, the thresholds we applied, and the advantages of the Wood model as discussed above, we assumed that these technical causes did not influence the results extensively.
Application Potential and Future Work
Information about the perturbation characteristics can be used for several purposes, including precision phenotyping of individual cows (e.g., for resilience and robustness;
). This can be a first step toward new genetic evaluation methods in which not only production levels, but also specific lactation dynamics, can be considered (
). To this end, inclusion of farm and herd effect should be considered as well. Moreover, the methodology to evaluate the perturbations' shapes can be used to identify animals that recover easily, and thus are less affected by external challenges (
Metabolic status is associated with the recovery of milk somatic cell count and milk secretion after lipopolysaccharide-induced mastitis in dairy cows.
). Therefore, characterization of perturbations is a way of compiling available data into useful information for farmers and veterinarians based on available data. To this end, looking further into correlations between the perturbations' characteristics caused by specific events can be interesting, and this provides an asset for monitoring and forecasting cure and the effect of treatment. Next to parity and lactation stage, the effects of herd, farm, seasonal effects, and diverse management practices on the characteristics and milk losses during perturbations also have to be further investigated and are identified as future work. To this end, not only the specific dynamics within a type of disease, but also individual and herd mate baselines can be considered, as well as time-series measures of complementary milk characteristics such as milk electrical conductivity or constituents (
). Moreover, it would be interesting to evaluate the effect based on the calendar date instead of lactation stage to detect challenges that affect the herd simultaneously, as is expected for heat stress. In this way, these data tools can be implemented for both individual and herd level monitoring applications, for example, to identify cows that are more or less resilient to challenges.
CONCLUSIONS
This study presents a new methodology to study perturbations in daily milk meter lactation curves of dairy cows managed on AMS herds. On average, 3.4 perturbations were detected per lactation, associated with an average milk loss of 92.1 kg. The perturbations were characterized using a third order Savitsky-Golay smoother, and their features were calculated. A distinction between the development and recovery phase was made to allow for studying specific dynamics. Mean values of the perturbation characteristics were compared between parities and lactation stage, and the differences were all found to be highly significant. Higher parity perturbations in early and mid-early lactation were shown to be more severe with faster development rates, higher milk losses, and slower recovery. This study shows the potential of using high-frequency milk meter data to characterize and identify the milk yield dynamics during challenges, which can be used for applications in on-farm monitoring and precision phenotyping of complex traits such as resilience and robustness.
ACKNOWLEDGMENTS
Ines Adriaens received funding from a KU Leuven postdoctoral mandate grant number PDM/19/132 (Leuven, Belgium). The data from this study were collected in the format of farm software back-up files by the authors, and the resulting database used is owned by the authors. Furthermore, this study is funded by VLAIO as an LA-trajectory (Brussels, Belgium), grant number HBC.2016.0774. Katherine Lumb (RAFT Solutions Ltd., Ripon, United Kingdom) contributed to the data collection of the UK farm data. The authors have not stated any conflicts of interest.
REFERENCES
Adriaens I.
Friggens N.C.
Ouweltjes W.
Scott H.
Aernouts B.
Statham J.
Productive life span and resilience rank can be predicted from on-farm first-parity sensor time series but not using a common equation across farms.
Modeling homeorhetic trajectories of milk component yields, body composition and dry-matter intake in dairy cows: Influence of parity, milk production potential and breed.
Characterization of best linear unbiased estimates generated from national genetic evaluations of reproductive performance, survival, and milk yield in dairy cows.
Milk yield and estrous behavior during eight consecutive estruses in Holstein cows fed standardized or high energy diets and grouped according to live weight changes in early lactation.
Metabolic status is associated with the recovery of milk somatic cell count and milk secretion after lipopolysaccharide-induced mastitis in dairy cows.
Quantification of antimicrobial consumption in adult cattle on dairy herds in Flanders, Belgium, and associations with udder health, milk quality, and production performance.
30961-9&id=gr1.jpg){kind=link}
30961-9&id=gr2.jpg){kind=link}
30961-9&id=gr3.jpg){kind=link}
30961-9&id=gr4.jpg){kind=link}
30961-9&id=gr5.jpg){kind=link}
30961-9&id=gr6.jpg){kind=link}