Document Type: Research Articles
Authors
^{1} Department of Animal Science, Faculty of Agriculture, University of Jiroft, Jiroft, Iran
^{2} Department of Animal Science, University College of Agriculture and Natural Resources, University of Tehran, Karaj, Iran
^{3} Department of Animal Science, University of Wisconsin-Madison, Madison, WI 53706, USA
Abstract
Keywords
INTRODUCTION
Structural equation models (SEMs) allow studying the causal relationships between phenotypes (Wright, 1934). These models were first introduced in genetics by Wright (1921). Calving is an important event on dairy farms, calving complications result in potential loss and/or impaired subsequent production and reproduction performance with implications on animal welfare, increased labor and veterinary costs with associated negative impact on revenue (Eaglen et al. 2012). Calving difficulty (CD) or dystocia is the consequence of inconsistency between the size of the calf and pelvic area of the dam and may be characterized as prolonged gestation and calving (Meijering, 1984). Phenotypically, CD influences on the lactation performance of cows in the subsequent lactation period (Dematawewa and Berger, 1997; Lopez de Maturana et al. 2007b; Atashi et al. 2012a). Fitting CD, solely as a fixed effect, would ignore the fact that both CD and lactation related traits may be genetically correlated. Then again, only considering the existence of genetic correlation between these traits via a standard multi-trait model would ignore the cause-and-effect relationships between them. Mokhtari et al. (2016) studied the phenotypic and genetic relationship between CD and fertility traits in first-parity Iranian Holsteins under standard and recursive models. Model comparison revealed that considering causal effects from CD on the studied fertility traits is of importance in terms of genetic evaluation of the considered fertility traits. Atashi et al. (2012b) reported that the incidence of calving difficulty in the first parity Iranian Holsteins (16%) is higher than other parities (8%). Therefore, the existence of potential cause-and-effect influence from CD on the subsequent lactation-related traits in first-parity Iranian Holsteins necessitates considering such causal effects for estimating more accurate genetic parameters of the lactation related traits in first-parity of Iranian Holsteins. To our knowledge, there has been no previous study for inferring causal effects of CD on the subsequent lactation curve traits in first-parity Holsteins. Although the influence of CD on the lactation performance in Holstein dairy cows have been discussed (Atashi et al. 2012a) but such influence has not been studied considering possible causal effects from CD on the lactation curve traits. Therefore, in the present research causal effects of CD on the lactation curve traits were studied in first-parity Iranian Holsteins. Furthermore, phenotypic and genetic parameters of CD and the lactation curve traits were estimated and compared under recursive mixed models (RMMs) and standard mixed models (SMMs).
MATERIALS AND METHODS
Data edition and studied traits
The milk yield records initially analyzed included 989582 test-day records of 160243 cows collected from 131 herds of Iranian Holstein dairy cows from 1995 to 2014 by the Animal Breeding and Improvement Centre of Iran.Test-day records with first milk recording measured between 5 and 60 days of lactation period, consecutive sampling intervals of 25-35 days and with not greater than 305 dayslactation length were considered. Wood function is considered for describing lactation curve in the studied data set of first-parity Iranian Holsteins in the present study. Therefore, edited test-day milk yield records were used for calculating the lactation curve traits applying Wood function as below (Wood, 1967):
y_{t}= A_{p} t^{Bp} e^{-Cp t}
Where:
t and y_{t}: represent day in milking and daily milk yield on day t of lactation period.
Ap: initial milk yield.
Bp and Cp: ascending and descending slope of the lactation curves, respectively.
Cows with atypical lactation curves i.e. with negative parameters Bpor Cp were discarded. Fitting wood function on the test day milk yield records was performed applying NLIN procedure and the Newton-Gauss method (SAS, 2004). Then individual cow lactation curve parameters including peak milk yield (Y_{m}), days to attain peak yield (T_{m}) and milk persistency (Pers) calculated using the following formula:
Only single records from the artificial inseminations were kept, and all suspect records, including records with out-of-range values or records with missing information were removed. Age of the cow at the first calving was limited to the range from 20 to 38 months. Detailed information on the considered editing protocol for CD records have been presented by Mokhtari et al. (2016). Finally, 22872 cows with records for both CD and the lactation traits were considered for analyses.
Statistical analyses
For addressing the significant fixed effects to be fitted in the final models least squares analyses were carried out using the general linear method (GLM) procedure of SAS software (SAS, 2004). Least squares means for all the studied lactation curve traits across the easy calving and difficult calving cows were compared applying t-test and a significant level of P-value equal or lower than 0.05 was considered. Because of the extreme category problems, for categorical traits such as CD genetic analysis under threshold animal models with direct additive and maternal additive genetic effects via the Gibbs sampling may result in biased estimates, poor mixing and "blowing up" of the Gibbs chains (Hoeschele and Tier, 1995). When a Gibbs chain blows up, the genetic variance will increase and soon reaches irrational amount, in this situation the inverse of the genetic variance matrix diminish to zero and Gibbs sampling stops (Luo et al. 2001). Therefore, most applications of threshold models for genetic evaluation of animals for categorical maternally influenced traits, such as CD in the present study, were based on sire-maternal grand sire models (Eaglen et al. 2013). In a sire-maternal grand sire model the effects of sire of the calf and sire of the dam are considered as additive direct and maternal genetic effects, respectively (Eaglen et al. 2012). The general form of the bivariate models considered for genetic analysis in matrix notation is as follow:
Λy_{i}= Xb_{i} + Z_{h} h_{i} + Z_{s1} s_{1} + Z_{s2} s_{2i} + e_{i}
Where:
y_{i}: imply vector of the records for i^{th} trait (the liability for CD and calculated records for the lactation curve traits).
b_{i}: vector of fixed effects for i^{th} trait containing sex of calves in 2 levels (for CD only), calving month in 12 levels and age of the cow at the first calving in 19 levels (20 to 38 mo).
h_{i}: vector of herd-year-season of calving in 769 levels.
s_{1}: vector of sires of the calves effects in 933 levels (fitted only for CD as sire effects).
s_{2i}: vector of sires of the cows effects for i^{th} trait in 563 levels (fitted for CD as maternal grand sire effects and for the lactation curve traits as sire effects).
e_{i}: vector of residual effects for i^{th} trait. Four calving seasons were defined: April to June, July to September, October to December and January to March.
X, Z_{h}, Z_{s1} and Z_{s2}: incidence matrices which associating the corresponding effects to vector of records.
The matrix of structural coefficients (Λ) is a square one; off-diagonal elements were determined according to causal structures between the considered traits (Valente et al. 2010). In this case, the matrix of Λ is a 2 × 2 matrix as below:
Implying that a recursive effect (represented by λ_{21}) from CD on each of the subsequent lactation curve traits. In the above model, if the structural coefficient matrix to be considered as an identity matrix structural equation model becomes the standard model. In other words, RMMs are multivariate models whose equations carry one-way causal links among the involved traits and the causal information carried by RMMs allows understanding how the value of one trait is affected (and not only associated with) the value of the other traits. Performing a multivariate analysis taking all traits into account was not feasible because of convergence not achieved in the present study. Therefore, six bivariate analyses including CD-Ap, CD-Bp, CD-Cp, CD Tm, CD-Ym and CD-Pers were performed under both SMMs and RMMs. Lopez de Maturana et al. (2007a) pointed out that the method described by Gianola and Sorensen (2004) is not straightforward enough to perform in a general manner and showed that recursive models could be handled by fitting parent trait (in this case, CD) as a covariate for the other trait(s) while genetic correlations between traits are considered in multivariate analyses. Therefore, this methodology was applied in the present study. Unabridged information about the methodology is given by Lopez de Maturana et al. (2007a).
Statistical inference
Bayesian MCMC implementation was performed applying the THRGIBBS1F90 program of Misztal et al. (2002), which implements Gibbs sampling to evaluate the posterior density of the parameter estimates. The length of the chain and the burn-in period were examined by visual inspection of the trace plots of posterior samples of the parameters. For each model, 300000 iterations were run and the posterior samples from each chain were thinned considering thinning intervals of 50 iterations after discarding the first 50000 iterations as burn-in period. Hence, 5000 samples were considered for computing features of the posterior distribution. Posterior analyses for calculating posterior means and posterior standard deviations were carried out applying the POSTGIBBSF90 program of Misztal et al. (2002). It was assumed that the effects of sires of the calves and sires of the cows followed multivariate normal distributions, a priori, with a null mean vector and a (co)variance matrix G A, where G and A are the genetic (co)variance matrix and numerator relationship matrix among bulls, respectively. A multivariate normal distribution was assumed for the effects of herd-year-season of calving with null mean vector and (co)variance matrix H I_{h}, where H is a (co)variance matrix for the herd-year-season of calving effects. The I_{h} matrix is an identity matrix of order equal to the number of herd-year-season of calving (769). Furthermore, it was assumed that the vector of residual effects followed multivariate normal distributions with a null mean vector and (co)variance matrix R I_{n}, where I_{n} is an identity matrix of order 22,872 and R is the residual (co)variance matrix; shows the Kronecker product. A multivariate normal distribution was also assumed for the fixed effects. The RMMs are not identifiable at the likelihood level due to the presence of extra parameters including structural coefficients. To achieve identification, it was assumed that residual correlations in the system were uncorrelated. In other words, R was assumed to be a diagonal matrix in RMMs for the identification purposes. An inverted Wishart distribution was assumed for the prior distributions of the genetic, herd-year-season of calving and residual effects, so that their fully conditional posterior distributions were also inverted Wishart. CD was analyzed applying threshold models, having 3 categories. For identification purposes, the first and second thresholds for CD were set to 0 and 1, respectively. The SMMs and RMMs were compared using deviance information criterion (DIC), the DIC considers the trade-off between model goodness-of-fit and corresponding complexity of model into account (Bouwman et al. 2014). Models with smaller DIC values are better supported by the data.
Direct and maternal genetic effects
The variances related to effects of sires of the calves and sires of the cows were transformed to the direct additive and maternal additive genetic (co)variances (Willham, 1972) as:
Where:
σ^{2}_{d}, σ^{2}_{m}, σ^{2}_{s1} and σ^{2}_{s2}: direct additive genetic, maternal additive genetic, sires of the calves and sires of the cows variances, respectively.
σ^{2}_{dm }and σ^{2}_{s1s2}: covariance between direct additive and maternal additive genetic effects and between sires of the calves and sires of the cows effects, respectively.
Additive genetic covariances between direct and maternal effects among each pair traits of i and j were computed according to Kriese et al. (1991):
Where:
: covariance between direct additive genetic effects for traits i and j.
: covariance between maternal additive genetic effects for traits i and j.
: covariance between direct additive genetic effects of trait i and maternal additive genetic effects of trait j.
: covariance between maternal additive genetic effects of trait i and direct additive genetic effects of trait j.
The phenotypic variance (σ^{2}_{p}) and direct heritability (h^{2}_{d}) were calculated as follow:
In which, σ^{2}_{hys} and σ^{2}_{e} denote the variances due to herd-year-season of calving and residual effects, respectively.
System parameters
The interpretation of parameters obtained under RMMs, the so-called system parameters in SEMs literature, is different from that of the analogous ones obtained under SMMs (Gianola and Sorensen, 2004). Therefore, further transformation is required to be able to compare (co)dispersion of random effects among two models fitted (i.e. RMMs and SMMs). Transformations for the estimated (co)variance matrices to the standard mixed model scale were carried out as:
G^{*}= Λ^{-1} G Λʹ ^{- 1}, H^{*}= Λ^{-1} H Λʹ ^{- 1}, R^{*}= Λ^{-1} R Λʹ ^{- 1} and P^{*}= Λ^{-1} P Λʹ ^{- 1}
The matrices G^{*} , H^{*}, R^{*}and P^{*}have (co)variance components for sires of the calves, sires of the cows, herd-year-season of calving, residual and phenotypic effects, respectively, which first obtained under RMMs and then were transformed to their equivalents under SMMs. R^{*}is a matrixwith non-zero off-diagonal elements.
RESULTS AND DISCUSSION
General considerations
Pedigree structure of the considered first-parity Holstein cows is presented in Table 1. After editing of CD records and elimination of cows with atypical lactation curves, cows with easy and/or difficult calving constituted 74.50% and 25.50% of total cows, respectively. Approximately, 24% of the cows were excluded because of atypical shape of lactation curve. By fitting wood function pearson's correlation coefficient between actual and the corresponding predicted milk records was obtained as 0.96, implying that lactation curve in the first parity Iranian Holsteins has been described properly by wood function. Wood is one of the most popular functions used for describing the lactation curve in dairy cows (Cilek and Keskin, 2008; Atashi et al. 2012a; Boujenane and Hilal, 2012). Descriptive statistics for the studied lactation curve traits were shown in Table 2.
Table 1 Pedigree structure of the considered first-lactation Holstein cows in Iran
Table 2 Descriptive statistics for the studied lactation curve traits of first-parity Holsteins
Ap: initial milk yield at the beginning of lactation; Bp: inclining slope of the lactation curve; Cp: declining slope of the lactation curve; Pers: milk persistency; Y_{m}: peak milk yield and T_{m}: days to attain peak yield.
SD: standard deviation.
Milk production started with an average value of 16.14 kg at the beginning of the lactation, increased with an average slope of 0.29 and reached average peak milk yield of 39.44 kg at an average time of 97 days after calving. After peak milk yield attained it decreased with an average slope of 0.003. Least squares means ± standard error for the studied lactation curve traits in cows with unassisted calving (CD=1) against cows with any type of assisted calving (CD>1) are presented in Table 3. Cows with difficult calving had significantly lower Ap and higher Bp than cows with easy calving (P<0.05). The value of Ap in difficult calving cows was significantly lower than the corresponding values in easy calving cows. While Bp in difficult calving cows was 0.03 higher than easy calving ones. In other words, dystocia significantly decreased initial milk yield and increased inclining slope of lactation curve in first-parity Iranian Holsteins. Atashi et al. (2012a) reported that the inclining slope of lactation curve in the first-parity Iranian Holstein cows was significantly higher for cows with difficult calving than for those with eutocia which is in agreement with our finding in the present study. In the present study, there were no significant differences between cows experienced dystocia and those with eutocia in terms of the other studied lactation curve traits including Cp, Pers, Y_{m} and T_{m}. Dystocia has been identified as the major cause of stillbirth, resulting in reduced lactation (Atashi et al. 2012a) and reproductive performance (Mokhtari et al. 2016). In a previous study the causal relationships between gestation length, calving difficulty and still birth of primiparous dairy cows also has been studied (Lopez de Maturana et al. 2009).
Structural coefficients for the causal effects of CD on the lactation curve traits
Features of posterior means and posterior standard deviations (PSD) of the structural coefficients of CD on the studied lactation curve traits are presented in Table 4. CD is a categorical trait and analyzed applying threshold model. Therefore, the causal effects of CD on the other studied lactation-related traits were expressed on the liability scale. The estimated structural coefficients were significant only for Ap and Bp; highest posterior density (HPD) intervals did not include zero. It means that initial milk yield of first-parity Iranian Holstein cows and inclining slope of their lactation curve causally affected by calving difficulty. Posterior mean of recursive effect from CD on Ap was negative value of -0.351; implying that each 1-unit increase of liability for CD would decrease Ap by 0.351. Posterior mean of the recursive effects from CD on Bp was 0.005; suggesting that each 1-unit increase of liability for CD would increase Bp 0.005. The estimated causal effects of CD on Cp, Pers, T_{m} and Y_{m} were not significant; highest posterior density (HPD) intervals included zero. It implies that these traits were not affected by CD.
Table 3 Least squares means ± SE for the studied lactation curve traits across easy and difficult calving cows
Ap: initial milk yield at the beginning of lactation; Bp: inclining slope of the lactation curve; Cp: declining slope of the lactation curve; Pers: milk persistency; Y_{m}: peak milk yield and T_{m}: days to attain peak yield.
Table 4 Posterior means and posterior standard deviations (PSD) for the structural coefficients of calving difficulty on lactation curve traits
CD: calving difficulty; Ap: initial milk yield at the beginning of lactation; Bp: inclining slope of the lactation curve; Cp: declining slope of the lactation curve; Pers: milk persistency; Y_{m}: peak milk yield and T_{m}: days to attain peak yield.
** (P<0.01) and * (P<0.05).
NS: non significant.
As structural coefficients measure recursiveness at the phenotypic level (Gianola and Sorensen, 2004) and are corrected for genetic effects it seems that they describe phenotypic relationships more accurately. In a previous study, Atashi et al. (2012a) quantified the effect of dystocia on the lactation performance traits of Iranian Holstein cows without considering a cause-and-effect relationship from CD on these traits. They found significantly (P<0.05) higher least square mean values for peak days in milk and milk persistency in difficult calving cows at the first lactation while peak milk yield was higher than in the easy calving cows.
Model comparisons based on DIC
The DIC values for CD-Ap and CD-Bp under SMMs were higher than the corresponding values obtained under RMMs (Table 5), implying the plausibility of considering a cause-and-effect from CD on Ap and Bp. The lower DIC values for the RMMs are partly because of a lower penalty for model complexity in the DIC for RMMs (Bouwman et al. 2014). The RMMs present sources of covariance from the causal associations, in these models the uncorrelated residuals for the RMMs were assumed (Bouwman et al. 2014). The DIC values for CD-Cp, CD-Pers, CD-Y_{m} and CD-T_{m} under RMMs were higher than those obtained under SMMs which denote that a standard bivariate mixed model may better explain relationship between CD and these traits than a recursive bivariate mixed model considering CD as a parent trait.
Genetic parameter estimates under RMMs and SMMs
Direct heritability estimates
The posterior means and PSD for direct heritability of the studied lactation curve traits under SMMs and RMMs are presented in Table 6. It should be noted that parameters presented as under "RMMs" are pertaining to the standard equivalent. Similar estimates were obtained for direct heritability estimates of the studied lactation curve traits under both models. Heritability estimates for the studied traits were generally low; under both SMMs and RMMs, the lowest heritability was obtained for Bp (0.039 under SMM and 0.038 under RMM) and the highest one for Y_{m} (0.189 under SMM and 0.185 under RMM). Generally, heritability estimates for the studied lactation curve traits of first-parity Iranian Holsteins denoted low additive genetic variations and high influence of non-genetic effects for these traits. There were no significant differences between estimated values of direct heritability for the studied lactation curve traits under RMMs and SMMs. The low heritability estimates for the studied lactation curve traits imply that these traits are mainly dependent on non-genetic variations (Boujenane and Hilal, 2012). The estimated values for direct heritability of the considered lactation curve traits in the present study were in general agreement with those reported by Boujenane and Hilal (2012) in Moroccan Holstein-Friesian cows. Rekaya et al. (2000) found heritability estimates of 0.14, 0.26, and 0.05 for persistency, Y_{m} and T_{m}, respectively in Holstein-Friesian cows via a Bayesian approach. Muir et al. (2004) reported heritability estimates of 0.18 and 0.09 for persistency and day in milk at peak milk yield in first-lactation Canadian Holsteins.
Genetic and phenotypic correlations
The posterior means and PSD for the direct genetic and phenotypic correlations of CD with the lactation curve traits under SMMs and RMMs are shown in Table 7.
Table 5 Deviance information criterion (DIC) values for the bivariate analyses under standard mixed models (SMMs) and recursive mixed models (RMMs)
CD: calving difficulty; Ap: initial milk yield at the beginning of lactation; Bp: inclining slope of the lactation curve; Cp: declining slope of the lactation curve; Pers: milk persistency; Y_{m}: peak milk yield and T_{m}: days to attain peak yield.
Table 6 Posterior means and posterior standard deviations (PSD) for direct (h^{2}_{d}) heritability estimates of lactation curve traits under standard mixed models (SMMs) and recursive mixed models (RMMs)
Ap: initial milk yield at the beginning of lactation; Bp: inclining slope of the lactation curve; Cp: declining slope of the lactation curve; Pers: milk persistency; Y_{m}: peak milk yield and T_{m}: days to attain peak yield.
** (P<0.01) and * (P<0.05).
Table 7 Posterior means and posterior standard deviations (PSD) for direct genetic (r_{d12}) and phenotypic (r_{p12}) correlations between CD and lactation curve traits under standard mixed models (SMMs) and recursive mixed models (RMMs)
CD: calving difficulty; Ap: initial milk yield at the beginning of lactation; Bp: inclining slope of the lactation curve; Cp: declining slope of the lactation curve; Pers: milk persistency; Y_{m}: peak milk yield and T_{m}: days to attain peak yield.
** (P<0.01) and * (P<0.05).
NS: non significant.
Table 8 Posterior means and posterior standard deviations (PSD) for herd-year-season of calving (r_{h12}) and residual (r_{e12}) correlations between CD and lactation curve traits under standard mixed models (SMMs) and recursive mixed models (RMMs)
CD: calving difficulty; Ap: initial milk yield at the beginning of lactation; Bp: inclining slope of the lactation curve; Cp: declining slope of the lactation curve; Pers: milk persistency; Y_{m}: peak milk yield and T_{m}: days to attain peak yield.
** (P<0.01) and * (P<0.05).
NS: non significant.
In general, most of the direct genetic and phenotypic correlations between CD and the studied lactation curve traits were statistically non-significant (95% HPD intervals included zero). All the direct genetic correlations between CD and the studied lactation curve traits were statistically non-significant (95% HPD intervals included zero) except for CD-T_{m}; which were 0.346 and 0.030 under SMMs and RMMs, respectively, suggesting that a more difficult birth at the first calving is moderately associated with a delay in time to peak in terms of direct genetic effects. Direct genetic correlation between CD and T_{m} varied significantly under SMMs and RMMs, since the posterior mean in one model did not fall within 95% HPD interval of the other model. It implies that selecting the appropriate model for these traits is of great importance. Muir et al. (2004) reported an estimate of 0.15 for genetic correlation between calving difficulty and peak days in milk in first-lactation of Canadian Holsteins which is lower than the estimated value under SMMs in the present study. Furthermore, Muir et al. (2004) reported a moderate estimate of 0.43 for genetic correlation between CD and persistency in first-lactation Canadian Holsteins which is not in agreement with our estimate in the present study. All the phenotypic correlations between CD and the lactation curve traits were statistically not significant (95% HPD intervals included zero) under SMMs (Table 7). Muir et al. (2004) obtained very low estimate of 0.01 for both phenotypic correlations of CD with persistency and times to peak milk yield in first-parity Canadian Holsteins. Under RMMs, all the phenotypic correlations between CD and the studied lactation curve traits were statistically not significant expect for CD-Ap (-0.054) and CD-Bp (0.052); 99% HPD intervals were not included zero. In other words, the phenotypic correlations vary significantly across both models for Ap and Bp; different magnitude of the phenotypic correlation estimates under SMMs and RMMs implied that the estimates may be interpreted differently according to the considered models.
Contemporary group and residual correlations
In Table 8, the posterior means and PSD for the contemporary group (herd-year-season of calving) and residual correlations between CD and the lactation curve traits under SMMs and RMMs are shown. All the estimated correlations between the contemporary group effects of CD with the lactation curve under both SMMs and RMMs were not statistically significant; 95% HPD intervals included zero. The estimated residual correlations for CD-Ap were significant, low and negative under both SMMs and RMMs (95% HPD intervals was not included zero). The residual correlations between CD and the other lactation curve traits were not significant under SMMs; 95% HPD intervals included zero. In other words, there are no correlations of CD with Bp, Cp, Pers, T_{m} and Y_{m }in terms of residual effects. The residual correlations between CD-Cp, CD-Pers, CD-T_{m} and CD-Y_{m} were not estimated under RMMs (Table 8) because there was found no causal relationships between CD and these traits in the present study. The genetic correlations between maternal genetic effects of CD and direct genetic effects of the lactation curve traits were non-significant (the estimates are not shown); 95% HPD intervals were included zero.
CONCLUSION
Causal effects of CD on the lactation curve traits in first-parity Iranian Holsteins were modeled applying recursive mixed models as a trait causally influences the subsequent lactation curve traits. After a difficult calving, the initial milk yield and inclining slope of lactation curve in first-parity Iranian Holstein cows were differently affected. Cows with calving difficulty have lower initial milk production but they recover by a quickly inclining slope of lactation curve (increasing production), and reaches the same peak milk production at the same time as cows without CD. Furthermore the milk persistency was identical in difficult and easy calving cows. In other words, a recovery and compensatory of cows with calving difficulty were observed. The existence of significant causal effects growing from CD on these traits indicated that inclusion of this trait in the genetic evaluation of Ap and Bp in first-parity Iranian Holsteins via RMMs is of importance. The application of the recursive model methodology gave clearance to distinguish between the effect of CD on Ap and Bp and the genetic relationship between CD and these two traits.
ACKNOWLEDGEMENT
The contribution of Animal Breeding Center of Iran, especially Mr. M.B. Sayyad Nejad, for providing the required data set is greatly acknowledged. The work was supported by the University of Tehran.