Genetic analysis of lactation consistency in U.S. Holsteins using temporal variation in daily milk weights

The ability of a dairy cow to perform reliably over time is an interesting trait to include in dairy cattle breeding programs aimed at improving dairy cow resilience. Consistency, defined as the quality of performing as expected each day of the lactation, could be highly associated with resilience, defined as animal’s ability to maintain health and performance in the presence of environmental challenges, including pathogens, heat waves, and nutritional changes. A total of 51,415,022 daily milk weights collected from 2018 to 2023 were provided for 255,191 multiparous Holstein cows milked 3 times daily in conventional parlor systems on farms in 32 states. The temporal variance (TempVar) of milk yield from 5 to 305 d postpartum was computed as the log-transformed variance of daily deviations be-tween observed and expected individual milk weights. Lower values of TempVar imply smaller day-to-day deviations from expectations, indicating consistent performance, whereas larger values indicate inconsistent performance. Expected daily milk weights were computed using 3 nonparametric and parametric regression models: 1) LOESS regression with a 0.75 span; 2) polynomial quantile regression using the median (0.5 quantile), and 3) polynomial quantile regression using a 0.7 quantile. The univariate statistical model included age at first calving and herd-year-season as fixed effects and cow as a random effect. Heritability estimates (standard errors) of TempVar phenotypes calculated over the entire lactation ranged between 0.227 (0.011) and 0.237 (0.011), demonstrating that cows are genetically predisposed to display consistent or inconsistent performance. Estimated genetic correlations calculated using a multiple trait model between TempVar traits and between lactations were high (>0.95), indicating TempVar is repeatable across lactations and robust to the model used to compute expected daily milk yield. Higher TempVar phenotypes reflect more variation in performance, hence greater inconsistency, which is undesirable. Therefore, correlations between predicted transmitting abilities (PTAs) for TempVar and milk yield of 0.57 indicate that high-producing cows exhibit more day-to-day variation in performance. Correlations with productive life and livability were −0.38 and −0.48, respectively. Correlations between PTAs for TempVar and those of postpartum health traits were also negative, ranging from −0.41 to −0.08. Given that health traits are derived from disease resistance measurements, and higher health trait PTAs are preferred, our results indicate that more consistent cows tend to have fewer health problems and greater longevity. Overall, our findings suggest that temporal variation in daily milk weights can be used to identify consistent animals that maintain expected performance throughout the lactation, which will enable selection for greater resilience to management and environmental perturbations.


