Genetic analysis of productive life length in Holstein dairy cows using Weibull proportional risk model

In recent decades, there has been a downward trend in length of productive life (LPL) in Holstein cows across industrial dairy herds. This study examined the factors that might influence LPL and estimated the genetic parameters of LPL in Holstein dairy herds in Isfahan province, Iran. LPL is defined as the number of days between the first calving and the end of recording. Data consisted of 35 137 records of productive life from registered cows that started first calving between 1991 and 2012. Cows that remained alive at the end of the study were considered right-censored. The average lifetime for culled and censored cows was 938 and 1003 days, respectively. A survival analysis was applied using a proportional risk model with a Weibull distribution. Milk production was divided into five groups, where the culling risk of cows with a milk yield of less than 1.5 standard deviations (SD) of the mean was 3.5 times greater than the culling risk of high-producing cows producing more than 1.5 SD above the mean. Results showed that culling risk increased almost linearly beyond the milk production groups. Furthermore, somatic cell count and age at first calving significantly increased the culling risk across the herds. The results for the combined effect of parity× stage of lactation showed a decrease in culling risk during the first calving, and an increase during the further parities. Moreover, a higher age at first calving was observed, reflecting a lower risk of culling. Estimated heritability were 0.074 and 0.18 based on a logarithmic scale and original scale, respectively. According to the results, use of Weibull models showed that the proportional culling risk was higher in low-production cows, but a higher risk ratio was revealed in high-milk-production cows. However, there were some fluctuations in genetic trends, but an overall increase was observed in LPL which will lead to a longer LPL of Holstein cows in Isfahan province.


