Impact of conformation traits on genetic evaluation of length of productive life of holstein cattle

1 Centro Nacional de Investigación Disciplinaria en Fisiología y Mejoramiento Animal, Instituto Nacional de Investigaciones Forestales, Agrícolas y Pecuarias, Ajuchitlán, Querétaro 76280, México. 2 Programa de Maestría y Doctorado en Ciencias de la Producción y de la Salud Animal, Universidad Nacional Autónoma de México, Facultad de Estudios Superiores Cuautitlán, Cuautitlán Izcalli, Estado de México. 3 Departamento de Genética y Bioestadística, Facultad de Medicina Veterinaria y Zootecnia, Universidad Nacional Autónoma de México, Ciudad Universitaria, CP 04510, D.F., México. 4 División de Ciencias de la Vida. Campus Irapuato–Salamanca. Universidad de Guanajuato. Ex Hacienda el Copal km 9 carretera Irapuato-Silao AP 311, CP 36500 Irapuato, Guanajuato, México.


INTRODUCTION
Functional longevity was defined by Ducrocq et al. (1988) as the ability of the cow to avoid culling for reasons other than low performance and it has been reported to be strongly related with length productive life (LPL) measured as the time from the first calving to the death or culling of a cow adjusted by production level (Chirinos et al., 2007).LPL is a trait of increasing importance in cattle breeding programs.In dairy, the economic advantage of LPL lies mainly in retaining productive and healthy cows for as long as possible in the herd and it's important because when herd life is increased, expenses for raising replacement heifers can be decreased (Boettcher et al., 1997); in addition, lifetime milk production can be larger (Van Raden and Wiggans, 1995).LPL of a cow in the herd is influenced by different factors, for example fertility, milk yield, health, management and other reasons of voluntary culling (Ducrocq and Sölkner, 1998a;Weigel et al., 2003).
Many of the current genetic evaluation models for LPL in dairy cattle are based on survival analysis, which allows to combine data on both dead (uncensored/ observed LPL) and alive (censored/unobserved LPL) individuals, and enables a proper statistical treatment of censored records and accounts for nonlinear characteristics of LPL data.Survival analysis also allows the estimation of random effects using the covariance structure among observations based on genetic relationships and the calculation of animal culling risks with a mixed model (Ducrocq and Sölkner 1998b;Caraviello et al., 2004).Different traits have been used in order to increase the reliability of the prediction of breeding values for LPL in cattle.For many years, conformation traits have been used as indirect selection criteria for herd life since they can be measured early in productive life (usually during the first lactation) and have moderate genetic correlations with LPL.Different studies have used conformation traits in order to predict LPL (Vukasinovic et al., 2002;Caraviello et al., 2004;Sewalem et al., 2004).In previous studies of dairy cattle in Mexico, longevity was calculated as stability at 48 months or as LPL at the third lactation, using linear mixed models (Valencia et al., 2004) and more recently, survival analysis was used to study LPL including the effect of milk yield level, age at first calving and the time dependent variable of lactation phase (Abadía et al., 2016), but the impact of indirect indicators such as conformation traits has not been evaluated.The Mexican Holstein association scores 24 conformation traits describing udder, feet and legs, rump and body structure systems in a linear scale from 1 to 9, and the explanation of each trait is presented in Appendix A. All traits are being scored according to the standards of the World Holstein Friesian Federation (WHFF, 2005).Considering that LPL is of high economic importance, it presents low heritability and that it is measured late in life, the use of indirect predictors measured early in the LPL of a cow to improve breeding value calculations for this trait is warranted.The objective of this study was to evaluate the effect of conformation traits on the genetic evaluation of LPL of Mexican Holstein Cattle in order to provide a suitable model for the studied population.