INTRODUCTION
The dairy industry has made enormous gains in production efficiency through improvements in genetic selection, nutrition, and management (Capper et al., 2009;Capper andCady., 2020, Guinan et al., 2023).Historically, genetic selection has emphasized increased milk production per cow, and more recently it has also focused on improving fertility, disease resistance and feed efficiency (Cole et al., 2021;VanRaden et al., 2021).Modern intensive farming systems have prioritized the average performance of an animal in optimal conditions, while ignoring the animal's ability to perform in variable or suboptimal conditions.Such suboptimal conditions could include fluctuating weather conditions, feed quality and composition, labor shortages, disease outbreaks, and animal movements, among others (Lacatera, 2019).There is a widespread recognition of the need to shift emphasis toward enhancing an animal's resilience and robustness in response to increased Genetic analysis of lactation consistency in U.S. Holsteins using temporal variation in daily milk weights Fiona L. Guinan,1* Robert H. Fourdraine,2 Francisco Peñagaricano,1 and Kent A. Weigel 1  environmental and nutritional variability (Friggens et al., 2017).
Resiliency is a measure of an animal's capacity to bounce back to normal functioning or maintain specific functions in the face of environmental or management disturbances (Scheffer et al., 2018, Colditz andHine, 2016).Resilience is a key contributor to the economic, social, and environmental sustainability of animal agriculture.Numerous approaches to assess resilience in animal agriculture have been explored, comprising 2 fundamental components: specific and generalized resilience.Specific resilience involves identifying animals resilient to particular challenges, such as heat stress.On the other hand, generalized resilience encompasses an animal's capacity to sustain functionality across diverse challenges, showcasing adaptability and robustness in various conditions.Work on specific resilience includes Nguyen et al. (2016) and Misztal (2017) who proposed selection for enhanced ability to cope with thermal stress.Variation exists among animals in terms of the onset and rate of milk production decline when faced with heat stress.However, the practical adoption of heat stress selection is challenging due to the requirement of computing numerous estimated breeding values (EBV) for each sire, corresponding to specific weather conditions in various regions.This emphasizes the importance of enhancing generalized resilience, rather than calculating distinct resilience EBV for every conceivable type of disturbance.Generalized resilience has been shown to be moderately heritable in different species, when calculated using metrics such as the log-transformed variance of daily deviations in milk weight (Poppe et al., 2020), egg production (Bedere et al., 2022), body weight (Berghof et al., 2019) and feed intake (Putz et al., 2019).
Resilience is a multifaceted trait that can be broken down into simpler components.Our hypothesis is that temporal variation in performance is an economically important phenotype that can serve as a metric of generalized resilience.Animals that demonstrate consistency, defined as a level of performance that varies minimally over time, relative to expectations, may tend to be trouble-free and require few management and veterinary interventions.Cows exhibit varying levels of consistency or inconsistency in daily milk yield relative to expectations.Figure 1 illustrates this difference, where 2 cows belonging to the same contemporary group are plotted.The consistent cow's actual daily milk weights closely align with the fitted curves, resulting in a temporal variance (TempVar) phenotype of 2.9.In contrast, the inconsistent cow shows significant fluctuations in milk production throughout the lactation period, leading to greater variation in production and a TempVar phenotype of 4.7, which we define as inconsistent performance.Consistency of performance is a key component of economic sustainability on dairy operations because it can help farmers reduce disease occurrences, optimize labor utilization, minimize milk losses, and avoid veterinary interventions.High frequency data, in the form of daily milk weights, are routinely collected on many farms.These records provide a unique opportunity to capture phenotypes that will enable selection for greater consistency while extracting added value from the data.
The aim of this study was to estimate genetic parameters of TempVar phenotypes on US dairy farms by computing log-transformed variances of daily deviations between observed and expected milk weights.Alternative methods to calculate expected daily milk yields were considered, and heritability and genetic correlation parameters were estimated within and between parities, as well as within and between different stages of the lactation.Lastly, relationships between sires' PTAs for TempVar of daily milk yield and their PTAs for production, health, and longevity traits were assessed to unravel potential correlated responses to selection for greater milk yield consistency.

Data
Data were provided by Dairy Records Management Systems (Raleigh, NC) and were extracted from PCDART on farm management software.Data were appended to a database built using the RSQLite package in R (version 3.6.0).Individual milk weights are stored for up to 100 d while daily milk weights are stored for up to 300 d on the farm, so historical data are limited, and it was necessary to upload data from participating herds monthly and aggregate these data over time to build our SQLite research database.Data were limited to cows milked 3 times per day from 2018 to 2023 and estimated daily milk weights corresponding to days with missing values were removed.Cows milked by automatic milking systems (AMS) were excluded from the analysis.Outliers were identified and removed after decomposing the seasonal trend of the lactation curve with the Multiple Seasonal Trend decomposition method using the function tsclean from the forecast package in R (Hyndman et al., 2023).Tsclean is a robust method to identify outliers in a univariate time series analysis using a modified Z score, in which outliers are identified based on their distance from the median (Hyndman et al., 2023).Herds were required to participate in Dairy Herd Improvement (DHI) milk recording and recorded breed of cow was restricted to Holstein.After summing individual milk weights, a to-  1).

Fitting lactation curves to predict daily milk weights
The first step in computing TempVar phenotypes was to fit flexible, nonlinear lactation curves to the daily milk weights of individual cows to obtain predicted values, which then formed the basis for computing daily milk yield deviations.We conducted a univariate analysis using DIM as the predictor and daily milk weight as the response variable.Cows were required to have at least 100 daily milk weights within a lactation period to fit a lactation curve.Three different lactation curves were fitted to investigate if the fitting method impacts either the estimation of genetic parameters or the estimation of breeding values.
The first model-based method used to fit lactation curves was LOESS (also known as local polynomial regression or locally weighted scatterplot smoothing) using a 0.75 span parameter and a quadratic polynomial using the loess function in R (R Project for Statistical Computing, Vienna, Austria).The span parameter determines the size of the local regression window and therefore the scope of the neighborhood considered for each data point.A larger span value corresponds to a wider window and includes more data points in the local neighborhood, whereas a smaller span focuses on a narrower region around each point and is therefore  more sensitive to short-term changes.The second and third model-based methods used to fit lactation curves were fourth-order polynomial quantile regression models using either 0.5 or 0.7 quantile parameters: where y t is the observed milk yield for a given cow on DIM t of a given lactation, and ε is the error term.The quantreg package in R (version 3.6.0)(Koenker, 2022), along with the poly function, were utilized for conducting the quantile regression analyses.The 0.7 quantile was used as an indicator of a cow's potential milk production in the absence of disturbances, thereby reflecting expected performance in an optimum production environment (Poppe et al., 2020).

