Genetic and phenotypic variation and consistency of cow’s preference and circadian use of robotic milking units

Many dairy herds use automatic milking stations (AMS) and in large herds, cows usually have access to two or more AMS, and consequently must make a choice between them whenever they go for milking. For this study, we assume that the degree of using the same milking box is an expression of preference. We used a total of 4,665 cows from a Holstein and a Jersey cohort covering 12 commercial dairy herds in Denmark to estimate heritability, repeatability and individual level correlations between preference choice score (PCS), Time_profile (preference for specific clock hours) and associated milking and behaviour traits. The average PCS was 0.43 (Holstein) and 0.41 (Jersey), respectively. Both PCS and Time_profile had low heritability (0.07/0.11 Holstein, 0.13/0.04 Jersey) and moderate repeatability (0.47/0.40 Holstein and 0.50/0.42 Jersey) and were weakly correlated to each other (Holstein 0.18 and Jersey 0.17) and weakly correlated to milk yield (r i range 0.0 to -0.10).


Introduction
Dairy cows have been subject to intensive genetic selection for production traits for decades, and many cows are living in housed environments interacting with the confinements of housing, and with various automated installations, including robotic milking machines. In these modern dairy herds cows are making a number of choices every day. For example, at milking cows form an order of entry to milking parlours, and through their use show preference to their favourite side of the milking parlour (Hopster et al. 1998). When accessing automated milking stations (AMS, milking robots), they may acquire a routine with preference for a specific one, or be more flexible. If such routines are long lasting, it indicates a degree of consistency, which may have a genetically determined component (Koolhaas and Van Reenen, 1999;Løvendahl et al., 2016). Similarly, cows may develop routines or preferences for being milked at specific hours of the day. The consistency of these routines are sometimes very strong, and the strength of it was shown to be both an acquired feature, but also to include some genetic variation (Løvendahl et al., 2016). However, little evidence is available on how different aspects of behaviour at robotic milking are related, and how these relate to key production traits at the phenotypic and the genetic level. The aim of the present study was to investigate preference consistency, in Holstein and in Jersey cows having access to two, three or four milking boxes, and their use of time of day to be milked, for repeatability over time and for genetic variation, and eventually for correlation to milk yield and time spent in the milking box. The quantitative genetic study was based on records from commercial dairy herds in Denmark.

Materials & Methods
Overall Design, Herds and Cows. This investigation was conducted as a cohort study using robotic milking data from 12 Danish dairy herds, having either Jersey or Holstein cows of mixed parities (1 -5) between October 2014 and September 2019. There were 2,407 Jersey cows with 28,952 records and 2,258 Holstein cows with 31,290 records. Cows were kept in groups having access to either 2, 3 or 4 milking boxes with free cow traffic, so that cows had a choice at every milking visit. The genetic relationships between cows within each breed was known and traced back 5 generations through the NAV database (NAV, SEGES, Skejby, Denmark). The Hostein pedigree consisted of 20,229 and the Jersey pedigree consisted of 14,543 animals. This study used non-invasive production data and was therefore exempt from animal experimental ethical approval. Milking data recording. At each milking the time of entry and exit from robot was recorded together with milk yield. From these data the time spent in the milking box (Boxtime, minutes) was calculated. The interval between milkings were expressed as milking frequency (i.e. milkings per 24 h, MilkFreq), because cows have voluntary access to milking boxes and thereby will have variable milking intervals. Milk composition data to calculate the energy corrected milk (ECM) (Sjaunja et al. 1990) was measured through the national recording organization RYK (RYK, Skejby, Denmark) at the Eurofins laboratory (Eurofins, Vejen, Denmark).

Milking box Preference Consistency Score calculations.
Preference consistency score (PCS) was calculated based on frequency of access's to each of the milking units in a given time period, here chosen at 15 days ("segment"). Cows accessed milking boxes voluntarily so it is assumed that the box being visited most frequently is the preferred one. If all boxes were equally used the frequency of first choice will be equal to 1/number of choices. The PCS is then calculated using:

