Culling reasons and risk factors in Estonian dairy cows

Background Culling is a major cost for dairy farms but also an essential part in managing herd productivity. This study aimed to identify the culling rates of Estonian dairy cows, identify the farmers’ stated reasons and risk factors for culling. This observational study used registry data of all cows from herds with ≥20 cow-years in 2013–2015. Cow lactation-level analyses included data of 86,373 primiparous cows from 409 herds and 177,561 lactations of 109,295 multiparous cows from 410 herds. Weibull proportional hazard regression models were used to identify risk factors for culling due to slaughter or death. Results The overall culling rate of Estonian dairy cows was 26.24 (95% CI 26.02; 26.46) per 100 cow-years. The most common reasons farmers stated for culling were feet/claw disorders (26.4%), udder disorders (22.6%), metabolic and digestive disorders (18.1%) and fertility problems (12.5%). Animal-level risk factors for culling were Holstein breed, older parity, lower milk yield breeding value, older age at first calving, longer previous calving interval, having assisted calving, stillbirth and birth of twins/triplets. Lower milk yield, somatic cell count over 200,000 cells/ml and fat/protein ratio over 1.5 at first test-milking after calving were associated with greater culling hazard during the lactation. Cows from larger herds, herds with decreasing size and higher milk yields had a higher culling probability. Conclusions This study emphasises the need for improved management of hoof health and prevention of mastitis and metabolic diseases. It is essential to ensure easy calving and good health of cows around calving in order to lower the culling hazard.


Background
Culling is defined as the departure of cows from the herd because of sale, slaughter, salvage, or death [1]. Culling is an important cost for dairy farms [1][2][3]. At the same time, culling is a way to increase herd productivity and profitability, as keeping diseased and unproductive cows might result in lower herd milk production and deteriorated reproduction. Keeping cows too long in a herd might also impair the herd's genetic improvement [3]. In order to maximize profitability, the proportion of voluntary culling (selling for dairy purposes or culling due to low production) should be highest among the total culling rate [1,3,4]. Previous studies indicate an ascending trend in the proportion of involuntary culling [4,5]. Culling rate, especially the proportion of involuntary culls, may also be considered as an animal welfare indicator [4,6,7].
During the last decades, culling rates have not increased; however, the longevity of cows has declined worldwide [8,9]. In Estonia, the average productive age of culled cows has decreased from 1113 to 1051 days between 2013 and 2018, respectively [10,11]. Longer productive lifetime would result in lower replacement costs, the opportunity to sell heifers [12,13] and the possibility of increasing the proportion of voluntary culling, bearing a desirable effect on herd profitability [14]. Return of the rearing costs up to first calving are covered roughly by the start of the second lactation [13,15]. Milk yield reaches the highest level after the third lactation [11], meaning that culling young cows is especially undesirable. Importantly, short lifespan is against consumers' expectations and is also associated with more detrimental environmental impacts [16]. Culling reasons have changed over the last decades, referring to the proportion of culling due to voluntary reasons (e.g., low milk yield), which has decreased and shifted into rather disease-related reasons [9].
Estonia is located in Northern Europe on the coast of Baltic Sea. Having a temperate climate condition the average temperature is + 5 C°being − 4 … -5 C°i n winter months and 15 … 18 C°in summer months. Estonian dairy cow population constitutes approximately 85,000 dairy cows, with roughly 67% of the dairy cows being housed in farms with more than 300 cows in 2018 (Estonian Agricultural Registers and Information Board, 2018). The average milk yield of Estonian dairy cows was 9785 kg in 2018 (Estonian Livestock Performance Recording Yearbook, 2018), taking second place in the EU (Eurostat, 2019). The majority of large Estonian dairy farms are loosehoused open-air barns using modern equipment and technology managed by hired labour. The Estonian dairy cow population is endemically infected with many cattle pathogens, including bovine herpesvirus 1, bovine respiratory syncytial virus and bovine viral diarrhea virus [17], digital dermatitis [18] as well as with main contagious mastitis pathogens Staphylococcus aureus, Streptococcus agalactiae and Mycoplasma bovis [19]. Due to these population characteristics, the results of previous studies investigating reasons and risk factors for dairy cow culling might not apply to these herds. Due to worldwide intensification of dairy production [20], the Estonian dairy cow population could serve as a valuable example, representing herds with high production level and mostly intensive keeping conditions.
In order to make use of opportunities that may be afforded by increasing cow longevity, there is a need for an overview about the main reasons as well as the risk factors that influence the probability of culling. This information would allow for the establishment of betterfocused hypotheses for further studies that aim to reveal causal associations between animal and herd-level risk factors and dairy cow longevity. The aim of the present study was to analyse the culling rates of Estonian dairy cows, identify the frequency of farmers' stated reasons for dairy cow culling, and investigate risk factors for cow culling.
The culling risk was highest during the early postpartum period (Figs. 1 and 2). By 30 days after calving the survivor probability of primiparous cows was 0.96 (95% CI 0.96; 0.97) and 0.94 (95% CI 0.94; 0.94) in multiparous cows. The average survival probability was 4.42% lower in multiparous cows compared to primiparous cows at 100 days of lactation. The survival probability dropped from 0.91 (95% CI 0.91; 0.91) to 0.89 (95% CI 0.88; 0.89) in primiparous cows and from 0.83 (95% CI 0.83; 0.83) to 0.78 (95% CI 0.78; 0.78) in multiparous cows between 200 to 305 days after calving (DAC), respectively. The survival probability declined by 16.1 and 23.5% between 305 and 600 DAC in primiparous and multiparous cows, respectively (Fig. 1). Around one year after calving, a small increase in the number of culling events was observed among primiparous cows and a stabilisation of continuous decrease of number of culling events over lactation was also identified in multiparous cows at that time (Fig. 2).