Calculating temporal variance phenotypes
TempVar phenotypes were calculated in 2 steps.First, by measuring the deviations between a cow's actual daily milk weights and her expected daily milk weights across the trajectory of her lactation, with expected values provided by the 3 lactation curve models described above.In the second step, the variances of these daily deviations were calculated.However, due to skewness of the distribution of variances across individual cows and lactations, a log-transformation was applied to variances to derive 3 TempVar phenotypes, namely LnVar_loess, LnVar_qr05, and LnVar_qr07, for each method and each cow in step 2, as follows: where i is the individual cow, j indicates parity and k represents DIM between 5 and 305.Thus, a consistent cow is defined by low temporal variation in actual versus predicted daily milk production throughout the lactation, and an inconsistent cow is defined by high temporal variation in actual versus predicted daily milk yield throughout the lactation.

Univariate analysis of temporal variance of milk yield in first lactation
Our initial analysis, which was restricted to first parity cows, included 20,787,272 daily milk weights from 102,216 cows in 213 herds in 30 states.Variance components and genetic parameter estimates were obtained using the AIREMLF90 software (Aguilar et al., 2018).
TempVar phenotypes were analyzed using the following model: where y ijkl is the TempVar phenotype, AFC i is the fixed effect of age at first calving (6 levels; < = 22, 23-24, 25-26, 27-28, 29-30, 30+), HYS j is the fixed effect of herd-year-season of calving (2,347 levels, with a minimum of 5 observations per level), a k is the random effect of cow with 102,216 levels distributed as ) and e ijk is the random residual effect distrib-

Univariate analysis of temporal variance of milk yield by stage of lactation
After fitting individual lactation curves for first lactation cows using LOESS and polynomial quantile regression, 3 different lactation stages based on DIM were considered to reflect early, mid, and late lactation.Early lactation ranged from 5 to 50 DIM, mid lactation ranged from 51 to 200 DIM, and late lactation ranged from 201 to 305 DIM.Individual cows were required to have daily milk weights spanning at least 50% of the period to be included in the analysis.In other words, at least 22 daily milk weights were required in early lactation, 75 daily milk weights were required in mid lactation, and 52 daily milk weights were required in late lactation.Consequently, a specific first parity cow could be included in the analysis for one, 2, or the 3 periods.After edits, 66,297 cows were used to estimate genetic parameters in early lactation, 85,445 cows in mid lactation, and 71,673 cows in late lactation.Variance components and heritability estimates for LnVar_loess, LnVar_qr05, and LnVar_qr07 were calculated within each lactation period using a univariate model, and relationships across periods were assessed using Pearson correlations of sire PTAs in different periods.

Analyses of temporal variance of milk yield across lactations
The edits described previously were subsequently applied to daily milk yield data of second and third parity cows to assess relationships in the TempVar of milk yield across parities.Cows were not required to have records in all 3 lactations because, given the structure of our database, the number of cows with records in multiple parities was limited.In the data set, there were 36,589 cows with records in both parity 1 and 2. Additionally, 25,702 cows had records in both parity 2 and parity 3. Furthermore, 5,496 cows had records in both parity 1 and parity 3. Finally, there were 4,904 cows with records in all 3 parities (parity 1, parity 2, and parity 3).We implemented a repeatability model, as well as 2 bivariate models (first and second lactation, second and third lactation), using a total of 51,415,022 daily milk weights records from 106,033 first parity cows, 89,315 s parity cows and 59,843 third parity cows.These data represented 222 herds from 32 US states (Figure 2).
Repeatability model: where y ijklmn is the TempVar phenotype, Parity i is the fixed effect lactation number with 3 levels, CA j is the fixed effect of calving age with 18 levels, HYS k is the fixed effect of herd calving year season with 2,856 levels, a l is the random effect of cow with 255,191 levels distributed as a ~N 0   where a i is the additive genetic effects for trait i, σ ai 2 is the additive genetic variance of trait i, σ aij is the additive genetic covariance between trait i and j, e i is the residual effect for trait i, σ ei 2 is the residual variance of trait i, and σ eij is the residual covariance between trait i and j.

