Genomic-based genetic parameters for milkability traits derived from automatic milking systems in North American Holstein cattle

The number of dairy farms adopting automatic milking systems (AMS) has considerably increased around the world aiming to reduce labor costs, improve cow welfare, increase overall performance, and generate a large amount of daily data, including production, behavior, health, and milk quality records. In this context, this study aimed to (1) estimate genomic-based variance components for milkability traits derived from AMS in North American Holstein cattle based on random regression models; and (2) derive and estimate genetic parameters for novel behavioral indicators based on AMS-derived data. A total of 1,752,713 daily records collected using 36 milking robot stations and 70,958 test-day records from 4,118 genotyped Holstein cows were used in this study. A total of 57,600 SNP remained after quality control. The daily-measured traits evaluated were milk yield (MY, kg), somatic cell score (SCS, score unit), milk electrical conductivity (EC, mS), milking efficiency (ME, kg/min), average milk flow rate (FR, kg/min), maximum milk flow rate (FRM, kg/min), milking time (MT, min), milking failures (MFAIL), and milking refusals (MREF). Variance components and genetic parameters for MY, SCS, ME, FR, FRM, MT, and EC were estimated using the AIREMLF90 software under a random regression model fitting a third-order Legendre orthogonal polynomial. A threshold Bayesian model using the THRGIBBS1F90 software was used for genetically evaluating MFAIL and MREF. The daily heritability estimates across days in milk (DIM) ranged from 0.07 to 0.28 for MY, 0.02 to 0.08 for SCS, 0.38 to 0.49 for EC, 0.45 to 0.56 for ME, 0.43 to 0.52 for FR, 0.47 to 0.58 for FRM, and 0.22 to 0.28 for MT. The estimates of heritability (± SD) for MFAIL and MREF were 0


INTRODUCTION
Automatic milking systems (AMS) are revolutionizing the worldwide dairy cattle industry, by improving animal welfare (Jacobs and Siegford, 2012), reducing the need of human involvement in milking practices, and enabling large-scale automatized data collection at the farm level (Dechow et al., 2020).Another advantage of using AMS is the continuous collection of milking information and related variables, including milk yield, milking speed, milking time (MT), and electrical conductivity (EC; Wethal and Heringstad, 2019).The AMS are equipped with sensors that can measure milk flow rate (FR) and MT, and they automatically store a variety of milking-related variables.Gathering a large amount of daily phenotypic records enables the derivation of novel phenotypes that can be used for genetically improving production efficiency in precision dairy farms.Repeated records can also yield more accurate genetic parameters and pedigree-based or genomically enhanced breeding values (Carlström et al., 2013).In addition, the availability of repeated records enables the use of random regression models (RRM) to estimate genetic parameters for each lactation day (Oliveira et al., 2019a).The RRM provide robust genetic modeling of time-specific effects for longitudinal traits, with key advantages compared with other statistical models, such as (1) prediction of genetic parameters, including genetic correlation among traits, for the time-series within the range of data collection (Schaeffer, 2004); (2) more precise use of data across time; and (3) no need to assume that all longitudinal records collected throughout lactation represent the same trait.
Only more recently, estimates of genetic parameters for milk yield using data derived from AMS have been reported for dairy cattle, as farmers, especially in Europe and North America, started to adopt more AMS (Mulder et al., 2004;Nixon et al., 2009).The genetic correlation between milk yield in AMS and conventional milking systems (e.g., milking parlor) has been demonstrated to be high (Carlström et al., 2014).Moreover, various novel traits derived from AMS have been genetically studied and shown to be of great selective potential and economic relevance (Wethal and Heringstad, 2019;Chen et al., 2020;Dechow et al., 2020;Piwczyński et al., 2021).
Indicators of FR or milking speed are more accurately measured using milking systems with sensors (Siewert et al., 2018), are genetically influenced by the cow temperament (Wethal et al., 2020), and, in turn, could influence milk production, milk quality, reproduction, and health (Chang et al., 2020).Electrical conductivity measures the ability of a solution to conduct an electric current between 2 electrodes, and it is usually measured in milliSiemens (mS).The AMS are equipped with sensors for measuring EC, which is significantly correlated with mastitis (Martin et al., 2021).Thus, AMS-derived traits can be target breeding goals for improving animal health and welfare (Westin et al., 2016).In addition, selection for traits related to the cow behavior during milking in AMS, including MT and milking failures (MFAIL) or milking refusals (MREF), could have a significant effect on the amount of milk harvested in the AMS or a decrease in the number of equipment breaking times (Dechow et al., 2020).
Despite the increased adoption of AMS technology, estimates of genomic-based genetic parameters derived from more sophisticated models such as RRM are still scarce in the literature.Evaluating the genetic background of milkability traits is paramount for the design and refinement of selection indices, especially those focusing on precision dairy systems.To our knowledge, this is the first report of genomic-based genetic parameters for key milkability traits in American dairy cattle.The use of genomic information to predict genetic parameters can provide more reliable estimates, which are essential for designing or refining selection indices in dairy cattle breeding programs.Furthermore, there are no reports of genetic parameters for MFAIL and MREF, which are considered new traits derived from AMS data of great economic importance as more frequent MFAIL would reduce the number of cows that can be milked per AMS station and a large number of MREF indicate an excessive number of visits to the AMS station.Thus, the main objectives of this study were to estimate genomic-based genetic parameters for milkability traits using RRM and novel AMS behavioral indicators in North American Holstein cattle.