MATERIALS AND METHODS
The data set obtained from the Mexican Holstein Association consisted of 37,870 lactation records, corresponding to 13,659 Holstein cows, calving for the first time between January 2000 and December 2014.Data files had corrected ME 305 d milk yield with an average ± STD of 10,050 ± 4,368 kg, age at first calving with an average ± STD of 24.76 ± 2.08 months, censoring indicator (34% of right censored records), LPL adjusted to 305 days per lactation, with accumulated days from the first to the fourth lactation with an average ± STD of 706.90 ± 330.82 days, information of the sire and maternal grandsire and score of 24 conformation traits: height to the withers (HW), stature(ST), size (SI), chest width (CW), body depth (BD), loin strength (LO), rump angle (RA), rump width (RW), foot angle (FA), claw uniformity (UN), heel depth (DH), bone quality (BQ), rear leg side view (RSV), rear leg rear view (RLW), fore udder attachment (FUA), front teat placement (FTP), teat length (TL) , median suspensory ligament (MSL), udder texture (TE), rear udder height (RUH), rear udder width (RUW), rear teat placement (RTP), udder depth (UD), and dairy form (DF).All conformation traits were measured in a 1 to 9 discrete scale, and classes with less than 50 observations were added to the immediate superior (classes 1 through 4) or inferior class (classes 6 through 9).Following the methodology described by Ducrocq et al. (1988) the production level was calculated based on the normal distribution of ME Milk yield, animals were classified into 10 levels of milk yield for each lactation, being 1 the lowest production level and 10 the highest.In order to better represent changes in culling risk due to within lactation reproduction and production stages of the cow, the time dependent variable lactation phase was included in the model.Three levels per lactation were considered, the first one from day 1 to 29, the second one from day 30 to 249 and the last one from day 250 to 305.The programs used for editing the data were FORTRAN 5.0 and SAS 9.3.

Model
A Weibull survival model was used for this study.The parameter estimation of the Weibull distribution and the prediction of genetic values were performed using the Survival Kit Software V3.12 (Ducrocq and Sölkner 1998a;Ducrocq 1994), using a sire-maternal grandsire model and including time dependent variables .The hazard function h(t) for a particular cow at time t was modeled as follows: where: is the probability of an animal of being culled at day t after the first calving, is baseline hazard function, assumed to follow a Weibull distribution with parameters ρ and λ, is the effect of herd-year of calving assuming that each herd has its own culling decision process within different calendar years, where i=1 is the score of the conformation traits included one at a time where m=1 to 9, is the random effect of the sire q, and is the random effect of the maternal grandsire r.Once each of the CT was included in the analysis, those which were statistically significant (p <0.05), were included in a model and in the end, two models were compared: -Model A) The base model without CT.
-Model B) The base model plus all statistically significant CT.

Heritability
This parameter was estimated in three ways for models A and B. Firstly, the log-linear scale was calculated by Ducrocq and Cassella (1996), using the following formula: Where is the sire variance, which is used to compute the variance of a log-gamma distribution.This parameter was calculated with the Digamma package of the Survival Kit because, the log-linear heritability lacks a biological interpretation and is not closely related with the reliability of genetic evaluations (Ducrocq and Cassella 1996;Ducrocq, 1999), new heritability scales have been developed.Heritability in the original scale was proposed by Ducrocq (1999), and it provides good results for the reliability of genetic evaluations when the parameter  is fixed to the value of two.Heritability on the original scale was calculated (Ducrocq 1999): Where ѵ is the Euler's constant (−0.5772), ρ is the shape parameter of the baseline Weibull distribution, and 2 log h is the heritability in the log-linear scale.This heritability method has been used successful in describing productive life in many populations (Ducrocq, 1999;Larroque and Ducrocq, 2001;Büenguer et al., 2001).