Correlations between temporal variance of milk yield and other economically relevant traits
Sires with ≥10 daughters with TempVar phenotypes (repeatability model, n = 2,572) were used to calculate approximate genetic correlations between milk yield TempVar and other economically relevant traits using the Calo's method (Calo et al., 1973;Blanchard et al., 1983).PTAs from economically relevant traits evaluated by the Council on Dairy Cattle Breeding (Bowie, MD) were extracted from the April 2023 genetic evaluation.The approximate genetic correlations were calculated as follows: where ˆ, r i j = approximate genetic correlation between traits i and j; ΣRel i and ΣRel j = the sum of reliabilities of traits i and j; Rel i and Rel j = reliabilities of traits i and j; and r i,j = Pearson correlation between PTA for traits i and j.

Modeling lactation curves
The differences observed among lactation models were minimal.The phenotypic means (SD) for LnVar_loess, LnVar_qr05 and LnVar_qr07 in first lactation cows were 3.59 (0.74), 3.59 (0.76) and 3.61 (0.76), respectively.During early lactation (5 to 50 DIM), the data exhibited a relatively weaker fit when using LOESS regression as compared with quantile regression at quantiles of 0.5 or 0.7.This discrepancy may be attributed to the use of a 0.75 span width in LOESS regression.It is possible that reducing the span width could improve the model fit between 5 to 50 DIM.Nonetheless, it should be noted that using a lower span width could potentially result in overfitting when modeling the entire lactation curve, and overfitting is of particular concern when analyzing deviations from the curve to assess consistency or resilience.When TempVar phenotypes were calculated across the entire lactation, correlations between sire PTAs computed from LnVar_loess, LnVar_qr05, and LnVar_qr07 were all > 0.99.In addition, and with the exception of correlations with LnVar_loess in early lactation (0.95), correlations between sire PTAs computed from alternative TempVar definitions were also ≥0.99 within each state of lactation.

Variance components and heritabilities for primiparous cows
Heritability estimates (SE) of TempVar phenotypes for primiparous cows were moderate and ranged from 0.227 (0.011) to 0.237 (0.011) (Table 2), suggesting that consistency or inconsistency of daily milk yield can be altered by genetic selection.Furthermore, the choice of whether lactation curves were modeled as LnVar_loess, LnVar_qr05, or LnVar_qr07 had little impact on estimated variance components or heritability.We found differences in the heritability estimates between lactation stages.Estimated heritabilities of TempVar phenotypes ranged from 0.129 to 0.154 in early lactation, 0.190 to 0.197 in mid lactation, and 0.159 to 0.163 in late lactation (Table 3).Within each TempVar phenotype, however, a similar pattern occurred in which estimated genetic variances were higher in early and mid-lactation and lower in late lactation, whereas residual variances were lowest in mid lactation and highest in early lactation, leading to the differences in heritability estimates described previously.Across all lactation stages, heritability estimates were highest for LnVar_loess and lowest for LnVar_qr07, where it should be noted that fitting a regression curve to the 0.7 quantile seeks to capture the cow's potential milk production in the absence of external perturbations.

