Estimation of variance components of milk, fat, and protein yields of Tunisian Holstein dairy cattle using Bayesian and REML methods

A multi-trait repeatability animal model under restricted maximum likelihood (REML) and Bayesian methods was used to estimate genetic parameters of milk, fat, and protein yields in Tunisian Holstein cows. The estimates of heritability for milk, fat, and protein yields from the REML procedure were 0.21± 0.05, 0.159± 0.04, and 0.158± 0.04, respectively. The corresponding results from the Bayesian procedure were 0.273± 0.02, 0.198± 0.01, and 0.187± 0.01. Heritability estimates tended to be larger via the Bayesian than those obtained by the REML method. Genetic and permanent environmental variances estimated by REML were smaller than those obtained by the Bayesian analysis. Inversely, REML estimates of the residual variances were larger than Bayesian estimates. Genetic and permanent correlation estimates were on the other hand comparable by both REML and Bayesian methods with permanent environmental being larger than genetic correlations. Results from this study confirm previous reports on genetic parameters for milk traits in Tunisian Holsteins and suggest that a multi-trait approach can be an alternative for implementing a routine genetic evaluation of the Tunisian dairy cattle population.


Introduction
In the past decades, dairy cattle management in Tunisia has been oriented toward increased milk yield.Importation of cattle from developed countries (USA, the Netherlands, and Germany) has been the main element for genetic improvement of the Tunisian dairy cattle population.Rekik et al. (2003) and Hammami et al. (2007) reported that 60 % of all inseminations of dairy cows in Tunisia used Holstein semen.In recent years, more attention has been placed on milk quality traits in breeding programmes.Estimates of genetic parameters for milk yield in dairy cows are abundant in the literature (Ben Gara et al., 2006, 2012;Hammami et al., 2008aHammami et al., , 2009a)).However, investigation of genetic parameters of milk components (quality traits) and the relationships between milk yield and quality traits is lacking.Multivariate models are of fundamental importance in applied and the-oretical quantitative genetics (Gianola and Sorensen, 2004).Precise estimation of genetic parameters is typically difficult in multiple-trait models due to a large number of genetic parameters and to insufficient statistical information (Rekaya et al., 2003).In animal breeding, two major methods were particularly applied, restricted maximum likelihood (REML) and Bayesian methods.REML has emerged as the method of choice in animal breeding for variance component estimation (Neumaier and Groeneveld, 1997).In recent years, Bayesian methods were broadly used to solve many of the difficulties faced by conventional statistical methods and extend the applicability of statistics on animal breeding data.Furthermore, Markov chain Monte Carlo (MCMC) has an important impact in applied statistics, especially from a Bayesian perspective for the estimation of genetic parameters in the linear mixed effect model (Sorensen and Gianola, 2002;Hallander et al., 2010).The aim of this research was to use a multi-trait Published by Copernicus Publications on behalf of the Leibniz Institute for Farm Animal Biology.repeatability model to estimate genetic parameters for milk, fat, and protein yields in the Tunisian Holstein cattle with REML and Bayesian approaches.

Data
Data were provided by the Tunisian Genetic Improvement Center, Livestock and Pasture Office, Tunis.Original data from the official milk recording database included 242 096 completed lactation records of parities 1 to 6 on Holstein cows from 1997 through 2014.The number of test-day (TD) records for milk, fat, and protein yields were not equal.Fat and protein yields were missing in some TD records due to technical reasons.Only records that included milk, fat, and protein yields were retained.Lactations having the date of first test > 50 days from parturition and/or average interval between successive tests > 50 days were excluded.Lactations were extended to 305 days for cows milked to or beyond this point.Cows without pedigree information were discarded and cows aged < 20 or > 40 months at first calving were deleted.After editing for unreasonable production to avoid possible erroneous data for daily milk yield (< 3 and > 60), fat content (< 1.5 and > 5 %) and protein percentage (< 1 and > 5 %), a total of 113 492 records remained.These records were of 54 105 cows sired by 3517 Holstein bulls.The pedigree file included the animal's identification number, the sire, the dam, and the date of birth and the herd of origin for each animal.Descriptive statistics for the edited data set used in the analysis are shown in Table 1.