Ethics Statement
No Animal Care Committee approval was necessary for the purposes of this study as all the information required was obtained from pre-existing databases.

Data Sets
Data information was obtained from a large commercial dairy cattle farm located in north central Indiana.The dairy currently has 36 milking robot stations (Lely), with capacity to milk 2,200 cows using the AMS.Daily records were used from the AMS management program, recorded for every visit to the AMS unit, and considered for daily milk yield (MY, kg), milking efficiency (ME, kg/min), average milk flow rate (kg/min), maximum milk flow rate (FRM, kg/min), MT (min), milk EC (mS), MFAIL, and MREF.Additionally, testday records of SCC were log-transformed to obtain SCS.A total of 1,752,713 daily records and 70,958 testday records (for SCS), from 4,118 genotyped Holstein cows, were collected from 2017 to 2021.

Genomic Information and Quality Control
A total of 4,118 Holstein cows were genotyped with 4 commercial SNP arrays,containing 11,405 (n = 210 individuals),17,557 (n = 755),18,819 (n = 956),and 30,754 (n = 2,197) SNPs.Genotype imputation was Pedrosa et al.: MILKABILITY TRAITS FROM SENSOR SYSTEMS performed using the FImpute software (Sargolzaei et al., 2014) and a reference population containing 22,090 individuals to converge all SNP panels available to a common set of 72,820 SNPs for all individuals.Quality control (QC) was performed to remove cows and SNPs with call rate lower than 0.90, SNPs located on nonautosomal chromosomes, with minor allele frequency lower than 0.05, and extreme deviation from Hardy-Weinberg equilibrium (P < 10 −8 ).The genomic QC was done using the BLUPF90 family programs (Misztal et al., 2018).After the QC, 57,600 SNPs and 4,118 animals remained for further analyses.