Farmers' stated reasons for culling
Culling reason was missing for 95 and 366 culling events in primiparous and multiparous cows, respectively. Main causes of culling due to death and slaughter were reported as feet/claw disorders (26.4% of culled cows), udder disorders (22.6%), metabolic and digestive disorders (18.1%), and fertility problems (12.5%). The importance of fertility problems as the reason for culling decreased after the second parity (Table 1).
In total, 50.7 and 46.5% of the culling events occurred during the first 100 DAC in primiparous and multiparous cows, respectively. Roughly one-third (33.9%) of the cullings took place at more than 200 DAC. Metabolic and digestive disorders were more prevalent up to 100 DAC (17.1 and 31.7% in primiparous and multiparous cows, respectively) compared to later stages of lactation (11.5 and 11.1% in primiparous and multiparous cows, respectively in 101-200 DAC). Feet and claw disorders (30.9 and 33.1% in primiparous and multiparous cows, respectively) and udder disorders (25.0 and 35.0% in primiparous and multiparous cows, respectively) were most common reasons for culling in the middle of the lactation period. Fertility problems were the primary reasons for culling in the late lactation stage (≥200 DAC), comprising 38.2 and 28.1% of all culls among primiparous and multiparous cows, respectively. Across all culling reasons, metabolic disorders were nearly twice as frequently the cause of culling in multiparous cows compared to primiparous cows. In primiparous cows, dystocia and low milk yield both constituted roughly 9% of all culling reasons during the first 100 DAC, being higher than that reported in multiparous cows (4.1 and 2.0%, respectively) ( Table 2).