PCS = (Actual_fraction -(1/number of choices)) / (1 -(1/number of choices)) (1)
Diurnal Milking Time Profile. The time profile of being milked at specific times of the day was here estimated based on frequency of milkings at each hour of the day, condensed in a single number indicating the variation. It was assumed that in case a cow has no hourpreference, the frequency in each hour-bin is equal to all other hour-bins. In contrast, if all or most milkings will fall into one or a few hour-bins, the deviation from equal use of hour-bins will be large. In between these extremes will be a continuum of intermediate cases. The deviation from equal frequency in all hour-bins was quantified using a Chi-square goodnessof-fit approach, and log-transformed to obtain an approximate normal distribution, better suited for analysis of variance. Statistical models and parameter estimates. Variance-covariance components were estimated within breed across lactation using uni-or bi-variate models in DMU software package (DMU version 6 Release 4.7; Madsen and Jensen, 2008): y = Xb + Zpeape + Zgag + Ie (3) where y is the vector of observed traits, b is a vector of fixed effects of herd-year-season, number of AMS units (2, 3 or 4), parity-group (first, second or older=3-5), segment (1 -20) by parity interaction. The matrix X contain the design of the fixed effect factors. The random effects are the additive genetic effect of cow (ag), the permanent (cow) environmental effect across segments (ape) and the random error effect (e).These have incidence matrices Zg being the relationship A-matrix and Zpe, being incidence cow within parity. I being the identity matrix and e the random error term. Random (co)variance components were assumed independent following a normal distribution N(0,σa 2 ), N(0,σpe 2 ), and N(0,σe 2 ). Genetic parameters were calculated from (co)variance estimates with associated asymptotic standard errors.

Results
The use of available milking boxes differed among cows, because some cows almost entirely used just one box while other cows spread their milkings over all available boxes. This was effectively quantified by the PCS, with values ranging between zero and 1.0, and an overall mean of 0.43 (sd = 0.28) and 0.41 (sd = 0.27) in Holstein and Jersey cows. This value is clearly larger than zero, where a zero value would indicate even use of all boxes. First parity Jersey cows had the highest PCS followed by first parity Holstein, and lowest values were found in the group of older Jersey cows. The PCS decreased with age in Jersey cows, but remained at similar level for all ages in Holsteins. The Time_profile varied among cows, and high values show that some cows make most of their milkings during few clock hours, whereas the variation show that other cows spread their milkings over several clock hours. The first parity Jersey cows, and to a lesser degree, the first parity Holstein cows condensed their milking events into fewer clock hours than older cows of both breeds. There was a general trend for the time profile to increase gradually during progression of lactation, and Jersey cows had always a slightly higher Time_profile score than Holsteins. The milking frequency increased during early lactation for young cows of both breeds, and then decreased gradually from about 100 DIM. In contrast, milking frequency decreased throughout lactation in all groups of older cows. The time spent in the milking box increased in early lactation for all breed and age groups, and then declined after reaching a peak at around 50 DIM. Heritability estimates for the behavioural traits and ECM are presented in table 1. The genetic and the individual level correlations between PCS and Time_profile (ri ~ 0.10 to 0.18) showed that cows who more consistently use a given milking box, also tend to have a stronger Time_profile, i.e. be milked more often at the same time of the day. Those cows that are more consistent also have somewhat lower milking frequency (ri ~ -0.18 to -0.16). However, consistency seem almost unrelated to total daily time spent in the milking box and neither to milk yield. Individual cows who show a stronger Time_profile were also spending less time in the milking box, and may have slightly lower milk yield (ri ~ -0.10). High milk yield was weakly related to Time_profile and even less to PCS, but positively with higher milking frequency and with more time per day spent in the milking box. Overall, the phenotypic correlations closely resembled the individual level correlations and correlations in the Jersey cows were almost similar to those obtained in Holsteins.

Discussion
Dairy cows entered the milking box on around 3 times per day in this study, and each time they had to wait for a vacant box, either being their favourite one, or another one of at least 2, and often 3 or 4 boxes. There were no obvious reasons for cows to prefer a specific box in general, as the available boxes within a group were used almost equally much, as also found in a previous study (Løvendahl et al., 2016). Assuming a specific cow had a preferred box she would have to wait until it was available, or to choose one of the other boxes. In the waiting area, she could also chase other cows away from the queue to get access to her preferred box. In this context, we interpret higher use of a specific box as preference and measure it over a short period of 15 d. In case only two possible choices are available this is reasonably simple, and a "consistency score" can be used (Hopster et al. 1998, Løvendahl et al. 2016. However, for this study some cows had 2 boxes, others had 3 or 4 boxes to choose between, so in order to adjust for number of choices an adjustment of the consistency score was implemented. This adjusted score, now called a preference consistency score PCS gave values between zero and 1.00 based on how many of the milkings were on the most frequently used box. One should note that the PCS itself resembles a repeatability coefficient, but here used as a response variable. This study showed that differences between cows in use of specific milking boxes are repeatable and to some degree also heritable. This extends earlier findings from a preliminary study only in Holsteins and only in herds with two milking boxes (Løvendahl et al, 2016), to include many more cows and two breeds, and situations also comprising 3 or 4 milking boxes. This study was primarily designed to quantify genetic and phenotypic co-variation in milking behavior traits as possible correlated response to selection for production, but not to include these traits in selection programs. Although more than 4,600 cows provided data, genetic correlations had large standard errors, whereas individual level correlations were used as proxies and gave similar results for the two breeds. Obtaining even larger cohorts or using genomic information may help resolve this problem in some future studies.