Temporal variance of milk yield across lactations
In lactation 1, mean TempVar phenotypes (SD) ranged from 3.59 (0.76) to 3.61 (0.76).In lactation 2, the means ranged from 4.05 (0.75) to 4.07 (0.75), and in lactation 3, the range was 4.22 (0.73) to 4.24 (0.74).Due to a scaling relationship between mean and variance, the mean of TempVar phenotypes (log-transformed variances) increased with lactation number, while the standard deviation (SD) decreased as lactation increased.Repeatability estimates (SE) ranged from 0.331 (0.003) to 0.341 (0.003), indicating TempVar phenotypes appear to be moderately repeatable over time (Table 4), noting that few cows had records in multiple lactations due to the structure of our data set.The assumption of the repeatability animal model is that the genetic correlation between records is equal to 1, indicating that TempVar traits are genetically identical across parities.
To check this assumption, we also analyzed TempVar phenotypes from first, second, and third lactations us- ing a multiple trait model applied to full lactation data.Table 5 shows estimated pairwise genetic correlations between parities for each TempVar phenotype.Genetic correlations (SE) ranged from 0.963 (0.010) to 0.999 (0.003), indicating that milk yield TempVar appears to be same across lactations (Table 5).

Correlations between sire PTAs for temporal variance of milk yield and PTAs for other economically relevant traits
Approximate genetic correlations, based on Calo method among sire PTAs, showed favorable relation-ships between milk yield TempVar phenotypes and several economically relevant traits in the US Net Merit selection index (Table 6).The correlations between TempVar phenotypes and milk yield were 0.57, highlighting that, presumably due to the scaling relationship between mean and variance, TempVar phenotypes tend to increase as milk production increases.Genetic correlations between milk yield TempVar and longevity traits, such as productive life, cow livability and heifer livability, were all negative, indicating that consistent cows (low TempVar) tend to have longer productive lives and lower mortality rates.Interestingly, correlations between milk yield TempVar and health-related traits were all negative.Specifically, for displaced abomasum, ketosis, mastitis, and metritis, we found correlations ranging from −0.21 to −0.41, which indicate that consistent cows tend to have fewer early postpartum health events.Correlations between TempVar phenotypes and somatic cell score were 0.26, which is also favorable, because higher values indicates greater incidence of mastitis.Similar correlations were evident for fertility traits, with correlations of approximately −0.42, indicating that consistent animals tend to have higher PTAs for daughter pregnancy rate, and cow and heifer conception rates.For feed efficiency traits, namely residual feed intake and feed saved, the approximate genetic correlations with milk yield consistency were weak, around 0.05 and −0.05, respectively, which seems to indicate that milk yield TempVar and feed efficiency are independent traits.

Modeling lactation curves
We investigated the use of LnVar of deviations in daily milk weights from their expectations as indicators of consistency in Holstein cows, using 3 different lactation curve models to compute expected yields.We evaluated the heritability and repeatability of milk yield TempVar, as well as the genetic correlations of TempVar phenotypes across methods, lactations, and stages within the lactation, as well as the relationships between milk yield TempVar and other economically relevant traits.Heritability estimates of TempVar phenotypes were similar for all 3 lactation curve models, and sire rankings were essentially identical across methods.This similarity can be partially explained by similarities between the models that were considered, coupled with vigorous outlier detection before fitting the lactation curves.Each of the 3 curve-fitting methods considered in this study are considered as robust to overfitting, which is critical in resilience studies that focus on evaluating the responses of individual cows to observed or unobserved perturbations.Therefore, using an outlier detection function that first decomposes the seasonal trend of the lactation curve and then detects outliers based on the remainder, along with robust methods to model time series data, is relevant.

Heritability estimates across the entire lactation for primiparous cows
Previous studies conducted in the Netherlands and the United States, namely Poppe et al. (2020) and Wang et al. (2022), respectively, have examined the use of deviations in daily milk weights as indicators of resil-ience in dairy cows.Different strategies, such as Woods and Wilmink curves (parametric) and rolling means and medians (non-parametric), have been used to model expected daily milk production.Larger differences (0.198 to 0.244) in heritabilities among lactation fitting methods were observed in Poppe et al. (2020).We did not observe such differences among the methods considered in our study when TempVar phenotypes were calculated using LOESS and quantile regression models applied to full lactation data (0.227 to 0.237).This outcome is likely attributed to a combination of factors.First, the choice of a relatively global span parameter of 0.75 in LOESS regression helped prevent the model from overfitting the data.Furthermore, this could be due to the utilization of polynomial quantile regression, which is a parametric model known for its low risk of being overly flexible (Poppe et al., 2020).Both Poppe et al. (2020) and Wang et al. (2022) reported similar heritability estimates, ranging from 0.198 to 0.250, when analyzing LnVar or log-transformed standard deviations of daily milk weights of Holstein cows.In poultry and pork production, similar approaches have been applied to investigate the LnVar of daily deviations in egg production (Bedere et al., 2022), feed intake (Putz et al., 2019), and body weight (Berghof et al., 2019) as potential resilience indicators.Heritability estimates (standard errors) of 0.10 (0.01) were found for LnVar of daily egg production in purebred laying hens (Bedere et al., 2022), and similar estimates of 0.10 (0.04) were reported for LnVar of weekly body weights in pigs (Berghof et al., 2019).Because of lower heritability estimates and unfavorable genetic correlations with fitness traits in previous studies that used autocorrelation and skewness as resilience indicators (Poppe et al., 2020;Wang et al., 2022;Berghof et al., 2019), we considered only LnVar when computing TempVar phenotypes.