Model
For the analyses a multi-trait repeatability model was used.The model equation is where Y = y 1 , . .., y t is the vector of records for the t traits; b and a are vectors of location effects and of additive genetic values, respectively.p is a vector of permanent environmental effects.X, Z and W are the corresponding incidence matrices.And e is a residual vector.Fixed effects included herd × year × season of calving and age × parity.
Subclasses for age at calving were determined for each parity (3-month intervals from 24 to 35 months and ≤ 36 for the first lactation; ≤ 39 months and 3-month intervals from 40 to 48 months for the second lactation; ≤ 51 months and 3-month intervals from 52 to 60 months for the third lactation; and ≤ 63 months and 3-month intervals from 64 to 96 months for the fourth lactation and later).Four seasons were defined (autumn, winter, spring, and summer).Heritability for each trait was calculated as additive genetic vari-  ance σ 2 ai divided by the phenotypic or total variance σ 2 phi , where σ 2 phi = σ 2 ai + σ 2 pi + σ 2 ei , and For the REML procedure, convergence of the iterative process was declared when the relative differences of consecutive parameters were lower than 10 −10 ( Misztal et al., 2002).Standard deviations of correlations and heritabilities were obtained as developed by Dodenhoff et al. (1998).For the Bayesian analysis, variance components were estimated using the GIBBS2F90 programme (Misztal et al., 2002).Posterior means of variance components, heritability, and correlation estimates were obtained using 100 000 samples.After a burn-in of 20 000 samples, 1 out of 10 iterations was then kept for subsequent analysis, resulting in an effective sample of 8000 iterations.Convergence of Gibbs chains was monitored by visual inspections of plots of samples for selected parameters by POSTGIBBSF90 programme (Misztal et al., 2002).

Genetic parameters
Variance components for each trait estimated from data using the REML procedure are shown in Table 2. Summary statistics (mean, mode, median, standard deviation, and 95 % highest probability density interval for genetic, residual, and permanent variances) that were estimated by Bayesian analysis are presented in Table 3. Posterior means and standard  deviations for heritabilities, genetic, and permanent correlations between milk, fat, and protein yields are presented in Tables 4-5.The additive genetic variance estimates for milk, fat, and protein yields by Bayesian method were 313 070, 397.16, and 216.33 kg 2 , respectively.The corresponding estimates by REML method were 217 800, 286.2, and 185.4 kg 2 , respectively.On the other hand, the residual variance estimates for 305 days of milk, fat, and protein yields obtained by REML method were larger than those obtained by the Bayesian procedure.Furthermore, residual variances estimated by both methods were the largest components, particularly that obtained by REML procedure, which was 712 200 kg 2 for a 305-day milk yield.However, for all traits studied, the permanent variances were slightly smaller compared with the other variance components.Estimates of heritability were moderately low for all traits (milk, fat, and protein yields).Heritability was 0.21 (±0.05) for milk yield, 0.159 (±0.04) for fat yield, and 0.158 (±0.04) for protein yield.Overall, heritability estimates for all traits, milk, fat, and protein yields obtained by REML method were lower than those found by Bayesian method: 0.21 vs. 0.273 for 305day milk yields, 0.159 vs. 0.198 for 305-day fat yields, and 0.158 vs. 0.187 for 305-day protein yields.All genetic correlations obtained in the current study were high and positive.Genetic correlations were higher between milk and protein yields than those between milk and fat yields (0.94 vs. 0.89).The largest genetic correlation was observed between fat and protein yields (0.95).Permanent correlations were large (around 0.9) among all studied traits.Environmental correlation estimates among milk yield components were consistently higher than estimates of genetic correlations: 0.9 vs. 0.89 for milk with fat yields and 0.95 vs. 0.94 for milk with protein yields.
Values obtained in this study for heritabilities for 305day milk and fat yields are comparable to those found by Carabaño et al. (1989) for Spanish data using the REML procedure.Results obtained in this study basically agree with those obtained by Alijani et al. (2012) in terms of the comparison between both methods.Heritability estimates in the Iranian Holsteins population ranged from 0.13 to 0.26, from 0.1 to 0.17, and from 0.15 to 0.21 for milk, fat, and protein yields, respectively, in the first three lactations by REML procedure.Respective estimates obtained in the same study using Bayesian analysis ranged from 0.19 to 0.29, from 0.17 to 0.21, and from 0.2 to 0.25 for milk, fat, and protein yields, respectively.Values of variance components estimated with this Bayesian method were different from those obtained by Ben Gara et al. (2006) in the same population.In fact, genetic variances associated with 305-day milk yields were consistently larger than those found by Ben Gara et al. (2006).However, the magnitude of genetic variance obtained in this study was low compared to estimates in other dairy cattle populations (Meyer, 1984;Misztal et al., 1992;Dedkova and Wolf, 2001).These differences were probably caused by difficulties encountered by daughters of superior sires to express their genetic potential under harsh climatic conditions and limited feed resources (Hammami et al., 2008b).The permanent environmental variances were more than 10 % of total variances for milk, fat, and protein yields, in agreement with previous reports by Ben Gara et al. (2006).The high values of the environmental variance would be explained by poor management practices, feeding fluctuations during the year, and stressful climatic conditions, which may result in an additional variation that is permanently associated with each cow (Hammami et al., 2008a).Residual variances estimated by the both methods were the largest components, particularly by REML procedure.Residual variances obtained in this study had standard deviations larger than those found by Ben Gara et al. (2006) implying elevated heterogeneity in estimates.This result can be explained by the use of the multi-trait model associated with a large number of genetic parameters and hindered by lack of information (Rekaya et al., 2003).Data in this study included records of years 2010 and 2011 where the civil unrest had a sizeable impact on herd management, data recording (reduced herd sizes and herds being recorded), and all activities -in particular those related to animal breeding.Heritability estimates for milk, fat, and protein yields in this study were also comparable with 0.25, 0.17, and 0.21 obtained in the Tunisian Holsteins by Hammami et al. (2008b) using a TD random regression model.The genetic correlations for all traits, milk, fat, and protein yields were high, ranging from 0.89 to 0.95.Genetic correlations were higher between milk and protein yields than between milk and fat yields (0.94 vs. 0.89).Genetic correlation estimates in this study were in accordance in terms of relations among milk traits, but obtained values were larger than most estimates reported in other studies (Meyer et al., 1984;Carabaño et al., 1989;Dedkova and Wolf, 2001).Nevertheless, occasionally above 90 % genetic correlation estimates were reported in the literature (Rekaya et al., 1999;Jakobsen et al., 2002;Hammami et al., 2008b).Carabaño et al. (1989) in the same study found different genetic correlations between milk and fat yields in US and Spanish Holsteins (0.63 and 0.69 for two US samples vs. 0.94 for a Spanish sample).In addition, fat and protein yields are derived from fat and protein percentage, respectively.Therefore, they might be essentially fat percentage, subject to sampling errors.Hammami et al. (2008b) iterated that the processes of sampling and chemical analyses in stressful climatic conditions in Tunisia could seriously affect data recording quality, especially for fat and protein yields.High genetic correlation estimates reveal that milk, fat, and protein yields are effectively controlled by the same genes.No notable differences were found among genetic and permanent correlations in terms of magnitude and sign.Consideration of different levels of production and genotype by environment interaction in the literature have shown that heritability differs significantly among cow populations and production levels (Meuwissen et al., 1996;Rekaya et al., 2003;Strabel and Jamrozik, 2006;Hammami et al., 2009b).Heritability estimates were lower for countries with low milk compared with countries with high milk production levels (Ben Gara et al., 2006;Hammami et al., 2009a).Gengler et al. (2005) using a TD model reported that heritability estimates for TD milk yield were higher (0.25) for high-yield herds and lower (0.15) for low-yield herds.Furthermore, it is important to take into consideration the constraining climatic conditions in Tunisia.
In fact, maximum temperature exceeded 32  Bouraoui et al., 2002).In addition, the Tunisian dairy herd is characterized by small sizes and limited production levels (Rekik et al., 2003).All the environmental factors might be a possible explanation for the low heritability estimates observed in this study.Furthermore, Rekaya et al. (2003) argued that small herd sizes may cause problems when creating contemporary groups for genetic evaluation.Ugarte et al. (1992) and Carabaño et al. (2004) reported that contemporary groups were required to have at least five records to be included in the analysis.Variances were found to be heterogeneous across herds, and numerous methods to ac-count for that variation were developed (Meuwissen et al., 1996;Rekaya et al., 1999;Gengler et al., 2005).Gengler et al. (2005) reported that a method based on transformed regressors for random regression effects can be used to adjust for heterogeneity of TD yield co-variances.