Animal-level risk factors for culling
Descriptive statistics and univariable associations between continuous and categorical predictor variables and culling are presented in Supplementary Tables 1A and 1B, respectively.
Several culling risk factors were common for primiparous and multiparous cows. Estonian Red and Estonian  were associated with a higher risk for culling compared to giving birth to a single female calf. In multiparous cows, the culling hazard was also higher in cows that gave birth to a male calf compared to a female calf (HR = 1.05, 95% CI 1.03; 1.07). There was a timedependent effect of assisted calving to culling hazardthe negative effect of assisted calving on culling hazard was higher during the first seven days post-partum compared to a later period. Also, higher age at first calving was associated with a higher culling hazard (HR = 1.07, 95% CI 1.06; 1.07 and HR = 1.02, 95% CI 1.02; 1.03 in primiparous and multiparous cows, respectively) but the association was not linear (Tables 3 and 4). In multiparous cows, the culling hazard increased with parity. Somatic cell count (SCC) over 200,000 per ml and lower milk yield at last test-milking before dry-off were factors associated with a higher culling hazard during the following lactation (HR = 1.21, 95% CI 1.18; 1.23 and HR = 0.99, 95% CI 0.989; 0.99, respectively). Robot milking system at last test-milking before dry-off was a protective factor for culling (HR = 0.83, 95% CI 0.78; 0.89) compared to other milking methods (Tables 3 and 4).
The graphical assessment confirmed meeting the proportional hazard assumption of tested categorical variables in a log cumulative hazard plot as well as the overall fit of the models assessed by plotting Cox-Snell residuals against the cumulative hazards of individual observations at their failure times.

Culling rates and farmers' stated reasons for culling
There is no single optimal culling rate that is applicable to all herds for all years due to a variety of economic factors, farm capacities, individual cow factors, morbidity and mortality rates within the herd, availability of replacements, biosecurity considerations, etc. [1]. Not all studies have presented the culling rates due to death/ slaughter excluding sales. In the current study, the average culling rate of Estonian dairy cows was 26.24 per 100 animal-years, excluding selling from the culling definition. In selected regions of the United States, the average culling rate was 31.6% in 1999 [3] but was 27.7% in Pennsylvanian herds in 2005 [21]. In Canada, Haine et al. [22] reported an average culling rate of 32% over the 2001-2010 decade and a dairy sell rate by 60 days in milk of 3.2%. The average culling rate of cows due to slaughter/death was 25.4%, ranging between 23% (in 2007) to 28% (in 2010) in Dutch dairy herds [23]. This shows that, on average, culling rates of the Estonian dairy cow population are mostly comparable to that reported in other countries.
Farmers' stated reasons for culling were analysed in the present study due to the absence of more sound  representative data. Care is needed when interpreting this data due to somewhat overlapping categories. Some cows might also have more than one reason for culling, whereas the ELPR system allows farmers to mark only one reason for each culling event. Also, the stated culling reasons might be the consequences of the primary disorder that might differ from what was reported by the farmers. Therefore these results are rather indicative and further studies including necropsy, meat inspection and laboratory data together with animal and disease history could reveal more reliable results. The four most common reasons for culling due to death and slaughter were: "feet/claw disorders", "udder disorders", "metabolic and digestive disorders", and "fertility problems", the order of frequency being broadly similar to those described by other authors [24]. The proportion of "feet/claw disorders" as the reason of culling was somewhat more prevalent compared to what was reported in studies performed in other European countries [6,25]. As small herds with less than 20 cow-years were excluded in this study and majority of the Estonian dairy cows were housed in large herds [26], the conditions of large dairy herds and their effect on the cow hoof health were probably over-represented in the present study compared to other studies. In large Estonian dairy farms, cows are mostly housed in freestalls, and the latter is shown to be associated with an increased risk for lameness relative to other housing systems, including tie stalls and straw yards [27]. Additionally, other factors accompanying freestalls such as prolonged standing time due to milking [28], overstocking [29] or poor stall design [28] can increase the risk of lameness. Metabolic disorders were a more common reason for culling in multiparous compared to primiparous cows. Roberts et al. [30] concluded that primiparous cows may have a different physiological response to post-calving metabolic challenges. The need to balance between energy demands for growth and milk production may have an effect on more effective fat mobilization before health and productivity are compromised [30].
A minor increase in the culling hazard occurred about a year after calving. As fertility problems were the   [31]. The Estonian milk recording register system allows farmers to report only one reason for each culling event. Therefore it might be possible that some proportion of cows that were culled due to other reasons than fertility during the last third of lactation were also non-pregnant due to suffering from chronic health disorders or according to the farmers' decisions. Loss of a cow at its first lactation is economically most devastating and therefore undesirable for the dairy farmer [1]. In order to lower the culling rate of primiparous cows, measures promoting the good health of feet and udder are of utmost importance. In primiparous cows, dystocia and low milk yield both constituted roughly 9% of all culling reasons during the first 100 DAC, being nearly twice as high as that reported in multiparous cows. According to Mee [32], the fetopelvic disproportion is the predominating risk factor for dystocia in primiparous cows, therefore bull selection as well as heifer nutrition and development might be critical factors lowering dystocia that leads to culling in primiparous cows.
We also identified some differences in culling reasons over lactations. Fertility, as a reason for culling, decreased in importance with each parity, indicating that probably more resilient cows in terms of breeding capability remained in the herds. Metabolic and digestive disorders as well as udder disorders were more frequently stated as the reason for culling in older cows. The fact that each calving event cumulatively adds the risk of suffering postpartum diseases, such as mastitis and ketosis [33,34], might also explain this identified association.