Heritability estimates for different lactation stages
Heritability estimates were highest between 50 to 200 DIM due to smaller residual variances during this period, followed by the period from 201 to 305 DIM.  5. Variance components and heritability estimates (SE) using a multiple-trait model applied to full lactation data of cows with records in lactation 1, 2 and/or 3. Phenotypes representing the temporal variance in daily milk yield were calculated as log-transformed variances of daily deviations from lactation curves predicted with LOESS regression (0.75 span parameter) or polynomial quantile regression (0.5 or 0.7 quantile)

Lactation
No. of cows Heritability estimates were lowest between 5 to 50 DIM due to large residual variances during this period.It is well documented that early lactation is the most challenging time for a dairy cow.Across all 3 TempVar phenotypes, estimates of the additive genetic variance were highest during the period from 5 to 50 DIM.This is an interesting and promising result, as it seems to indicate that genetic differences in milk yield consistency are expressed more fully under challenging conditions, albeit with the complication of greater residual variances.Mulder et al. (2013) previously indicated that greater genetic variation in resilience would be observed when an animal is exposed to environmental challenges.Increased variation in daily milk production during this period may reflect the challenges of decreased voluntary feed intake, coupled with the physiological demands of rapid increases in energy requirements for milk production (White, 2015).Previous authors have described challenges such as negative energy balance (Collard et al., 2000), hyperketonemia (Duffield et al., 2009), and resumption of ovarian cyclicity (Gaillard et al., 2016) during this period, all of which can contribute to an increase in the environmental variance.In contrast, the high heritability estimates for milk yield TempVar between d 51 and 200 reflect low estimates of the estimated environmental variance during this period.Greater heritability estimates in mid lactation, as compared with early and late lactation, were expected due to several factors that tend to reduce the environmental variance.First, metabolic disorders typically occur during the first 60 d postpartum (Pryce et al., 2016), but their incidence diminishes by mid lactation.Second, during the lactation period between 50 to 200 DIM, the cow is able to maximize dry matter intake, maintain a relatively stable body weight, and recover from negative energy balance.Third, cows within their allocated pens have established a social hierarchy, and therefore, milk production, dry matter intake, and rumination time are less affected by establishing new dominance relationships (von Keyserlingk et al., 2008;Schirmann et al., 2011).Heritability estimates for milk yield TempVar during late lactation, from 201 to 305 DIM, were lower than in mid lactation.This difference could be attributed to factors such as an increase in energy requirements for a growing fetus, or by management strategies such as regrouping of cows in late lactation based on nutritional needs or pregnancy status, all of which could add environmental noise and where r 2 is the estimated correlation, and N (2,572) represents the number of sires with ≥10 daughters using a repeatability model applied to full lactation data of cows with records in lactation 1, 2 and/or 3. decrease heritability estimates.However, Poppe et al., 2021 reported heritability estimates for LnVar of daily milk yield deviations of 0.13 during the period from 11 to 110 DIM, 0.12 during the period from 111 to 210 DIM, and 0.15 during the period from 211 to 340 DIM.The differences between studies could reflect differences in definitions of the lactation periods, differences in the number of milk weights required in each period, or differences in grouping strategies of cows in Dutch herds versus US herds.Of special interest, Pearson correlations of TempVar phenotypes across periods of the lactation suggest that TempVar is not a genetically identical trait in different stages (Figure 3).However, TempVar PTAs computed from daily milk yields at 50 to 200 DIM had relatively high correlations with PTAs derived from 5 to 50 DIM (0.66 to 0.67) and those computed from 201 to 305 DIM (0.72 to 0.73), suggesting that selection for milk yield consistency in mid lactation would lead to improvements in generalized resilience throughout the lactation.