Introduction
Length of productive life (LPL) is a highly desirable functional trait that significantly affects overall profitability in dairy herds (Van Arendonk, 1991).Longevity of dairy cattle is a trait composed of many different factors that accumulate over time (Sørensen et al., 2008).With cow's increased longevity, milk performance of the herd could increase for two reasons.First, a greater proportion of the culling decisions are based on production.Second, the proportion of mature cows, which produce more milk than young cows, is increased (Allaire and Gibson, 1992;VanRaden and Wiggans, 1995).With herd sizes and costs of labour increasing and production and management systems becoming more intensive, fitness traits such as health, fertility and survival rates will become more and more crucial for the dairy industry (Pedersen and Christensen, 1989).LPL is currently one of the most important functional traits in dairy cattle breeding schemes.A longer productive life decreases the costs of herd replacement and leads to higher production from larger herds as older cows produce more milk (Vukasinovic et al., 1997).The definition of two herd-life traits was reported by Ducrocq et al. (1988): true herd life (THL), or actual days from birth to disposal, and functional herd life (FHL), rep-resenting the ability of the cow to avoid involuntary culling (functional herd life adjusted for 305-day milk production).Otherwise, FHL should reflect involuntary culling (Ducrocq et al., 1998).Functional herd life was defined as the length of life from first calving to death, culling or censoring.The LPL trait has been defined as the number of days from first calving to culling, death or censoring (Ducrocq et al., 1988).In this context, censoring means that the animal was either still alive at the time of analysis or its status was unknown, based on the reported disposal code for each cow.In recent research, several factors affecting LPL have been identified, including milk production, parity × stage of lactation, somatic cell score and age at first calving.
The level of milk production has been reported as one of the most important factors affecting the LPL (Vukasinovic et al., 1997;Sewalem et al., 2005a).An LPL by Raguz et al. (2011) showed that milk production influenced productive life significantly.Previous reports have found that, as the risk ratio decreased, milk production increased (Vukasinovic et al., 1997;Chirinos et al., 2007;Raguz et al., 2011).Van der Linde and De Jong (2003) reported that cows with low milk production are exposed to a culling risk 4 to 5 times higher than average milk producing cows.
Several studies have examined the interaction of LPL with parity × stage of lactation (Van der Linde et al., 2006;Chirinos et al., 2007).These studies argued that the risk of culling increases across and along parity.Ducrocq et al. (1988) showed an increase in the trend of risk ratio from the second parity.However, many reports have indicated a drop in risk ratio between the parities (Bonetti et al., 2009;Raguz et al., 2011). M'Hamdi et al. (2010) highlighted that dairy cows have a higher culling risk at both the beginning and end of a parity period.
A Belgian study showed that between days 5 and 14 after the first calving culling risk increased by 11 % due to a high somatic cell score (SCS) and by 32 % due to udder problems (De Vliegher et al., 2005).A study on Canadian Holstein cows showed the effect of high SCS on longevity was greater than any other factors (Sewalem et al., 2005a).The cows with higher SCS are approximately 4 times more prone to culling risk and a higher culling rate was observed in cows with deep udders (Van der Linde and De Jong, 2003).
Numerous studies comparing age at first calving and LPL have found that the risk of culling increases with increasing age of cows (Vukasinovic, 1999;Chirinos et al., 2007;M'Hamdi et al., 2010).However, Ducrocq (1994) analysed the data and concluded no significant effect of the age at first calving on LPL.M' Hamdi et al. (2010) observed a linear increase in relative culling risk proportional to the increase in age at first calving.
Longevity traits should be analysed using appropriate models.Records can contain data on cows that ceased to have production records (uncensored) or are still alive (censored).The usual approaches adopted for LPL records which are based on the threshold model are linear and survival mod- els (Boettcher et al., 1999;Caraviello et al., 2004;Sewalem et al., 2005a).Random regression linear models are not suitable for analysing LPL as censured data cannot be treated appropriately and the non-linearity of the factors cannot be fully considered (Caraviello et al., 2004).Survival analysis models are preferred for genetic evaluations of LPL.The implementation of survival analysis techniques and genetic evaluations for longevity enables proper statistical treatment of censored records and appropriate analyses of longevity data (Vukasinovic et al., 1997).With this method it is possible to use censored data, while environmental and management changes are used for a specific period (Allaire and Gibson, 1992).These models are referred to as proportional risk models, which can be subdivided into two models: the Cox model (semi-parametric) and the Weibull model (parametric).This study aimed to use survival analysis and examine the factors affecting functional herd life (FHL) of cows in Isfahan, Iran, using the Weibull proportional risk model to estimate breeding values for functional herd life.

Materials and methods
Records of Holstein cows born between 1991 and 2012 were obtained from the Vahdat Industrial Agriculturists and Dairymen Cooperative, Isfahan province.The total number of records was 50 090, consisting of 34 273 daughters sired by 1390 sires.Cows that reached parity six onward were considered right-censored.Herds fewer than 20 evaluated cows were omitted.Cows of at least the first parity were used, and age at first calving ranged between 20 and 40 months.In total, FHL records of 35 137 cows were considered.FHL was defined as the age at first calving until the last collected record, culling/death or data censorship.A summary of statistics for FHL (Table 1) shows the percentage of the cows culled due to calving difficulties on the first day; therefore, a FHL of 1 day was considered in this study.The last parity was considered for those cows whose culling date was missed.In this case, cows were assigned as culled only if the time interval between the end of the last parity and the date of recording exceeded 365 days.Three types of cows were excluded in this study: (1) sold cows, (2) cows without records, and (3) cows transferred to other herds.Sires with one daughter in a herd were removed from the pedigree.
The analysis for functional herd life was carried out using the Weibull risk sire model.
where h ij klnmp (t) is the risk function of a cow t days after calving; h 0 (t) is the Weibull baseline risk function with the shape parameter of ρ and scale of λ; hys i (t) is the random effect, which is dependent on ith time of herd-year-season assumed under log gamma distribution with gamma parameter (Y, Y ); p j (t 1 + t 2 ) is interaction between j th milk parity and stage of lactation; t 1 is days of the first milking records; t 2 is days in the stage of lactations; AFC k is the fixed timeindependent effect of the age at first calving; m l (t) is the fixed time-dependent effect of the relative milk production; D m (t) is the time-dependent effect of SC; and α n is the random sire effect.Genetic parameters and genetic evaluations were predicted based on the Weibull method (Ducrocq and Solkner, 1998).Risk ratios for given effects were computed using Survival Kit 6.0 software (Ducrocq et al., 2010).Following the designed algorithm in this software, the records with known longevity and low FHL limit were used.Hence, the records were considered uncensored data if the cows were either culled or died for any reason.Therefore, censoring the records represented the cows were sold, exported or leased to other herds.Both Cox and Weibull models were implemented in Survival Kit, and they could be used for continuous and discontinuous (time-dependent) variables.
The heritability on the log scale was calculated (Vukasinovic, 1999) as where (π2/6) is the variance of the standard extreme value distribution corresponding to residual variation in a normal distribution model and the heritability on the original scale (Yazdi et al., 2002) as where Var(s) is the sire variance.

Results and discussion
The SC, milk production, age at first calving and milk period showed significant effects on LPL (P < 0.001) according to the likelihood ratio test (similar to the type III sum of squares in linear models).The results are presented below according to the proportional risk pattern and the relations of each factor with the FHL was described.Test statistics for the likelihood ratio test (−2 LL), DF (degrees of freedom) and Maddala's R 2 (measure of explained variation) statistics for the data set are shown in Table 2.The Maddala's R 2 statistics for the model ranged from 0.11 to 0.37.The factor showing the largest contribution to the likelihood was the parity and stage of lactation effect, followed by the milk production effect and the age at first calving.

Relationship between SC and proportional risk of culling
In this study, one of the major factors of culling was mastitis.Hence, improving udder health can reduce the involuntary culling in dairy herds.Results indicated a direct relationship between the increase in SC and proportional culling risk (Fig. 1).Proportional culling risk at class 6 was 4 times greater (more than 1000 SC in 1 mL of milk) than class 1 (lower than 200 SC in 1 mL of milk).The difference between classes 3 and 4 was not significant.Caraviello et al. (2005) found that the cows with 700 000 SC were 3 times more exposed to culling risk than cows with the average SC.De Vliegher et al. (2005) revealed that this phenomenon was more obvious in cows at the first milk period.According to Chirinos et al. (2005), a direct relationship exists between SC and FHL, where a gradual increase in SC enhances culling risk.

Relationship between milk production and proportional culling risk
Milk production has an important role in the length of productive life.Results in Fig. 2 showed culling risk among cows whose production was less than 1.5 SD was 3.5 times greater than for highly productive cows at more than 1.5 SD above the mean.This study was similar to others where culling risk among highly productive cows is greater than that of to the average group with milk production between −0.5 and 0.5 SD (Boettcher et al., 1999;Haile-Mariam et al., 2003;Sewalem et al., 2005a).Groups 4 and 5 were not significantly influenced by the proportional culling risk and thus it is deduced that cows that produce more milk were not subject to culling (with the exception of the report of Raguz et al., 2011).Vukasinovic et al. (1997) showed that lowproduction cows were at higher risk than herd-mates with an average milk production, while high-production cows were less likely to be culled.Another study by Vukasinovic et al. (2001) showed that cows with less than 80 % of the herd milk average within each parity were at 3 to 4 times higher risk than herd-mates with an average production.Similar findings have been published by Vollema and Groen (1998).

Interaction between parity × stage of lactation and proportional culling risk
The effect of parity × stage of lactation on functional herd life is shown in Fig. 3.In the current study the highest culling risk observed in parity one occurred in the second stage of lactation and repeated in the last stage of the parity.This culling risk could be due to a peak in milk production when a negative energy balance is developed, leading to either metabolic disorders or diseases such as mastitis and acidosis (Roche , 2006).In addition, for the chance of culling at parity six, culling was increased with the stages of lactation, e.g. higher than the risk observed in first lactation (with the exception of the second stage of lactation; Fig. 3).The risk of culling in the second lactation and the stage of lactation by 150 days was twice as high in comparison with the first lactation at the same stage (Fig. 3).Generally, the combined effect of parity × stage of lactation showed a decrease in risk during the first parity and an increase in risk during later parities, with an overall decrease over parities.Conversely, an increase in risk over the parity was observed by Van der Linde et al. (2006) and Chirinos et al. (2007).Vukasinovic et al. (1997) expected that a cow finishing parity was at a much higher risk of culling than one placed in the early or middle stages of parity.According to M'Hamdi et al. ( 2010), cows at the beginning and end of parity as well as at the end of succeeding parities become more vulnerable to culling.Strapakova et al. (2014) stated that the subsequent risk ratio increased moderately at the end of the first parity and that the trend of culling risk was similar across and along each parity.Potocnik et al. (2011) reported that risk level at the end of parity can increase significantly 271 days after the start of milking and during dry periods, which is in agreement with our findings.Vukasinovic (2001) found that cows in the first parity were more vulnerable to culling than those at parity four and also increased gradually after that point, which indicated different selection criteria in the first parity; culling risk also increased after this parity.The interaction between the period and stage of milking indicated that proportional culling risk increased in a linear manner with the length of each period of 300 days after calving.The intensive selection for culling of non-pregnant cows corresponds to dry periods.Novaković et al. (2009) reported that an increase regarding the proportional culling risk was observed only in the first milk period (30 days), indicating an increase in the intensive selection of this period.

Relationship between age at first calving and proportional culling risk
This study indicated that, in general, a higher age at the first calving showed a lower risk of culling (Fig. 4).This is in agreement with other studies (Vollema and Groen, 1998;Vukasinovic et al., 2001).However, age at the first calving in this study was shown to have a lesser effect on LPL.The proportional culling risk among cows with the first calving after 34 months of age was twice that of the cows with first calving at 20-22 months.Cows calved between 26 and 34 months showed a similar proportional culling risk.Proportional culling risk in the current study was higher in those cows with low milk production and high age at first calving, which was usual in the herds that go through changes in their size on an annual basis (M'Hamdi et al., 2010).In this way, the results of this study were in agreement with Cox and Oakes (1984) and Vollema and Groen (1998).Strapakova et al. (2014) confirmed an increasing risk ratio with respect to the increasing age at first calving.Vukasinovic et al. (2001) observed no notable difference between groups of different ages at first calving.This same study reported a slight increase in the proportional culling risk for old cows with first calving at 37 months of age, and although the first calving age was not significantly different for different groups of cows, results revealed a relatively higher proportional culling risk among older cows compared to those first calving at more than 37 months (Vukasinovic et al., 2001).

Genetic parameters
Estimating genetic parameters with Weibull models could improve expected genetic progress if the selection to increase herd life is assigned in breeding objectives.Table 3 shows the estimated heritability based on the log and original scales for productive life.As discussed previously, ρ = 1 indicates consistency in culling risk as FHL increases.Furthermore, ρ < 1 indicates a reduction in proportional culling risk and ρ > 1 indicates an increase in proportional culling.The ρ obtained in this study (Table 3) showed that proportional culling commonly existed across the herds.Caraviello et al. (2004) and Vukasinovic et al. (2001) reported ρ in the range of 0.56 and 2.25.
Heritability may change depending on the estimates of ρ and scale (λ).In the present study, heritability estimates were 0.074 and 0.18, based on log and original scales, respectively.Studies on other populations using survival analysis showed higher heritability estimates compared to those resulting from linear models.For example, in Canada, heritability estimates using survival analysis and linear models were 0.04 and 0.12, respectively (Sewalem et al., 2005b).Based on linear and Weibull methods, Caraviello et al. (2004) reported that heritability estimates for FHL were 0.07 and 0.18, respectively.Strapakova et al. (2014) reported heritability on the original scale as being 0.13 in a Slovakian Holstein population.Similar or lower estimates were reported by Caraviello et al. ( 2004).

Genetic trends
An attempt to estimate the genetic trend for sires was made by grouping sires according to their year of birth.Figure 5 shows considerable fluctuations in the breeding values of two different years, and no significant changes were observed even in very low regression coefficients.Hence, based on the variation of the obtained breeding values, it is possible to increase the lifetime of cows by selecting cows of higher breeding value.Figure 5 indicates that the genetic merit for  1991 1992  1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008  this trait has been improved across the herds in the current study, and this will lead to a longer FHL of Holstein cows.Since the lower breeding value is an indication of longer LPL, the mean of the breeding value shows that LPL has changed across the times.However, most breeders increase the FHL of cows by choosing high-production cows.The previous results showed that culling risk of low-production cows was 3.5 times greater than that of high-production cows.Dürr et al. (1999) reported that, from 1977 to 1988, under Canadian production circumstances the genetic trend of herd life dropped, which was favourable because ETAs (estimated transmitting abilities) are expressed as relative culling rates.

Conclusions
Weibull models have been considered for genetic analysis of productive life in Isfahan province because of censored data, time-dependent variables and the rate of culling.In addition, survival analysis methodology can be examined to evaluate the impact of environmental and genetic factors on the risk of culling and to estimate the genetic parameters for longevity.The proportional culling risk in low-milk-producing cows was considerably higher than in high-milk-producing cows.Furthermore, the risk increased with higher age at first calving.The highest chance of culling occurred at the beginning of the first parity and at the end of the following parities.Therefore, both parity and stage of lactation should be used as a combined effect in the model to correct for different behaviours of the first and following parities.For the present data set the estimated heritability of functional longevity on the logarithmic scale was 0.074, while on the original scale it was estimated as 0.18.Although there were some fluctuations in genetic trends, an overall increase was observed which leads to a longer FHL of Holstein cows in Isfahan province.

Data availability
The data used for this paper were longevity data in Isfahan province, last updated on 12 January 2013 and obtained from Vahdat Industrial Agriculturists and Dairymen Cooperative.The intellectual property rights for the data belong to this company.Please contact Mehrdad Kamali via http: //vahdat-co.ir/modules/liaise/index.php?lang=english to discuss obtaining a copy of the data.

FigureFigure 1 .
Figure 1.The relationship between somatic cell and proportional culling risk

FigureFigure 2 .
Figure 2. Relationship between proportional culling risk and milk production

FigureFigure 3 .
Figure 3. Interaction between parity × stage of lactation and proportional culling risk

FigureFigure 4 .
Figure 4. Relationship between age at first calving and proportional culling risk

Figure
Figure 5. Averaged breeding values of proportional culling risk in sire

Figure 5 .
Figure 5. Averaged breeding values of proportional culling risk in sire.

Table 1 .
A summary of statistics for length of productive life of Holstein cows.

Table 2 .
Test statistics for the likelihood ratio test.

Table 3 .
Genetic parameters estimated for productive life.