Trait Definition
Continuous Traits.Daily MY was defined as the total amount of milk produced in 24 h, where milk records from multiple milkings throughout the day were summed up to compute MY.The SCC was log-transformed to obtain SCS based on the formula (Piwczyński et al., 2021): tivity is a measure of the resistance of a specific substance to an electric current.During milking, EC measures the increased blood capillary permeability resulted by mastitis, which affects the ion concentrations in milk, enabling the measurement of EC in milliSiemens.
Electrial conductivity was measured individually, and the values were obtained from the daily AMS report from the Lely system.Milking efficiency was calculated by dividing the average daily milk production by the daily MT.The FR was calculated as an average of milk in kilograms per minute of MT, obtained from the AMS milking report.The FRM was calculated based on the daily average of the 2 highest flows obtained within each milking activity.Finally, daily MT was calculated as the total amount of time considering all milking events per cow within 24 h.Categorical Traits.Milking failure is a trait that represents the total number of failed milkings on a day.
The record of an animal with one or more failures within a day was set as 1, while records from cows that did not have any failures on a given day were set as 0. This was done because less than 0.7% of the MFAIL occurred 2 or more times in a day for a particular cow.The main reasons for MFAIL are (1) teats not found by the robot; (2) connection time, where teats were detected but the milker robot was unable to connect one of the teats; (3) the AMS stopped because a cow left the AMS unit before the gate was closed; (4) dead MT, in which teat cups were successfully connected but there was no milk flow from one (or more) of the quarters.
Milking refusals are the total number of refused milking by the AMS unit on a day.Milking refusals indicates a consistent flow through the AMS unit and therefore, should not be interpreted as a negative aspect of the cow.Constant refusals occur by an adequate adaptation of the cows to AMS, which makes them want to visit the system more often than necessary.A cow with one or more refusals within a day was set as 1 in the final data set, and females who did not present refusals on a day were set as 0. A certain number of refusals per cow during the lactation season are normally caused by: animal adaptation to the AMS; high feed dispensure speed, leading to leftovers of concentrate in the AMS; and, low feed amount provided at the AMS feeder.
Phenotypic data QC was performed using the R software (R Core Team, v 4.1.2),with a threshold of 3.5 standard deviations from the mean to exclude potential outliers for the continuous traits.The descriptive statistics for each trait in the final data set (after QC) are presented in Table 1.
Statistical Models and Genetic Parameter Estimation.Variance components and genetic parameters for MY, SCS, ME, FR, FRM, MT, and EC were estimated using the AIREMLF90 software (Misztal et al., 2018).The RRM fitting a third-order Legendre orthogonal polynomial (LP) were used to describe the covariance structure of the data as a function of DIM from 5 to 305 d postpartum (de Boor, 1980;Meyer, Pedrosa et  ( ) is a regression function ac- cording to DIM j using a third-order LP.The residual effect (e ijkl ) was assumed to be homogeneous throughout lactation and the random effects were assumed to be normally distributed with zero mean.A heterogeneous residual effect throughout lactation was previously tested but resulted in no convergence (results not shown).
A threshold Bayesian model using the Gibbs sampler and Markov Chain Monte Carlo algorithm was considered for MFAIL and MREF, as follows: where y corresponds to the phenotypic records, μ is the overall mean; YS is the concatenation of calving year and season (March to May, June to August, September to November, December to February) of calving, fitted as a random effect; LN is the fixed effect of the lactation number (1, 2, or 3+); DIM is the linear covariate fixed effect of days in milk (5-305 d); pe and a are the permanent environmental and the additive genetic effects, respectively; e is the residual effect.The variance components and genetic parameters for MFAIL and MREF and genetic correlations between these 2 traits with all other traits included in the study were estimated using the THRGIBBS1F90 software (Misztal et al., 2018).Convergence achievement was verified using Geweke (1992) and Heidelberger and Welch (1983), and visual criteria using the "BOA -Bayesian Output Analysis" package (Smith, 2007) implemented in the R software (R Core Team, v 4.1.2).No pedigree information was available, and therefore, all the analyses were based on a genomic BLUP approach.

Phenotypic Information and Variance Components
Phenotypic averages measured in the first, second, and third lactation according to DIM for traits with daily records are presented in Figure 1.The pattern of the phenotypic curves was similar when comparing the second and third lactations for all traits, except for EC and MT where only in the first 35 DIM the curve of the first lactation differed from the others.In addition, for MY, ME, FR, and FRM the curve in the first lactation was substantially lower compared with the second and third lactation curves.The opposite was observed for SCC and EC, where for the first lactation, higher values were observed in the beginning of DIM, dropping as the lactation progresses and stabilizing in the final part of the lactation.The additive genetic variance and permanent environment variance according to DIM for MY, SCS, EC, ME, FR, FRM, and MT are shown in Figure 2. Higher levels of additive genetic variance compared with the permanent environment variance were observed for EC, ME, FR, and FRM.Furthermore, for SCS and MT, it was observed that at the beginning of the lactation the permanent environment variances were higher, with a reduction along the lactation.Contrariwise, for these 2 traits, the additive genetic variances were lower in early lactation, with a gradual increase as DIM increased, which resulted in an inversion of the curves positioning at certain times of the lactation.

Heritability Estimates
The heritability estimates for the traits analyzed by RRM showed low variation throughout the lactation, except for MY in which the heritabilities increased along the DIM (Figure 3).The daily heritability estimates ranged from 0.07 to 0.28 for MY, 0.02 to 0.08 for SCS, 0.38 to 0.49 for EC, 0.45 to 0.56 for ME, 0.43 to 0.52 for FR, 0.47 to 0.58 for FRM, and 0.22 to 0.28 for MT.The shape of the heritability curves for EC, ME, FR, and FRM was similar, with lower values around the first DIM, increasing rapidly in subsequent days reaching the highest heritability values around the lactation peak, and stabilizing right after across the DIM.For SCS and MT, the shape of the heritability curve presented lower values until the middle of the lactation with a smoother gradually increasing wave from the middle until the end of lactation.Finally, the heritability estimates for MFAIL and MREF, evaluated based on a Bayesian threshold model are presented in Table 2.The heritability (± SD) estimates for MFAIL and MREF were 0.02 ± 0.01 and 0.09 ± 0.01, respectively.

Genetic Correlations
The genetic correlation estimates across the DIM per trait for the longitudinal traits is shown in Table 3.Most of traits (FR, FRM, ME, and MT) were highly genetically correlated (≥0.80) along the whole lactation for each trait, indicating that selection based on lactation stage specific records would result in improvement throughout the whole lactation.For EC, even though a high genetic correlation was observed in the most part of lactation, a stronger decrease trend in the genetic correlation was evidenced, especially on more distant DIM, as in the case of DIM 5 versus DIM 305.Moreover, major differences between DIM were observed only for MY and SCS, where the genetic correlations between distant DIM (5 d vs. ≥185 d) were low.This result might be related to the intrinsic variability existing along DIM for MY and SCS but could also indicate that a different set of genes are involved in the phenotypic expression of those traits across lactation stages.
The genetic correlation between traits analyzed using the RRM approach is presented in a heatmap in Figure 4, showing the genetic relationship among traits across DIM.Slight differences were observed across DIM within the same traits.For most correlated traits lower genetic correlation was observed at the beginning of the lactation, gradually increasing positively along DIM.However, different pattern was observed for MY-SCS and MY-EC, in which higher and positive genetic correlations were detected in the first third of lactation and decreasing later until achieving negative values.The genetic correlation between SCS-EC was moderate and ranged from 0.25 to 0.37 with a slight fluctuation across DIM.Interestingly, the genetic correlation between EC and the other continuous traits (except SCS) were close to zero with a small oscillation across low values.Moderate to high positive genetic correlations were observed between SCS-ME, SCS-FR, and SCS-FRM, mostly ranging from 0.50 to 0.70, indicating that milking speed is genetically related to SCS across different stages of lactation.Also, moderate to high and  observed between FR-FRM.Low standard errors for all genetic correlations were observed, ranging from 0.01 to 0.03 (Figure 4).The genetic correlation between MFAIL and MREF evaluated based on a Bayesian threshold model, as well as their genetic correlation with the other studied traits are presented in Table 2.The genetic correlation (± SD) between MFAIL and MREF was 0.25 ± 0.02, indicating a low to moderate genetic relationship between these 2 traits.For MFAIL and the other studied traits, most genetic correlations were low to moderate (−0.05, 0.18, 0.14, −0.03, and −0.06 between MFAIL and MY, SCS, EC, ME, and MT, respectively).However, moderate and negative genetic correlations were observed between MFAIL and FR (−0.58 ± 0.02) and between MFAIL and FRM (−0.56 ± 0.02), indicating that cows with more MFAIL tend to have lower FR.The genetic correlations between MREF and the other analyzed traits were low (Table 2).Positive and low associations were observed for MREF-MY (0.11 ± 0.01), MREF-ME (0.11 ± 0.01), MREF-FR (0.10 ± 0.02), MREF-FRM (0.12 ± 0.02), and MREF-MT (0.14 ± 0.02).Negative, but still low, genetic correlations were observed between MREF-SCS (−0.12 ± 0.02) and MREF-EC (−0.10 ± 0.01), demonstrating a slight genetic tendency of a relationship between higher numbers of MREF resulting in a decrease in the number of SCS and EC, and vice versa.

DISCUSSION
Sensor-based systems in dairy farming have enabled the accurate collection of data for numerous variables and traits, including herd-management variables, feed efficiency indicators, animal behavior, animal health, and the exact time assessment in the machine (Wethal and Heringstad, 2019).In this scenario, the use of AMS has been a key technological advancement in the dairy industry in the past few decades (Cogato et al., 2021).Milking robotic systems are considered one of the new technologies increasing the progress in the phenomics era, and, for instance, great contribution can be made not only on dairy farming, as well as in the animal breeding side.As an example, milkability traits related to efficiency as ME, FR, FRM, and MT can be precisely measured under AMS, and consequently, only properly genetic estimated in herds belonging to farms containing milking robotics (Carlström et al., 2013).In addition, indicator traits of cow behavior such as MFAIL and MREF, measured in AMS, have become increasingly important due to their relationship with animal welfare and ME.Selecting for favorable cow behavior in AMS could enhance ME and improve milk productivity (Dechow et al., 2020).Finally, sensorbased systems, especially those monitoring production, ME, and behavior, provide extensive information about individual animals at each milking event.Thus, enabling the application of more complex models to longitudinal measures to be used for a more accurate estimation of the dairy information collected.The advantages of using RRM to investigate longitudinal traits in dairy cattle have been discussed in the literature (Henderson Jr., 1982;Schaeffer, 1994;Jamrozik et al., 1997;Oliveira et al., 2019a).From a genetic perspective, among the main benefits of the application of RRM is the increase on accuracy of breeding values as well as improving estimates of genetic parameters, especially because of the possibility of fitting in the model random genetic and random environmental effects along time (Veerkamp et al., 2001;Oliveira et al., 2019a).More recently, the use of RRM with the addition of genomic information for the genetic evaluation of longitudinal traits was recommended as a feasible alternative for a more precise genetic estimation based on the complete pattern of the production curve compared with the parent average (PA) only, when genomic information is not available (Koivula et al., 2015;Oliveira et al., 2019b).Thus, using the same RRM approach with the addition of genomic informa-tion to longitudinal data measured in AMS should provide more accurate estimates of genetic parameters that allow point-to-point estimation along the lactation curve, which has not yet been explored in AMS.

Variance Components
Almost constant additive genetic variance was observed across DIM for milking speed-related traits (ME, FR, and FRM) and MT, where similar shape was observed over the entire curve with a slight oscillation along the DIM.In addition, it was noted that the environmental effects were less influential for EC, ME, FR, and FRM, even during the different stages of lactation, and distinct environmental occurrences along the year.It is also relevant to highlight that some studies have demonstrated that the use of RRM could better capture the additive genetic variance compared with repeatability or 305-d adjusted models (Wang et al., 2009;Kramer et al., 2013).In this context, the use of RRM enabled the identification of the pattern of genetic variance for the milking speed-related traits, highlighting the consistent additive genetic variance for ME, FR, and FRM throughout the lactation.
For SCS and MT, the permanent environment variances were higher in the beginning of the lactation, with a reduction along DIM associated with a constant increase in the additive genetic variance, especially for SCS.Higher values of SCC are constantly observed in the final third of the lactation for different species as goats (Tedde et al., 2019), buffaloes (Sodhi et al., 2021), and cattle (Graesbøll et al., 2016), mainly because of physiological differences across lactation stages, and variations in milk cell count (Souza et al., 2012).Lactation is associated with several biological changes, where fat and glucose synthesis greatly increase along DIM contributing to the variations at the end of the SCS lactation curve (Wankhade et al., 2017;Lázaro et al., 2021).With this, it becomes evident that the permanent environment effects play a crucial role in the final part of lactation for SCS, with an additional variation in the additive genetic component, which consequently influence its variation throughout the lactation of dairy species, including dairy cattle.Although all the data used in this study was obtained from a single farm, it is one of the largest herds using AMS in the United States with over 5,000 lactating cows.In addition, the farm has agreements with AI companies and semen from a wide variety of bulls have been used over time and they also own dams (and close relatives) of widely used AI bulls.Therefore, it is reasonable to assume that the variance components obtained in this study are a good representation of the North American Holstein population parameters.However, future studies using larger data sets are suggested.

Heritability Estimates
Estimates of heritabilities using RRM of longitudinal traits enabled the estimation of time-dependent components at different stages of lactation, which is a great advantage in terms of interpretation of the estimates because some heritability variation can be identified across DIM.Even though the heritability estimates found in our study revealed low variation throughout the lactation for most traits analyzed in RRM, the curve of heritability for MY along the lactation have demonstrated substantial differences in the estimates ranging from 0.07 to 0.28 (Figure 3).Similar variation for MY was reported in other studies (Brito et al., 2017;Bohlouli et al., 2019;Atashi et al., 2021).However, most RRM studies used test-day information, which might indirectly contain a bias associated with the nondaily data collection system.One of the great advantages of our study is that daily production data are available for a more accurate RRM analysis, because with the AMS real-time information is available from all adjusted points in the regression throughout lactation.
The daily heritabilities for SCS ranged from 0.02 to 0.08 with a gradual increase after the middle of the lactation.The heritability curve pattern for SCS was nearly equal to those observed for the additive genetic variance, with the highest estimates appearing in the final part of the lactation, indicating that higher values of heritability are well connected to greater values of the additive genetic variance presented in advanced DIM.Our results reflect the usually low heritability of SCS, which is in accordance with several studies (Koeck et al., 2012;Negussie et al., 2013;Kheirabadi, 2018;Soumri et al., 2020), in which the reported estimates ranged from 0.02 to 0.12 when using similar methodologies.In our case, a factor that may have influenced the low heritability values for SCS is the strict sanitary control adopted by the property where the data were collected (i.e., low environmental variability).As shown in Figure 1, the mean phenotypic values of SCC ranged from 80 to 220 (× 1,000 cells/mL), which demonstrates how well the SCC was controlled in this farm.
The SCC is an important and widely used indicator of mastitis diagnostics, especially for the detection of subclinical mastitis.However, EC, which is another mastitis indicator trait has been used mainly in farms under AMS with a successful identification of mastitis occurrence, presenting highly estimates of heritability compared with SCS (Norberg et al., 2004b;Wethal et al., 2020).The use of EC as a strategy for detecting mastitis was proposed in the 1980s, with an increase in its usage in the following decades (Fernando et al., 1982;Nielen et al., 1992;Norberg et al., 2004a).Nevertheless, the great leap in the use of this tool on a large scale for dairy occurred after the widespread implementation of AMS, which allowed real-time inference of EC in each milking event.The heritabilities across DIM for EC ranged from 0.38 to 0.49 with higher values around the peak of lactation.These results are close to those reported in the literature for EC, in which moderate to high heritability magnitudes are described for dairy cattle, normally ranging from 0.23 to 0.38 (Norberg et al., 2004a;Egger-Danner et al., 2015;Wethal et al., 2020).The slightly higher values presented here in comparison to the literature estimates can be justified by the different statistical approach used, where all previous studies evaluating EC used a repeatability model versus the RRM fitted in this study.RRM usually presents a lower error associated in comparison to repeatability models.This is because a repeatability model does not consider that the correlation across phenotypes recorded at closer times is higher than the correlation between more distant stages, affecting the genetic and environmental variances over time (Van Der Werf et al., 1998;Oliveira et al., 2019a).Due to the ease measurement, in addition to the high heritability estimates, it can be recommended EC as a genetic selection criteria of dairy cattle breeding programs to control mastitis, particularly in AMS herds.
Estimates of heritabilities for milking speed-related traits ranged from 0.45 to 0.56 for ME, 0.43 to 0.52 for FR, and 0.47 to 0.58 for FRM.The shape of the heritability curves for ME, FR, and FRM presented lower values in the beginning of lactation, increasing in subsequent DIM, reaching the highest heritability values near the lactation peak.Adaptation and milk production effectiveness in AMS are well dependent on the cow's ME, and thus, the relation of milking speed and MT is crucial for a profitable process in AMS.Consequently, it is important that milking speed-related traits present high heritability estimates for higher selection responses in breeding programs aiming to breed for more efficient and adapted animals for this type of milking system.Moderate to high heritability estimates for ME measured in AMS were previously reported in Holstein and Norwegian Red cattle in Europe, ranging from 0.27 to 0.50, where higher estimates were obtained using an RRM compared with other studies that have used repeatability models (Carlström et al., 2013;Bakke and Heringstad, 2015;Vosman et al., 2018;Wethal and Heringstad, 2019).Also, high heritability estimates were reported for FR and FRM in Holstein and Swedish Red cows raised in the United States and Europe, varying from 0.37 to 0.55, confirming the great magnitude of heritability estimates for milking speedrelated traits measured in sensor-based systems (Pretto et al., 2014;Carlström et al., 2016;Dechow et al., 2020).
For MT, the shape of the heritability curve presented lower values until the middle of DIM with a gradually rising curve from the middle until the final part of the lactation, ranging from 0.22 to 0.28.Milking time usually followed the same pattern of milk production over the years, increasing in almost the same proportion as milk yield increased in the last decades, mainly because MT was only recently considered as a selection criterion (Miglior et al., 2017).However, the production efficiency in AMS is related to the increase in milk production in a shorter period.Still, milk production was prioritized in the objectives of selection of dairy cattle in the past years, which, in turn, did not contribute to the increase in production efficiency because selection for MT was not previously considered (Cogato et al., 2021).Therefore, MT can be considered as an economically relevant performance indicator, also clearly related to management and animal behavior, especially in dairy sensor-based systems (Lyons et al., 2014).Other studies have proved the selective potential for this trait, showing heritabilities varying from 0.13 to 0.42 in Holstein, Jersey, and Norwegian Red cattle (Carlström et al., 2014;Wethal and Heringstad, 2019;Løvendahl and Buitenhuis, 2022).These finding corroborates the concept that MT should be incorporated into genetic breeding programs to increase the efficiency of dairy herds.
Finally, the heritability estimates (± SD) for MFAIL and MREF, evaluated based on a threshold Bayesian model were 0.02 ± 0.01 and 0.09 ± 0.01, respectively.Despite the low heritability, MFAIL and MREF are under genetic control.The low heritability estimates presented for these traits may be related to the low frequency of occurrences of MFAIL and MREF in the investigated herd, as shown in Table 1.It is expected a learnability (cows learning to use the AMS over time) associated with refusals, but especially with failures in dairy cattle AMS (Siewert et al., 2019), which indicates that cows with greater lactation number tend to have a lower number of failures in later lactations.As in this study the most advanced lactations were evaluated together with the first lactation, a reduction in the number of failures may have resulted in a lower proportion of MFAIL and MREF across lactations.Other traits with similar pattern of measurement as kick-off during milking and incomplete milkings have also presented low heritabilities, as seen in a study in Holstein (0.02 ± 0.03) for incomplete milkings (Carlström et al., 2016), and Norwegian Red cattle (0.06 ± 0.01 and 0.01 ± 0.01; Wethal and Heringstad, 2019), for kick-off and incomplete milkings, respectively.Both MFAIL and MREF are new behavioral indicators for AMS in dairy cattle and can be further explored in future studies, to assess the possibility of including one or both of them in genetic selection schemes.

Genetic Correlations
High and positive genetic correlations along DIM within the same trait were observed for FR, FRM, ME, and MT, with estimates ≥0.80 across the whole lactation.These results demonstrate an intimate genetic relation for the milking speed-related traits along the entire lactation.In this sense, animals with a high genetic merit for milking speed traits in the beginning of lactation tend to also demonstrate high breeding values for those traits until the end of lactation, which facilitates the selection process.High genetic correlations for FR (0.95-0.99) and MT (0.95-0.99) using RRM in a Swedish Holstein cattle population were reported in the literature (Carlström et al., 2014), reinforcing that we would select most of the same cows regardless of the lactation stage.The EC across DIM also presented a high genetic correlation in most parts of the lactation with a decrease in the genetic correlation between more distant DIM (e.g., 5 vs. 305 DIM).This indicates that there are differences between the gene actions at different stages of lactation for EC.Thus, different genetic interpretation must be considered at different stages, especially between early and late lactation.Finally, the largest differences between DIM within trait were observed for MY and SCS, where the genetic correlation between the far apart DIM (5 vs. ≥185 d) were low, demonstrating a polygenic nature of those traits across DIM.A similar pattern was observed in previous studies in Holstein cattle milked in different systems as parlor, either for MY or SCS (Bohlouli et al., 2019;Salimiyekta et al., 2021).
Considering the genetic correlation between different longitudinal traits, small differences were observed across DIM within the same pair of traits.Similar pattern was previously identified for several milk related traits in European Holstein populations using RRM analyses, suggesting that a reranking of animals is not expected across DIM (Calus and Veerkamp, 2003;Paiva et al., 2022).However, RRM are important for the estimation of genetic parameters because negative genetic correlations between MY-SCS are usually presented in studies using repeatability models (Abdalla et al., 2016;Bobbo et al., 2019).As observed here and previously demonstrated in other experiments, the additive genetic correlations between MY-SCS vary throughout lactation, representing magnitude differences not captured by other analytical models (Yamazaki et al., 2013;Atashi and Hostens, 2021).
The genetic correlation between mastitis indicators, SCS-EC, was moderate, ranging from 0.25 to 0.37, slightly fluctuating across the DIM.Unfortunately, clinical mastitis information was not available in our data set, and therefore, it was not possible to include this trait in the set of analyses performed.Studies in the literature have already reported moderate to high genetic relationships between mastitis and SCS (0.92-0.99), as well as mastitis and EC (0.65 to 0.80) in dairy cattle (Norberg et al., 2004a;Urioste et al., 2012;Martin et al., 2018;Kirsanova et al., 2019).Nevertheless, few studies have demonstrated the genetic correlation between SCS and EC, which is important because, depending on the situation as in farms with no AMS, only one of the traits can be recommended as an indirect indicator for the selection of animals more resistant to bovine mastitis.Wethal et al. (2020), working with 77 herds of Norwegian dairy breeds, have also indicated moderate genetic correlation between SCS-EC from 0.34 to 0.37, similar to the estimates obtained here.Considering the moderate genetic correlation between SCS-EC, it is not possible to directly recommend that one trait be replaced by the other as an indirect indicator of selection for mastitis in dairy cattle.However, as both traits have the same objective of being indirect indicators of mastitis, it can be recommended EC for cattle milked in AMS due to the abundant availability of data collected at each milking event in AMS, and also because of its higher heritability compared with SCS, as shown in our results.
Moderate to high positive genetic correlations were observed between SCS-ME, SCS-FR, and SCS-FRM, mostly ranging from 0.50 to 0.70, indicating that milking speed may be genetically related to SCS across lactation stages.Cows with higher milking speed tend to be more productive causing higher physiological stress, which in turn, decreases their immunity leading to more SCS in their milk content (Piwczyński et al., 2020).Another physiological association might be related to teat canal conformation, where more open sphincters could allow the entrance of a larger number of bacteria, indirectly resulting in the increase of SCS.For this reason, it is necessary that genetic selection for an intensive increase in dairy production be accompanied by selection aimed at controlling SCC, in addition to a strict sanitary control.Also, moderate to high genetic correlations were observed between MT and other traits such as SCS, ME, FR, and FRM, but in this case negative values, ranging from −0.86 to −0.48, with higher values in the final lactation stage.An inverse relationship between MT and milking speed-related traits was already expected, because selecting for more efficient animals tends to produce herds with higher milking speed and at the same time produce more milk in a shorter time.Another study with US Holstein cattle identified similar results, showing genetic correlations between MT and FR varying from −0.85 to −0.34, highlighting the importance of increasing AMS efficiency without compromising cow health and well-being (Dechow et al., 2020).Healthy udders should be able to allow increased milking speed without damaging the conformation of the mammary system throughout lactations over time, demonstrating ME as well as adaptation to sensor-based systems.
In this context, the stronger genetic correlations among all longitudinal traits were noted between traits related to ME and milking speed, where positive values from 0.94 to 0.99 were obtained.Tremblay et al. (2016) addressed the impact of FR in dairy herds raised in AMS, indicating this trait as one of the most important factors affecting milk production.Animals identified with high milking speed occupies the AMS unit for a shorter period, promoting a more efficient use of the automatic milking units, thus increasing the farm profitability (Piwczyński et al., 2020).Selection for traits related to milking speed can also indirectly affect cows' longevity, increasing lifetime yield and life index, which also have a great influence on long-term efficiency (Sewalem et al., 2010).It is evident that genetic selection for milkability traits is essential for improving efficiency, adaptation, and welfare of animals raised in precision farms.In this context, the selection must go far beyond the simple increase in production, but rather have a direct relationship with the best adaptation to robotic milking through selection for a variety of traits collected by this type of milking system.
A low genetic correlation (± SD) between MFAIL-MREF was observed (0.25 ± 0.02), indicating a low genetic influence of MFAIL on MREF and vice versa.The low genetic relationship between these traits makes sense as MFAIL are related to undesirable behavior of animals within AMS.Milking refusals may be indicative of the good adaptation of females to AMS, and, therefore, can sometimes be linked to good behavior when it is not exaggerated.Both traits, MFAIL and MREF, can be considered new indicators of behavior in AMS, because, as far as we know, there are no studies on the estimates of genetic parameters for these traits.Other characteristics related to cows' behavior in AMS have already been reported, such as box time, incomplete milking or kick-off machine (Carlström et al., 2016;Wethal and Heringstad, 2019;Dechow et al., 2020), which added to the findings reported here for MFAIL and MREF and provide new knowledge associated with dairy cows' behavior in AMS.
The low genetic correlations between MFAIL and other traits such as MY (−0.05),SCS (0.18), EC (0.14), ME (−0.03), and MT (−0.06), demonstrate that there is a small genetic relationship between the number of failures and milkability traits.However, moderately negative genetic correlations were observed between MFAIL-FR (−0.58 ± 0.02) and MFAIL-FRM (−0.56 ± 0.02), indicating an inverse and favorable genetic relationship between MFAIL and FR/FRM, where cows with a lower incidence of MFAIL are expected to have higher breeding values for FR.By definition, one of the causes of MFAIL is dead MT, where teat cups are effectively connected to the udder system, but a milk flow problem occurring in at least one of the quarters cause a forced break in AMS.The presence of a nonlactating quarter or a low milk FR were previously associated with a decrease in MY, longer MT spent in the AMS unit, and occurrence of failures in the AMS in US Holstein herds (Wieland et al., 2017).
Finally, low positive genetic correlations were obtained for MREF-MY (0.11 ± 0.01), MREF-ME (0.11 ± 0.01), MREF-FR (0.10 ± 0.02), MREF-FRM (0.12 ± 0.02), and MREF-MT (0.14 ± 0.02).France et al. (2022) reported that high-yielding cows visited an AMS unit more often, not only because of the greater milk production, but also for the high demand of feed of-fered in the AMS units.This behavior can explain in part that high-yielding cows are those with a genetic trend to slight present higher numbers of MREF.Contrary, negative low genetic relationships were observed between MREF-SCS (−0.12 ± 0.02) and MREF-EC (−0.10 ± 0.01), demonstrating a small genetic tendency that animals with higher breeding values for MREF present negative breeding values for SCS and EC.This genetic relationship follows the propensity that more productive animals, on average, present negative values for SCS and EC, mainly in the final stage of lactation (Piwczyński et al., 2021).A positive aspect of the low genetic correlation between MREF and the other evaluated traits is its genetic independence from the other milkability traits.Thus, intense genetic selection for MREF can be performed without indirectly affecting other traits of economic interest related to AMS.

Potential Implications and Limitations
Novel estimates of genetic parameters for traits derived from AMS in dairy cattle were reported, confirming previous estimates of heritability and genetic correlations, but also presenting new important findings for North American Holstein cattle.The use of RRM for analyzing traits measured in AMS showed great potential because the high number of phenotypic information provided by AMS allows the implementation of robust models for longitudinal traits, as seen in this type of milking system.Additionally, the use of genomic information for the prediction of estimates of genetic parameters in sensor-based systems proved to be effective.With this, genetic models that jointly consider the application of RRM and genomic information are fully feasible, being able to improve the reliability of the estimates of parameters that are essential for AMS dairy breeding programs.Other studies evaluating the production and control of the mammary system per quarter in AMS can be carried out, verifying the potential productive differences existing not only between the quarters, but also between the front and rear compartments of the mammary gland.Such studies can be useful for a better targeting of the selection of animals that present better adaptation to AMS, in addition to allowing a better understanding of the genetic relationships with the physiology of milk production in the complex of the mammary system.In the end, it is important to highlight that many estimates obtained here are limited to robotic milking systems, which, despite showing a great growth trend in their use in the global dairy scenario, are still present in a limited number of herds.Estimates of genetic parameters may vary across populations, methods of analysis, but also Pedrosa et al.: MILKABILITY TRAITS FROM SENSOR SYSTEMS according to the breeding system adopted, which leads to an interpretation of the results according to different dairy realities.

CONCLUSIONS
The heritability estimates for longitudinal AMS traits using RRM revealed low variation across DIM for most traits evaluated.The novel AMS behavior indicators MFAIL and MREF presented low estimates of heritabilities, with a higher selective potential for MREF than for MFAIL.Slight differences in the genetic correlations across DIM were observed within the same longitudinal trait.For most genetic correlations among traits, lower estimates were observed at the beginning of lactation, gradually increasing along DIM.Moderate to high positive genetic correlations were observed between SCS-ME, SCS-FR, and SCS-FRM, indicating a potential genetic connection between milking speed traits and SCS.Remarkable high and negative genetic correlations were observed between MFAIL-FR and MFAIL-FRM, indicating that cows with a genetic tendency to have more MFAIL are those with lower milk FR.Genetic selection for milkability traits can directly benefit farmers who seek to improve milking efficiency of herds managed in AMS.

Figure 1 .Figure 2 .
Figure 1.Phenotypic average estimated in the first, second, and third lactation for milk yield (MY), SCC, electrical conductivity (EC), milking efficiency (ME), milk flow rate (FR), maximum milk flow rate (FRM), and milking time (MT) in North American Holstein cattle according to days in lactation, obtained based on random regression models.
Figure 3. Heritability (h 2 ) estimated for daily milk yield (MY), SCS, electrical conductivity (EC), milking efficiency (ME), milk flow rate (FR), maximum milk flow rate (FRM), and milking time (MT) in North American Holstein cattle according to DIM, obtained based on random regression models.

Figure 4 .
Figure 4. Genetic correlation (rg) estimated for milk yield (MY), SCS, electrical conductivity (EC), milking efficiency (ME), milk flow rate (FR), maximum milk flow rate (FRM), and milking time (MT) in North American Holstein cattle according to DIM, obtained based on random regression models.

Table 1 .
al.: MILKABILITY TRAITS FROM SENSOR SYSTEMS Summary statistics for milk yield, SCC, electrical conductivity, milking efficiency, milk flow rate, milk flow rate maximum, milking time, milking failures, and milking refusals in North American Holstein cattle

Table 2 .
Estimates of h 2 and genetic correlation (rg) ± SD for daily milking failures and daily milking refusals, as well as their genetic correlation with other variables 1 in North American Holstein cattle measured in automatic milking systems negative genetic correlations were observed between MT and other traits such as SCS, ME, FR, and FRM, ranging from −0.86 to −0.48, with higher values in the last third of lactation.Finally, stronger genetic correla-tions were observed between all traits related to ME and milking speed with estimates ranging from 0.94 to 0.99.The genetic correlations across DIM for ME-FR and ME-FRM were similar, with a small variation only Pedrosa et al.: MILKABILITY TRAITS FROM SENSOR SYSTEMS

Table 3 .
Estimates of genetic correlation within trait across DIM for daily milk yield (MY), SCS, electrical conductivity (EC), milking efficiency (ME), milk flow rate (FR), milk flow rate maximum (FRM), and milking time (MT) in North American Holstein cattle measured in automatic milking systems

Table 3 (
Pedrosa et al.:MILKABILITY TRAITS FROM SENSOR SYSTEMS Continued).Estimates of genetic correlation within trait across DIM for daily milk yield (MY), SCS, electrical conductivity (EC), milking efficiency (ME), milk flow rate (FR), milk flow rate maximum (FRM), and milking time (MT) in North American Holstein cattle measured in automatic milking systems Pedrosa et al.: MILKABILITY TRAITS FROM SENSOR SYSTEMS