Temporal variance phenotypes across lactations
Because first lactation cows are immature and still growing, many traits are genetically or physiologically distinct between primiparous and multiparous cows.Therefore, we decided to estimate heritability, repeatability, and genetic correlations of TempVar phenotypes in first, second, and third parities.Heritability estimates were slightly lower when calculated using the repeatability model, but relatively few cows had TempVar phenotypes for multiple parities, because daily milk yield records from cows that calved before initiation of our research database were unavailable.Genetic correlations among TempVar phenotypes were > 0.95 between parities, indicating that TempVar phenotypes are genetically quite similar throughout the cow's life, and suggesting that we should consider milk yield consistency as the same trait across lactations.Poppe et al. (2020) reported high genetic correlations (around 0.75) between resilience indicators based on daily milk yield deviations and 305-d milk production.These authors estimated the correlations using the multiple across country evaluation procedure (Interbull, 2017).Our estimates between PTAs for milk yield and PTAs for milk yield TempVar were approximately 0.57.Apart from different methods to calculate correlations, one difference between our study and Poppe et al. ( 2020) is that we required each cow to have at least 100 daily milk weights per lactation to be included in the analysis.Another is that we excluded data from automatic milking systems, in which differences in daily milking frequency could contribute to differences in milk yield TempVar.A third is the average size of commercial dairy herds in the US versus those in The Netherlands, because larger farms in the US are more likely to group animals into pens by age, milk yield, pregnancy status, or nutritional needs (Contreras-Govea et al., 2015).Overall, our results suggest that it will be possible to select cows that can achieve both high average milk production and consistent day-to-day performance.

Correlations between temporal variance of milk yield and other economically relevant traits
Approximate genetic correlations between longevity traits, such as productive life and livability, and milk yield TempVar were favorable.Our results suggest that sires whose daughters exhibit consistent daily milk yields tend to have higher PTAs for longevity.Similar correlations with longevity traits were also observed by Elgersma et al., 2018.Fertility traits in the present study also displayed favorable correlations with TempVar phenotypes.Sires whose daughters are more consistent in daily milk yield tend to have above average PTAs for female fertility traits such as daughter pregnancy rate and cow conception rate.
Favorable correlations observed between TempVar phenotypes and health traits in this study indicate that consistent cows may exhibit greater disease resistance and lower somatic cell scores.This is logical, because consistent cows will tend to have fewer disease events, fewer visits to the hospital pen, and fewer management interventions that may cause fluctuations in daily milk yield.The strongest correlations observed in this study were between milk yield TempVar and mastitis, presumably because mastitis is a common disease that causes large decreases in milk production (Liang et al., 2017;Seegers et al., 2003).Overall, lower temporal variance (consistent performance) was associated with superior health, longevity, and fertility.It should be noted that, while we required a minimum of 100 daily milk weights to compute TempVar phenotypes in the present study, we recognize that fragile cows with health or fertility problems may get culled at a higher rate than resilient cows, and we should therefore recognize that improvements in generalized resilience will likely come from a selection index containing EBVs for milk yield TempVar, longevity, and early postpartum health disorders.Correlations of TempVar traits with feed efficiency were near zero, suggesting that selection for lower residual feed intake and greater feed saved will increase farm profitability (Parker Gaddis et al., 2021) with no adverse impacts on resilience.

Future work on resilience indicators
Recent studies of milk yield consistency represent novel contributions to the field by shifting the focus from maximum performance under ideal conditions to sustainable performance under variable conditions.In this study, we defined consistency as the ability of dairy cows to achieve predictable daily milk yield regardless of changes in management and environmental circumstances during the lactation.We contend that milk yield consistency is a potentially valuable indicator of resilience that can be used to improve animal health and welfare while decreasing labor costs associated with management interventions, but resilience is a complex topic that requires deeper study.Resilience is often measured as an animal's response to a known disturbance or perturbation, such as exposure to thermal stress, an infectious disease, or a nutritional challenge (Ravagnolo and Misztal, 2000;Putz et al., 2019;Friggens et al., 2016).Our proposed approach involves using milk yield temporal variance as an indicator of generalized resilience that can be utilized with or without knowledge of the timing and severity of specific management or environmental perturbations.This approach can be expanded to accommodate additional resilience indicators, perhaps in a selection index, if specific perturbations are detected and quantified.Lastly, a novelty of our data set is the availability of daily pen moves and locations.To bridge the gap between consistency and resilience indicators, further research will be carried out to identify calendar dates with greater phenotypic variability (Garcia-Baccino et al., 2021), and we hypothesize that quantifying individual and group responses to perturbations at the pen level will aid in capturing additional resilience indicators in future studies.