Animal-level risk factors for culling
In the current study, several common risk factors for culling in primiparous and multiparous cows were identified. Holstein breed cows had a significantly higher culling hazard compared to Estonian Red and Estonian Native breed cows. Holstein breed cows have a higher milk yield in Estonia [11]. Concomitantly, pure Holstein breed cows are more susceptible to production-related diseases [35] and have poorer reproductive performance than crossbred cows [36], thus being more prone to culling.
We identified that higher individual milk yield breeding value was a protective factor for culling, which may be explained by farmers trying to keep cows with good genetic merit. A milk yield breeding value could be calculated for cows who have at least two test-milking results available and whose sire has obtained a milk yield breeding value. The category "missing" also included cows that were culled during early lactation, which might be the cause of high culling hazard among primiparous cows who had no breeding value in the dataset.
In agreement with other studies, it is essential to pay attention to predisposing conditions and transmission of infectious diseases that might be associated with the incidence of stillbirth or abortion, as these are important risk factors for cow culling and longevity [37,38]. Additionally, higher culling risk due to giving birth to a male calf compared to a female calf indicates that birth weight of an offspring might be an important factor in terms of culling via increasing the probability for dystocia. This is known to be associated with a higher mortality hazard in cows [33] as well as an increased risk of post-partum diseases eventually leading to culling [39]. Furthermore, calving for the first time at a higher age was associated with a higher culling hazard. Interestingly, the association was also present in multiparous cows, suggesting a possible long-term impact. Heifers might calve at an older age due to management factors, health disorders, feeding management, or due to herd breeding strategies [23,40], and the reason for the old age at first calving may be more important for the culling risk than the age per se.
In multiparous cows, a longer previous calving interval was associated with a higher culling probability at next lactation. Prolonged calving interval might be related to negative energy balance and diseases associated with the early post-partum period, delaying conception. Due to the possible recurrent propensity of post-partum diseases at the cow level, the undesired impact of a longer calving interval might manifest at the next lactation. In addition, a longer calving interval allows cows to gain more weight, which may be a risk factor for developing post-partum diseases, increasing the culling hazard [41].
High somatic cell count at last test-milking during the previous lactation and/or at first test-milking after calving, indicating the presence of subclinical intramammary infection, was a factor associated with higher culling risk during lactation. The presence of clinical or subclinical mastitis is a known risk factor for dairy cow culling [39]. Lower milk yield at the end of the previous lactation or at the first testmilking of ongoing lactation was also associated with a higher culling probability. It has been shown that, in general, farmers are more eager to cull lowproducing cows [39]. Still, low milk yield soon after calving might be associated with an underlying disease. The high milk fat/protein ratio is a valuable indicator of ketosis in the early post-partum period [42], and it is concomitantly related to increased probabilities of developing displaced abomasum, retained placenta, metritis, clinical endometritis and clinical ketosis as well as higher culling probabilities [43,44].