Because the aforementioned estimated
2 o h only included one of the two Weibull distribution parameters, and both parameters are strongly related, different combinations of these parameters may lead to a similar fit of the data.For this reason, Yazdi et al. (2002) developed a formula to estimate the effective heritability, which does not depend on the Weibull parameters and includes the random effects of the gamma distribution.Thus, the effective heritability ( ) was calculated as: ) 1 ) var( ( Where is the sire variance and is the trigamma function.In dairy populations, this formula has been shown to represent heritability correctly (Büenger et al., 2001;Roxström et al., 2003;Ducrocq, 2005).

Reliability
Sire genetic breeding values reliabilities (R) were calculated as a function of sire variance for models A and B as: Where is the number of daughters with observed LPL (uncensored records) and is the sire variance (Yazdi et al., 2002).

RESULTS AND DISCUSION
From each hazard function, it is possible to semiparametrically estimate a baseline survivor function (S) and if a Weibull model is adequate, a plot of ln (-ln Ŝ) versus ln (t) should give a straight line with a slope equal to  (Kleinbaum, 1996).In this study, the test graph of survival analysis, gave a straight line with a slope close to ρ (1.89) (graph not shown); this results indicates that a Weibull model is proper for the data.Similar ρ values have been estimated in other studies (Dürr et al., 1999;Vollema et al., 2000;Chirinos et al., 2007;Schneider et al., 2005).
All base model effects (AC, PL and LP) were statistically significant (p <0.005), and when CT were tested one at a time, only one body and four udder CT were statistically significant (CW, TL, MSL, TE and UD).Table 1 shows the contribution of variables to the log likelihood function, including all CT scored in Mexico.

Age at first calving (AC)
Cows calving at 21 months of age presented a relative culling risk of approximately 0.7 and it decreased as age at first calving increased until 25 months of age, when cows showed the lowest risk to be culled.Previous studies in the Mexican Holstein population also reported higher culling risks when cows calve early (Abadía et al., 2016), and it could be due to the fact that heifers that calve before 24 months of age, have not reached the appropriate size and weight to calve and this could influence the relative risks.In this study, cows calving after the 26 months showed an increased relative culling risk until the 30 months (Figure 1) possibly because they start their productive life later and reach later parities at older ages.These findings agree with those reported in other populations, where culling risk increased with age (Chirinos et al., 2007, Mészáros et al., 2008).

Production level
The influence of production level on relative culling rates is shown in Figure 2. Low producing cows are more likely to be culled than high producing cows, an indication of the influence of voluntary culling for low production.Low producing animals (level 1) have 16.5 times more possibilities to be culled at any given time than animals that have higher production levels (level 7 and 8).Similar results were found in other studies (Vukasinovic et al., 1999;Pasman and Reinhardt, 1999;Terawaki et al., 2006;Weigel et al., 2003).High risks in low production levels have been associated to voluntary culling for milk yield or to health problems, which in turn lower production (Chirinos et al., 2007, Vukasinovic et al., 2001).Animals with extremely high milk production (level 10) have a slightly higher culling risk than those with moderate production (6, 7 and 8).This increase in relative risk was explained previously in this population by Ruiz et al. (1994) and in other populations by Ducrocq et al. (1988)    and Weigel et al. (2003) who indicated that cows with high production commonly are under stress, particularly in large herds.

Lactation phase (LP)
Figure 3 shows the effect of relative risk in each lactation phase (including 1st, 2nd, 3rd and 4th lactations).
Relative culling rates increased from the beginning to the end of each lactation.This trend is in accordance with results presented by Dürr, et al. (1999); Vukasinovic et al. (1999) andRoxtrom et al. (2003) who explained the changes in relative risks as a result of selection pressure changes during the lactation or more intensive culling of non-pregnant cows near the dry period, because cows are culled when they are not pregnant and milk production is finished (Dürr et al., 1999;Vukasinovic et al., 1999 andTerawaki et al., 2006).The relative culling risks across lactations also increased with the cow´s age as reported in other studies in the Mexican dairy cattle Holstein population (Abadía et al., 2016) and in other populations (Dürr et al., 1999;Ducrocq, 1999;Terawaki et al., 2006;Chirinos et al., 2007).

Conformation traits (CT)
Five CT were statistically significant as predictors of LPL.One of them is from the capacity and structure system (CW) and the others from the mammary system (TL, MSL, TE and UD).These findings agree with the results of Dadpasand et al. (2008), which concluded that mammary system traits have a strong relationship with functional longevity.The relative culling risks within scores of CW, TL, MSL, TE and UD traits are presented in Table 2. CW low scores show a high relative risk of culling which decreased as CW scores increased (Table 2, CW).This is probably because cows with low scores, anatomically have narrow chests, which could indicate not enough space for housing vital organs, as the heart and lungs whereas high CW scores indicate a wide chest and enough thoracic capacity.Similar CW culling risk patterns were reported in the Canadian Holstein Population, although, this trait had a low contribution to the likelihood function (<5%; Sewalem et al., 2004) and other populations where other authors have reported moderate correlation between CW and LPL (-0.24) in dairy cattle (Zavadilová et al., 2009).TL presented the lowest culling relative risk in intermediate classes (4, 5 and 6) (Table 2, TL), probably because it is difficult to attach milking units on short teats (scores 1 to 3) and the vacuum will not be optimal, whereas long teats (score 7 to 9) could be in contact with the feet and legs and might catch an infection.Similarly to CW, TL has been reported to have a low contribution to the likelihood function (Sewalem et al., 2004;Morek-Kopeć and Zarnecki, 2012) and its genetic correlation with LPL was still lower (-0.16)than CW (Zavadilová et al., 2009).MSL low scores (≤ 4) did not present a clear relative risk of culling trend, but for the upper five classes (5 to 9) it was observed that animals with high scores, have lower risks to be culled (Table 2, MSL) because animals with strong MSL present better supported udders.In other populations, MLS has been reported with medium contributions (~17%) to the likelihood function (Sewalem et al., 2004).TE is a qualitative CT highly associated to longevity in Holstein cattle (Sewalem et al., 2004).In the present study, TE showed that animals with the highest score have lower culling risks than animals with low scores (Table 2, TE).UD is a CT strongly related with true and functional longevity in many populations (Caraviello et al., 2004;Morek-Kopeć and Zarnecki, 2012;Sewalem et al., 2004;Zavadilová et al., 2009).As it was found in Canadian and German Holsteins UD presented a medium optimum in the Mexican Holstein cattle which suggests that cows with the udder floor slightly above the hock have less risk to be culled than cows with udder floors below or above the hock (Sewalem et al., 2004) (Büenger et al., 2001),.These results agree with those reported in Polish Holstein cattle (Morek-Kopeć and Zarnecki, 2012) for true longevity.

Heritability
As mentioned earlier, this parameter was calculated in different ways.LPL heritabilities calculated with model A were 0.06, 0.10 and 0.09 for logarithmic, original and effective scales respectively while for model B the values were 0.08, 0.12 and 0.14, respectively.Previous studies used model A for predicting longevity in the Mexican Holstein population using more life time records (36,507) because animals included in that study were not limited by the presence of conformation traits but only for milk yield.The estimated logarithmic, original and effective scale heritabilities were higher than the ones presented with a similar model (0.08, 0.13 and 0.12, respectively (Abadía et al., 2016).However, results of the current study are in the range of values reported for LPL in other populations (from 0.02 to 0.11 for logarithmic scale from 0.04 to 0.22 for original scale and from 0.03 to 0.19 for equivalent scale (Ducrocq, 1999;Vucasinovic et al., 2001;M'hamdi et al., 2010;Wiebelitz et al., 2014).When the statistically significant CT were included in the model (Model B) for LPL, heritability estimation improved for all scales and the heritability in the effective scale was higher compared to values obtained by Abadía et al. (2016).The heritability values obtained in the current study were also higher values than values obtained in Australian Brown Swiss and Simmental cattle (Sölkner et al., 1999).

Estimated breeding values and reliability
The mean and standard deviation of BV, expressed, as relative risk ratios were 1.06 ± 2.89 and 1.31 ± 4.22 for models A and B respectively.BV distributions for both models are showed in Figure 4. Results suggest that model B, which includes CT has a better fit for LPL, because the included CT increased the proportion of the explained genetic variance and the somewhat larger BV range allows to select the better animals to improve LPL.Additionally, Model B increased the average reliability by 14 percentage points, compared to model A, because reliability depends on the estimated sire variance (Yazdi et al., 2002) which was 0.03 and 0.04 for models A and B, respectively.Results of this study agree with the findings of Vukasinovic et al. ( 2002), which concluded that the use of CT improves reliability of longevity.

Conclusion
Survival analysis is adequate to calculate length of productive life, to estimate heritability and to predict breeding values in the Mexican Holstein cattle registered population.
The time dependent variables included in the analysis were good predictors for length of productive life, and five conformation traits were statistically significant in order to improve the length of productive life model.One conformation trait (chest width) was from the structure and capacity system, and the others were related with udder composition (teat length, median suspensory ligament, udder texture and udder depth).The inclusion of these five conformation traits improves the estimation of length of productive life, the prediction of breeding values and its reliability.The only challenge of using conformation traits in the longevity prediction for the Mexican Holstein population is that this information (conformation traits) must be available for all animals included in the analysis, which limits, especially, the number of uncensored data.Nevertheless, the reliability gain and the improvement of breeding value estimation justify the inclusion of the conformation traits.Therefore, inclusion of conformation traits in the length of productive life model of the Mexican Holstein cattle is recommended.
to 911, is the effect of age at first calving in months, where j *Corresponding author.E-mail: ruiz.felipe@inifap.gob.mx .Tel: +52 55 38718700 ext 80216.Author(s) agree that this article remain permanently open access under the terms of the Creative Commons Attribution License 4.0 International License 14 classes from 17 to 30 months, is the effect of production level adjusted in each lactation, where k=1 to 10, and is the effect of lactation phase with changes at the 29, 249 and 305 days in each lactation, where l=1 to 12, three phases for each of the first four lactations.

Figure 1 .
Figure 1.Relative culling risks associated to age at first calving of Holstein cattle in Mexico.

Figure 2 .
Figure 2. Relative culling risks associated to milk production levels of Holstein cattle in Mexico.

Figure 4 .
Figure 4. Length of Productive Life (LPL) breeding values distribution for the two models evaluated: Model A) without conformation traits and Model B) including five conformation traits (Chest width, Teat length, Median suspensory ligament, and udder texture and depth.

Table 1 .
Variable contributions to the likelihood function and Chi square test.

Table 2 .
Relative culling rate of the Mexican Holstein cows for the statistically significant (p <0.05) conformation traits.TL=Teat length, MSL= Median suspensory ligament, TE= Udder texture, UD=Udder depth.Classes with less than 50 observations were not included in the analysis.