Genetic evaluation for longevity of Croatian Simmental bulls using a piecewise Weibull model

Genetic evaluation of sires for functional longevity was conducted using survival analysis techniques. The data set consisted of 49 659 Simmental cows with first calving from 1997 to 2008. A piecewise Weibull sire model was used to estimate breeding values of 251 bulls for functional length of productive life of their daughters. The model was stratified by parity i.e. a separate baseline hazard was computed for each stratum. Besides the random sire effect, the model included the fixed time independent effects of age at first calving, herd size and region as well as the time dependent effects of relative milk production and year*season of first calving. The highest impact on longevity was found for relative milk production. Cows with the lowest milk yields were at approximately 2.7 times higher risk of culling compared to cows with average milk production. Effects of age at first calving, herd size and country region had lower impact on longevity. Sire variance was 0.023 which results in a heritability of 0.06 for functional length of productive life. The average approximate reliability of estimated breeding values was 0.49. Genetic trend showed no clear tendency by year of birth of bulls.


Introduction
Longevity is the most important functional trait in the majority of selection indices all over the world.In dairy cows it is more relevant to speak about so called functional longevity (Ducrocq 1987).Functional longevity represents the ability of a cow to avoid involuntary culling due to low production.It is an important trait to measure the general fitness of a cow, particularly when direct measures of the functional traits (health and reproductive performance) are not available.In Croatia, conformation traits and somatic cell count are the only non production traits recorded at national level, making the genetic evaluation of functional longevity very important in order to improve functional characteristics and animal welfare of cows.
Increased longevity in dairy cows affects the average herd production due to an increase in the number of mature cows, which produce more milk than young cows.Also, it influences economics of the farm by reducing expenses for replacement costs and veterinary treatments.Furthermore, a higher herd life gives the potential to increase selection intensity on milk production.In general, the relative economic value of longevity compared to milk production largely depends on the actual level of herdlife in the respective population (Essl 1998).According to Burnside et al. (1984) and Allaire & Gibson (1992), the range of the economic value of longevity is between 25 to 70 % of the value for production.For Austrian conditions, Essl (1982aEssl ( , 1982b) calculated that a decrease of average herd life from 5 to 4, 4 to 3 and 3 to 2 lactations is economically equivalent to a decline of 263, 572 and 2 246 kg fat corrected milk per lactation, respectively.According to other studies, the relative economic value of one additive genetic standard deviation of herd life compared to that of milk yield (set to 1) was within range of 0.25 to 1.01 with a rough average of 0.5 (Veerkamp et al. 1995, Van Raden & Wiggans 1995, Miesenberger 1997).The relative values of yield and herd life derived from several studies indicate that yield is 1 to 8 times as important as herd life (Van Arendonk 1991, Dekkers 1993, Harris & Freeman 1993).In spite of this rather wide range of the relative economic weight, herd life is the most important functional trait in all studies.
For longevity analysis, Weibull proportional hazard models (Ducrocq & Sölkner 1998) are used by most of the authors, and these provide heritability estimates for longevity from 0.04 to 0.21 (Ducrocq 1999, Schneider & Miglior 1999, Vollema et al. 2000, Páchová et al. 2005, Sewalem et al. 2005, Mészáros et al. 2008, Bonetti et al. 2009).Although this approach is considered as the most appropriate way to analyse longevity data, certain problems related to the existence of censored records still remain.When used for genetic evaluation and estimation of breeding values of sires, the predictive power of survival analysis directly depends on the proportion of censored records.An increased proportion of censored records results in lower reliabilities particularly for young sires (Vukasinovic 1999).
In Croatia, genetic evaluation of longevity in dairy cattle is still not considered in breeding programmes.Raguž et al. (2011) determined factors affecting longevity in Croatian Simmental and Holstein cows.The time dependent effects as the level of production within herd and parity had the highest impact on the functional length of productive life.
The aim of this study was to analyse effects affecting length of productive life as well as to estimate genetic parameters and breeding values for functional longevity in Croatian Simmental cows.

Material and methods
Longevity data and production records were provided by the Croatian Agricultural Agency and collected on 11 December 2009.After editing, records from 49 659 Simmental cows with first calving from January 1997 through May 2008 remained for analysis.All cows were required to have valid sire, paternal grandsire and granddam identification as well as age at first calving between 20 and 40 months.The total number of herds over the whole period, assigned to 4 geographical regions, was 10 754 comprising herds with only 1 to maximum 66 cows due to deleting the animals with unknown sire pedigree.The minimum number of daughters per sire was 10 with a total number of sires of 251.
The average length of productive life was 3.10 years for all animals, 2.95 for the uncensored records and 3.55 years for the censored records.Average age at first calving was 27 months.Milk production in first lactation was in average 3 891 kg with an SD of ±1 036 kg while the average lifetime production was 10 186 kg with an SD of ±6 254 kg (Table 1).Longevity was defined as the number of days from first calving until culling or censoring.
For the animals with missing culling date, the date of the last known lactation was used.In case of missing information on culling, the cow was considered as culled if the number of days between the end of last lactation and the date of data collection exceeded 365 days.Otherwise, the record was considered as right censored.If the cow reached more than 7 lactations, it was considered as censored at the end of 7th lactation because of low number of animals in higher lactations.In total, the censoring procedure resulted in 11 559 (23.28 %) right censored cows.
The following proportional hazards model was used: where λ(t) is the hazard function (instantaneous probability of culling) for a given cow at time t; λ 0 (t) is the Weibull baseline hazard function with scale parameter λ and shape parameter ρ; ys i is the random time dependent effect of year*season following a log-gamma distribution; rp j is the fixed time dependent effect of the relative milk production; afc k is the fixed time independent effect of age at first calving; reg l is the fixed time independent effect of region; hs m is the fixed time independent effect of herd size; s n is the random time independent effect of sire.The model was stratified by parity, reseting the baseline to 0 for each stratum.
For estimation of the standard lactation milk yield, the Test interval method approved by ICAR ( 2009) was used.The relative milk production in each lactation was calculated after their adjustment to the first lactation, as the milk production of a cow was compared to the herd average in the given year.The resulting difference was expressed in standard deviations and divided into 9 classes.The details about calculations of relative production were described in Mészáros et al. (2008).To account for culling for production, the effect of within herd production level was included, which resulted into estimation of the functional productive life.The lower and upper bounds for these classes are given in Table 2.
The effect of the herd size was taken into account using 4 classes.The highest proportion of the cows was from small herds less than 6 animals.Other animals were divided into classes ≥6 and <10, ≥10 and <15 and ≥15 as shown in Table 4.
The cows were grouped according to the slightly different climatic conditions in four groups covering the four Croatian regions: Slavonia, North-West Croatia, Upland Croatia and Mediterranean Croatia (Table 5).According to the age of first calving, cows were divided in 7 different classes (Table 3).The year*season effect was taken into account through three seasons, where the periods January to April, May to September and October to December represented seasons 1 to 3, respectively.
The effect of the parity was integrated in the model as a stratification variable where a different hazard line is estimated for each parity, and thus no results for risk ratios where obtained for this effect.
The heritability was calculated as proposed by Yazdi et al. (2002): where h 2 is the coefficient of heritability for functional length of productive life; σ s 2 is the genetic variance of sires and p is the proportion of uncensored records.
Breeding values for functional length of productive life were computed using these parameter estimates.Average breeding values were calculated compared to the mean and standard deviation of the group of base sires, with birth year between 1995 and 1999.The following equation for breeding value computation was used: where BV is the breeding value of a sire; estimate is the regression coefficient estimate for a particular bull; a is the mean of estimates in the group of base sires; sd is the standard deviation of estimates of base sires.The minus sign represents the alternation of breeding values to express longer productive life with higher values.Reliability values for each bull were computed according to Ducrocq (1999): where R is the value of the reliability; n is the number of uncensored daughters of a bull; h 2 is the estimated heritability of functional length of productive life.All analyses were conducted using the software Survival Kit v. 6.0 (Ducrocq et al. 2010), while the initial textual file was prepared using SAS v. 9.1 (SAS Institute Inc., Cary, NC, USA).

Results
All effects in the model were found to be highly significant (P<0.01).The Weibull parameter ρ was estimated separately for each parity due to using piecewise Weibull model and accounting parity as a stratification variable in the model.The lowest ρ was estimated for the 6th parity (2.96) and the highest (4.56) for the 1st parity.It means that the increase of the culling risk is the highest for animals in the 1st parity.The results for fixed effects and the random effect of the year*season were all presented as relative culling risks, defined as the ratio between estimated risk of being culled under the influence of certain environmental factor [exp(β i ) for level i] and the reference class with risk ratio set to 1 (β i =0).Values of risk ratio >1 indicate higher culling risk associated with that effect in comparison with the reference class and vice versa.For example, if the relative culling risk for a given class is 4, it means that a cow in that class has four times higher chance to be culled compared with a cow in the reference class (set to 1) for a given effect.So if the relative culling risk for a given class is 0.3, then a cow in that particular class has a 70 % less chance of being culled then a cow in the reference class.The culling risks for relative milk production are presented in Table 2.The cows with the highest amount of milk production were at lowest risk to be culled.In contrary, the risk continuously increased with decreasing relative milk production.Cows with a milk yield ±0.2 standard deviations relative to the herd mean were assigned as a reference class.In comparison with cows in the reference class, the cows with the highest milk yields (more than 1.5 standard deviations aboveherd mean) were at approximately 2.4 times lower risk to be culled.Cows with milk yield 1.5 standard deviations below the herd mean were 2.7 times more likely to be culled than cows with average milk yields.Age at first calving had weak effect on culling risks (Table 3).The highest risk ratio was found for the cows that first calved at approximately 24 months of age.Generally speaking, cows having calved for the first time at lower age were more exposed to be culled.Later calvers were at slightly lower risks of culling.Table 4 shows the results for the effect of the herd size.The highest selection pressure was perceived in bigger herds (>15 animals) resulting ina risk ratio of 1.54.The lowest risk of culling was found for the class of the herd size between 6 and 10 animals which is also assigned as a reference class.The effect of region affected longevity in a lower degree (Table 5).Cows from Eastern, Upland and Mediterranean Croatia had slightly lower chance to be culled compared to North-West region.Our statistical model for the analysis of functional longevity resulted insire variance of 0.023, which is a heritability of 0.06. Figure 2 shows the distribution of estimated breeding values of the sires calculated using equation ( 3).
The values ranged from 67 to 165.Sires born between 1995 and 1999 were taken as the base population with an average breeding value of 100.
The genetic trend for functional longevity comprising reliability values for each year is presented in Figure 3.The minimum number of sires born per year was 5.
Reliability of breeding values of sires was calculated according to equation ( 4) with an average reliability of 0.49 and a standard deviation of 0.28.The average reliability was achieved from bulls with 63 uncensored daughters.To achieve higher reliabilities as 60, 70, 80 and 90 %, 97, 153, 264 and 658 uncensored daughters were needed, respectively.Maximum reliability of 0.95 was found for a sire with 1 286 uncensored daughters.

Discussion
Our model contained seven effects including time-dependent (milk production within herd, parity, year*season) and time-independent effects (age at first calving, herd size, region and sire).Many authors reported milk production as one of the most important effects affecting longevity (Vollema et al. 2000, Vukasinovic et al. 2001, Canavesi et al. 2004, Sewalem et al. 2005, Mészáros et al. 2008, Bonetti et al. 2009, Raguž et al. 2011).According to our data structure, milk production was taken into account through 9 classes ranging from <−1.5 to >+1.5 SD relative to herd mean.However, this effect was considered through different number of classes according to different authors ranging from only 3 (Sewalem et al. 2005) to 9 (Ducrocq et al. 1988, Canavesi et al. 2004).Our results report a clear tendency in risk ratios going from the highest value (2.7) for animals with low milk yields to lowest value (0.4) for animals with highest production.That kind of trend is obvious in most of other studies (Páchová et al. 2005, Chirinos et al. 2007, Mészáros et al. 2008, Bonetti et al. 2009) and corresponds to our results.Vukasinovic et al. (2001) found that cows producing less than 80 % of the herd average are at 3 to 4 times higher risk than their herd-mates with average production.High producing cows are less likely to be culled.Páchová et al. (2005) and Mészáros et al. (2008) found similar results for high producing animals being at approximately 5 times higher risk to be culled than the herd average.Similar findings were reported by Vollema & Groen (1998), Weigel et al. (2003) and Sewalem et al. (2005).In many studies, fat and protein yield was investigated showing the same trend, but with lower range of culling risk compared to milk production (Vukasinovic et al. 2001, Egger-Danner et al. 2005, Sewalem et al. 2005, Chirinos et al. 2007).
Age at first calving is generally considered as an effect of minor importance for functional longevity.The risk of culling by age at first calving was almost unchanged for all classes except class 2 for cows calving for the first time at 24 months.This class had the highest risk ratio value of 1.12.Many authors found this effect as a significant one, showing a slight increase in culling risks towards later calvers (Dürr et al. 1999, Vollema et al. 2000, Caraviello et al. 2004, Páchová et al. 2005).It could be speculated that the higher age at first calving is a signal of low fertility or other health problems of a cow.In our case, younger cows were at higher risk of culling probably due to insufficient development for the first calving in respective age.The risks stabilize around the reference and later classes.Ducrocq (1994) and Ojango et al. (2005) reported no significant effect of the age at first calving.
In the Croatian population of Simmental cows, the highest proportion of animals belongs to small herds, that is herds containing <6 animals.From this reason, effect of the herd size change could not be taken into account in contrary to other studies (Dürr et al. 1999, Vollema et al. 2000, Forabosco et al. 2006, Chirinos et al. 2007).Only the effect of herd size was considered.The results were expected showing higher culling risk for animals from herds bigger than 15 animals.Probably, the selection pressure in such herds is much higher compared to small and less intensive family farms with only few animals.
Due to different climatic condition, Croatia is geographically divided into 3 main regions.For the purpose of our study, one more region (North-West Croatia) was added for a special reason.This region has the most developed dairy farming in Croatia and thus it should be considered separately.The results confirmed our assumption showing that animals from the mentioned region were at higher culling risk compared to other regions.On theother hand, cows from the Mediterranean region were at lowest risk of culling which is also expected due to less intensive farming and smaller herd size compared to other regions.Several studies investigated the effect of country region (Smith et al. 2000, Ducrocq 2005, Chirinos et al. 2007, Raguž et al. 2011) where Chirinos et al. (2007) obtained the similar consistency of the results among the regions as in our case, concluding that the proposed model could serve as a basis for the development of the model used at the national level.In contrast, Smith et al. (2000) found different percentages of cows leaving the herds in different regions with providing the reasons of disposals.Ducrocq (2005) found systematic differences between regions in France, although the magnitude of the effect was not extremely large.
In contrast to results of some previous studies where lactation number and stage of lactation were estimated jointly (Ducrocq 1994, Dürr et al. 1999, Vukasinovic et al. 2001, Roxström & Strandberg 2002, Mészáros 2008), in our research only the effect of parity was taken into account as a stratification variable because of the dataset specificity.In fact, there was no cows in earlier lactation stages but only cows over 150 days in lactation.Due to this reason, effect of the lactation stage couldn't be taken into account.
In case of taking the parity as a stratification variable using a piecewise Weibull model, a different hazard line wasestimated for each of the 7 parities.The approach using piecewise Weibull hazard function with stratification within year and parity was applied by Ducrocq (2005) and Van der Linde (2004) in an improved model for the French and Dutch genetic evaluation of dairy bulls on length of productive life of their daughters.
A slightly decreasing trend was observed for the effect of the year*season.The lowest risk ratios were found in years near the end of the study which could be explained by higher proportion of censored animals in those years.In the last few years an unclear increase of risk ratios was determined.Namely, it could be explained by a higher proportion of animals that ended with first lactation in 2009 but did not start with the second lactation and thus were culled.It resulted in an increase of culling risks in a particular year.Raguž et al. (2011) found an increasing tendency in risk ratios for this effect using different models.This trend could be partially explained by more intensive selection in recent years, also due to the fewer number of observations in first few years of the research.
To estimate breeding values, a base group of sires born between 1995 and 1999 was defined.For comparison to all other bulls, the average and standard deviation of these sires were used.Younger bulls had lower average breeding values than the older ones probably due to higher number of still living daughters.In genetic trend from 1986 to 2001 no tendency was observed due to lack of any selection for longevity.Distribution of breeding value reliabilities over the sires' birth years indicates clear trend.Lower values from later years are as expected due to higher number of censored daughters.Sires born from 1986 to 1988 had lower reliability values probably due to social events in early 1990's in Croatia and eventually decreasing number of daughters per sire.A peak value from 1993 could be the consequence of considerable import of animals from western countries resulting with large increase in a number of daughters.
The reliability of breeding values depends on the number of uncensored records and is very low for young sires having only few uncensored daughters.Increasing the reliability of breeding values for young sires will be one of the most important tasks for the population of Simmental bulls in Croatia in the future.In a case when a young sire has no uncensored daughters, the reliability of his breeding value is based on pedigree information only (Vukasinovic et al. 2001).The reliability could be increased by taking weighted information on censored daughters into account as well as by the inclusion of information from early predictors (type traits, somatic cell count and fertility).Posavi et al. (1999) found high heritability estimates of 0.50, 0.43 and 0.42 for teat placement, length of rump and rear legs side view respectively in Croatian Simmental cows.Moderate and lower heritability estimates were found for back length, rump height, rump width, muscularity and some of the udder traits as suspensory ligament, teat length and teat thickness.In an analysis of the relationships between type traits and longevity using survival analysis, Jovanovac & Raguž et al. (2011) found high linkage between 12 type traits and true longevity.It indicates that type traits could be used in further studies with aim to improve early predictions of longevity breeding values, particularly for young sires with higher number of censored daughters.Inclusion of type traits using MACE techniques (Schaeffer 1994) could also be promising and requires further research.
In conclusion, a piecewise Weibull proportional hazard model was used to evaluate the impact of several factors on the culling risk and to estimate the heritability and breeding values of sires for functional longevity.Effect of parity was included as a stratification variable where different baseline hazard was estimated for each parity.The factors with highest impact on the length of productive life were milk production within the herd and year*season.Cows with the lowest milk production level within herd had 2.7 times higher chance to be culled than animals with average milk production.On the other hand, animals with the highest production level were at approximately 2.4 times lower risk to be culled in compare to cows with average milk yields.The effect of year*season showed respectable impact on cows longevity with a slightly decreasing trend until 2006 followed by unclear tendency in the last few years.Considering the effect of herd size, animals from bigger herds were at relatively higher risks of culling (1.5) than the animals from smaller herds.Effects of age at first calving and country region were effects of minor importance for longevity.A sire variance of 0.023 and heritability estimate of 0.06 were calculated for the data set considered.Average reliability of the breeding values was 0.49 with the maximum value of 0.95.In genetic trend, no clear tendency by year of birth of bulls was observed.Further analyses of the length of productive life should be conducted with the aim of improving the reliabilities of breeding values of young sires.

Figure 1
Figure 1 presents the culling risks for the random time dependent effect of the year*season of calving.A slightly decreasing trend between 1999 and 2006 was determined followed by unclear tendency in the last few years.Overall, in 1998's and 2007's season 3 the highest risk ratios were found while in 2008's season 2 and 3 the lowest risks of culling were found.

Figure 1 :
Figure 1: Estimates of risk ratios for the year*season effect Figure 2 Distribution of es timated bulls' breeding values

Table 2
Estimates of risk ratios for classes of the relative milk yield 1 expressed as difference in standard deviations compared to herd average

Table 3
Estimates of risk ratios for classes of age at first calving

Table 4
Estimates of risk ratios for classes of the herd size

Table 5
Estimates of risk ratios for classes of the country region