Herd-level factors associated with culling
On average, the cow culling hazard was higher in larger herds. The association between herd size and the health and welfare of dairy cows is complex and includes the impact of several factors, including facilities, management and operational factors [20,45]. According to previous studies, the positive association between the incidence of metabolic diseases and herd size have been found [46]. Also, the presence and dynamics of infectious diseases and different biosecurity management could explain higher culling risk in larger herds [23]. According to Gieseke and co-authors [47], housing conditions and management practices have a greater effect on cow welfare than the herd size itself and more research is needed to identify factors in large farms that affect animal health. Simultaneously to the European Union, the milk price dropped suddenly in the autumn of 2014 in Estonia [48], resulting in a 5.2 and 11% reduction of the number of dairy cows and herds, respectively, in Estonia during the year of 2015 [49]. Due to this, change of herd size was controlled as a factor in the statistical models to account for its confounding effect.
Higher herd average milk yield was associated with increased culling hazard of individual cows in our study. Several previous studies have found that higher milk yield has an adverse effect on the cows' resistance to diseases, as it is correlated with presence of the clinical mastitis, reproductive diseases [50], and other postpartum disorders [43]. Although high milk yield and the high genetic potential for milk production are often been blamed for the short longevity of dairy cows, this does not always seem to affect cow longevity [51]. On the other hand, due to uncertain causality, herds with higher milk yield might also have better reproductive performance, allowing more cows to be culled.
The current study showed that herds with a longer average calving interval had a lower risk of culling. A longer calving interval may be the consequence of fertility problems in the herd or the result of a voluntary decision of the farmer to delay with breeding after calving and thus to extend the lactation period [52]. Although a shorter calving interval is considered as economically optimal [53], an economic benefit in extending lactations in high-yielding cows was also found [54]. In the study by Allore and Erb [55], a lower risk of culling for reproductive failure was present in herds with an extended voluntary waiting period. Still, due to a cross-sectional study design, it is impossible to draw causal inferences, and the identified association might result from farm lower thresholds for culling non-pregnant cows [39].
In the current study, the cow culling hazard was also negatively associated with herd average lactation number that cumulatively aggregates the individual animal culling hazards. Using robot milking system was a protective factor for culling compared to other milking methods. Additionally, cows that were milked three times a day at first test-milking had a lower culling probability throughout the lactation compared to cows that were milked twice a day. Farms with automatic milking systems differ from those with other milking systems in many aspects, such as environmental conditions, feeding management, grouping policies, etc. Unfortunately, it was not possible to discriminate milking method, but only the milking frequency, meaning that the identified associations could be affected by other factors.

Validity and limitations of the study
The present study included lactation-level records of all cows from all herds that had at least 20 cow-years in years 2013-2015 in Estonia and participated in the milk recording system. According to the ELPR [26], 94.1% (in 2013) to 95.4% (in years 2014 and 2015) of the Estonian dairy cow population was enrolled in the milk recording system, probably leaving out smaller farms that produced milk for their own consumption. After skipping herds with < 20 cow-years, we emphasize a good external validity of this study for medium or large sized herds, whereas the study results should not be extrapolated to small dairy holdings.
When analysing farmers' reported reasons for culling, a reporting bias might be present. In Estonia, farmers are allowed to report only one reason for each animal exit to the milk recording register. Still, studies have shown that in many instances, farmers report more than one reason of culling when allowed [56,57]. To our knowledge, there are no studies that investigate the farmers' behaviour in reporting culling reasons. We assume that farmers report the main and most obvious reason of culling at the time the cow is leaving the herd. However, this might not be the primary or ultimate disease or disorder which leads to culling.
In the EU, there is a harmonized mandatory registration and reporting of animal births, movements, and deaths [58] making registry data reliable for research purposes. As participation in the milk recording register is voluntary, the reporting in that system might not always be precise. Although the animal registry and milk recording registry makes crosschecks in their data, small discrepancies were found when comparing the cows' exit dates in the two registries. Also, some cows had very long lactations and low number of culls occurred far away from the latest calving (Fig. 2) in which a new calving date might have remained unreported in the registry. Still, due to the high sample size of this study, these aberrations probably have no effect on the overall results and conclusions.