CONCLUSIONS
This study aimed to examine the genetics of milk yield consistency within and between lactations of US Holstein cows.Our findings suggest that TempVar is moderately heritable, which may allow selection to focus on cows with smaller fluctuations in daily milk yields throughout lactation than their contemporaries.TempVar phenotypes appear to be robust to the choice of lactation curve models, and genetic rankings seem to be consistent across lactations.Cows displaying superior milk yield consistency tend to be genetically superior for productive life, female fertility, and resistance to early postpartum health disorders relative to their inconsistent contemporaries.Definition of consistency phenotypes and characterization of their genetic basis is an important initial step in developing resilience indicators that will allow selection for consistent performance in unpredictable conditions.Improving resilience will lead to improvements in dairy farm profitability, reduce animal health and welfare risks associated with management and weather disturbances, and improve the social and environmental sustainability of dairy farming.

Figure 1 .
Figure 1.Example of 2 primiparous cows with similar milk production levels from the same herd that are consistent (low temporal variance of daily milk yield) and inconsistent (high temporal variance of daily milk yield) with the following fitted curves: LnVar_loess (blue); LnVar_qr05 (yellow); LnVar_qr07 (green).
Guinan et al.: LACTATION CONSISTENCY IN U.S. HOLSTEINS is the random perma- nent environmental effect distributed as pe ~N 0 2 , , Iσ pe ( ) and e ijklmn is the random residual effect distributed as e ~ N 0 2 , .Iσ e ( )y ijkl = CA i + HYS j + a k + e ijk ,Multiple trait model.where all model terms are as described previously.The assumptions of both bivariate analyses, which were carried out using first and second lactation TempVar phenotypes or second and third lactation TempVar phenotypes, were as follows: Figure 2. Number of cows per state with records for either lactation 1, 2, and/or 3.
Guinan et al.: LACTATION CONSISTENCY IN U.S. HOLSTEINS Table

Figure 3 .
Figure 3. Pearson correlations between sire PTAs for the temporal variance of daily milk yield using a univariate analysis within and across different lactation periods using different lactation curve fitting methods.
Guinan et al.: LACTATION CONSISTENCY IN U.S. HOLSTEINS Guinan et al.: LACTATION CONSISTENCY IN U.S. HOLSTEINS tal of 51,415,022 daily milk weight phenotypes recorded between 5 and 305 d in milk (DIM) remained for our analysis of milk yield TempVar (Table

Table 1 .
Summary statistics of the data sets used for the estimation of genetic parameters Guinan et al.: LACTATION CONSISTENCY IN U.S. HOLSTEINS

Table 2 .
Guinan et al.:LACTATION CONSISTENCY IN U.S. HOLSTEINS Variance components and heritability estimates (SE) 1 for temporal variance (TempVar) of daily milk yield in first parity Holstein cows over the entire lactation.TempVar phenotypes were calculated as log-transformed variances of daily deviations from lactation curves predicted with LOESS regression (0.75 span parameter) or polynomial quantile regression (0.5 or 0.7 quantile)

Table 6 .
Guinan et al.:LACTATION CONSISTENCY IN U.S. HOLSTEINS Correlations (SE) estimated using the Calo's method between sire PTAs for temporal variance (TempVar) of daily milk yield and other economically relevant traits.PTAs for TempVar were obtained for sires with ≥10 daughters using a repeatability model applied to full lactation data of cows with records in lactation 1, 2 and/or 3. PTAs for production, longevity, fertility, health, and efficiency were retrieved from the April 2023 National Genetic Evaluation generated by the Council on Dairy Cattle Breeding.