Conclusions
Variance components of 305-day milk, fat, and protein yields were investigated by REML and Bayesian procedures using a multi-trait animal model.Moderate heritability estimates were found from the Tunisian data compared with those from other studies on Holstein populations probably because of the reduced additive genetic and the important environmental and residual variances observed in the Tunisian population.The large environmental variances can be explained by the Tunisian constraining conditions and the heat stress effects permanently associated with lifetime performance of the cow population.Arguably, this study suggests that results from both methods were reasonably similar to suggest both methods can be used.However, time of computing to store all observations with Bayesian procedures was greater than the corresponding REML.Genetic parameter estimates in this research might be included as integral elements in a routine genetic evaluation via a multi-trait repeatability model to consolidate ongoing breeding practices to improve the Tunisian Holstein population.Nevertheless, efforts should be made to improve data quality mainly for fat and protein daily records.

Table 1 .
Characteristics of data used in the analysis.
* Contemporary group was defined as (herd-year-season) and had to have at least five observations.

Table 2 .
Estimates * c 2 was calculated as the ratio of the permanent environmental variance to total variance.

Table 3 .
Summary of marginal distributions of the variance components for 305 days of milk, fat, and protein yields.HPD indicates high posterior density region.e3 are the residual variances associated with 305 days of milk yield, 305 days of fat yield, and 305 days of protein yield, respectively.

Table 4 .
Genetic (above diagonal) and permanent environmental (below diagonal) correlations (SD in brackets) and heritabilities (diagonal) for 305 days of milk, fat, and protein yields by Bayesian analysis.

Table 5 .
Genetic (above diagonal) and permanent environmental (below diagonal) correlations (SD in brackets) and heritabilities (diagonal) for 305 days of milk, fat, and protein yields by REML analysis.