Conclusions
Culling rates of Estonian dairy cows housed in farms with more than 20 cows were comparable to that reported in other countries. More attention should be paid to prevention of feet/claw disorders and mastitis, because they accounted for roughly half of the reported culling reasons. The early post-partum period was the period bearing the highest risk for cow culling, with metabolic diseases and feet/claw diseases being the most frequently stated reasons for culling at that time. Our analyses confirmed the importance of easy calving and ensuring good health around calving to avoid the loss of a cow. Due to the intensification of dairy production and the identified higher culling hazards in larger herds, further studies should concentrate on factors that contribute to cow culling and longevity in large commercial herds. Along these same lines, a widened list of herd factors, including animal housing, environment, management as well as farm workers' attitude and motivation, should be studied helping to understand the complex problem of cow culling and combat the problem of reduced cow longevity.

Datasets
All data used for the present study was based on recordings in the Estonian Farm Animal Register kept by the Estonian Agricultural Registers and Information Board and data of Estonian Livestock Performance Recording Ltd. (ELPR, national voluntary animal production system involved with monthly test-milkings including data of around 95% of dairy cows in Estonia, called "milk recording register" hereafter). Cow lactation level records for the period between January 1, 2013 and December 31, 2015 were collected from the milk recording register. Data was required of all Estonian dairy herds that participated in the milk recording register and had a herd size of at least 20 cows at the start of the study period. The initial datasets included lactation level data of 86,459 primiparous and 109,314 multiparous cows (177,712 lactations) from 409 and 410 dairy herds, respectively.
The lactation level observation period started at the day of calving and ended at the date the cow left the herd, had a new calving, or up to the end of the study period (December 31, 2015). Due to the aim of reducing the impact of voluntary culling in the analysis, the definition "culling" established by [1] was modified in the sense that sold cows (cows that were sold to another farm for productive life) were considered as right censored observations, thus excluding sales from the culling definition. Therefore the definition of culling in this study includes on-farm mortality (unassisted death and euthanasia) and slaughter (sending cows to an abattoir).
Farmers are obliged to report all movements and exits of cattle to the Estonian Farm Animal Register within 7 days [59]. In the Farm Animal Register, it is specified for each exit event whether the animal died on-farm, was euthanized, got lost, was slaughtered or was sold/ exported. Farmers report the reason of exit in the milk recording register, choosing one reason from the list of 24 pre-defined categories, including selling, old age, low milk yield, mastitis, metabolic disorders, abortion, etc. (the complete list is given below in the section "Data editing"). The definition of these categories was not provided to the farmers by the registry. In order to analyse the proportion of deaths (including euthanasia) and slaughter as well as include a reason for each exit event, the data of two registries was merged by cow identification numbers. Due to the fact that the analyzed datasets were on a lactation level and Farm Animal Register data on cow level and one animal might have more than one exit during its lifetime (e.g., sold and slaughtered/died), the exit dates of the two datasets were compared in order to connect a correct exit type (death / euthanasia / selling / slaughter) and farmers' stated reason for animal exit. Reporting reasons of cattle exit in the milk recording register by farmers is voluntary, therefore it was expected that the dates of cattle exit in two datasets would slightly differ. When joining the datasets of two registries, differences of exit dates up to seven days were allowed and were considered to be the same exit event.
The data from Farm Animal Register consisted of farm identification number together with the owner's name, location of a herd, animal identification number, birth date, sex and breed of an animal, exit date and type of exit. Milk recording register data included animal, lactation and herd level data. A more detailed description of the data provided in the milk recording system is described in [60].

Data editing
Due to inexplicable short calving intervals (calving interval < 290 days and no stated abortion or stillbirth in the registry), 86 and 151 lactation records were removed from the datasets of primiparous and multiparous cows, respectively.
The category "other reasons" included farmers' stated reasons that were rarely reported or not reasonably associated with death and slaughter, thus considered to be recording mistakes. In order to determine whether the cow was purchased or not, the birth farm identification number was compared with the farm number in which the cow was located. Milk yield breeding value was not available for all cows (in case the sire of a cow had not received a milk yield breeding value or the cow had less than two test milking results available), thus the variable was categorised into roughly equally sized groups with a separate category for cows with missing information. Due to a low number of Estonian Native breed cows, these observations were merged with Estonian Red breed cows. Counties were combined into four regions: Northeast, Southeast, Southwest, and Northwest regions. Four seasons were created based on the calving dates of the cows: winter (December-February), spring (March-May), summer (June-August) and autumn (September-November). The variable "change of herd size" was created by comparing the number of cows in the herd between 2013 to 2015. Based on that, five categories were created: no considerable change in the number of cows (± 5%), decrease > 5 to 15%, decrease > 15%, increase > 5 to 15% and increase > 15%.
Data related with test-milkings was categorized when a biological threshold for discriminating normal from pathological conditions was available. The variable "milk somatic cell count" had two categories separated at the level of 200,000 cells/ml [61]. A threshold of 1.5 was used to dichotomise "milk fat/protein ratio" [43] and 6.78 mmol/L was the limit value for "milk urea content" [62].
For the herd-level production and reproduction data, three-year averages were calculated (Supplementary  Table 1A).

Data analysis
A descriptive statistical data analysis was performed to identify farmers´stated reasons for culling. The analysis was performed separately per parity as well as by stage of lactation.
Risk factor analysis for cow culling (on-farm death, euthanasia or slaughter) was performed separately for primiparous and multiparous cows due to different sets of variables used. A complete description of the statistical methods used can be found in [60]. Briefly explained, the Weibull proportional hazard regression model specifying herd as a random factor (gamma distributed frailty effect) was applied specifying culling as a failure event. Lactation level observation period started at the day of calving or at the day of purchase of a lactating cow. In order to account for the left truncation in the analyses [63] the "start date" was accounted as January 1, 2013 for observations that started before that date by specifying the "stset" options in STATA MP 14 software (StataCorp LP, College Station, USA). Observations that ended with new calving, selling, or the end of the study period were right censored. Due to several lactations of multiparous cows in the dataset during a study period, the animal identification number was specified as the multiple-record ID variable when setting time-toevent data. A Kaplan-Meier survival curve was created to compare the survival probability of primiparous and multiparous cows over lactation by plotting estimates survival curves. Culling rates together with 95% confidence intervals were calculated by using the Mantel-Haenszel method. Estimated cumulative survivor hazard function was calculated to describe the culling probabilities at different time-points over lactation. In order to identify the difference of culling hazards between primiparous and multiparous cows, a univariable proportional hazard Weibull regression model was composed including cow parity (primiparous / multiparous) as a fixed effect and herd as a frailty term.
After univariable screening of the predictors in their association with culling probability all variables with a liberal p-value of < 0.05 and not collinear to each other were considered as candidates in multivariable analysis. A multivariable Weibull model was built by removing insignificant predictor variables from the model by manual backward elimination technique. Due to the variable "Assisted calving" having a time-dependent effect, the lactation-level observation period was split on the seventh day post-calving using the 'stsplit' command. Due to the high power of the analyses coming from a large sample size, the predictor variables included in the final model had to be associated with the culling hazard at a 1% significance level.
Due to differences in the beginning of the observation period, a separate model was conducted to analyse the associations between first test-milking results and culling during lactation controlling for confounding variables ("days in milk", "breed", "parity", "calving year", "herd