Genetic parameters by Bayesian inference for dual purpose Jaffarabadi buffaloes

Knowledge of genetic parameters is essential for improved reproductive management and increased yield. Quantitative analysis of genetic parameters is lacking for many breeds of buffaloes. This article provides the first estimate of genetic parameters for dual purpose (meat and milk) Brazilian Jaffarabadi buffaloes, using Bayesian inference. Data on milk yield (MY), lactation length (LL), weight at 205 days (W205) and 365 (W365) days of age, and average daily gain (ADG) from 205 to 365 days of age were collected in two herds. Bivariate analyses (using the program MTGSAM) were performed with the Gibbs sampler to obtain estimates of variance and covariance. Average lactation milk yield and lactation length were 1 620.2±450.9 kg and 257.6±46.8 days, respectively, and the mean values for weight traits (kg) were 181.6±63.3 (W205), 298.04±116.1 (W365), and 0.73±0.35 (ADG). Heritability estimates (modes) were 0.16 for MY, 0.10 for LL, 0.43 for W205, 0.48 for W365 and 0.32 for ADG. There was a high genetic correlation (0.96) between milk yield and lactation length and very high genetic correlations (0.99) between the three growth traits. Our data suggest that both milk production and growth traits have clear potential for yield improvement through direct selection in this dual purpose breed. The selection for weight at an early age would be successful and selection for MY can be performed in the first lactation.


Introduction
The water buffalo (Bubalus bubalis) was originally bred in Asia and has become a widespread dairy animal throughout the world (Jain et al. 2007, Borghese 2010. The buffalo has a variety of agricultural uses, as a working animal and for the production of meat and milk. It is also remarkably hardy, being well-adapted and productive in low-quality pastures (Pereira et al. 2007). Dual purpose (meat and milk) cattle systems are traditional in developing countries, and are characterised by low inputs and the use of crossbred cattle (Ortega et al. 2007). Similar data on agricultural use of buffaloes in South America are lacking, although it seems likely that the majority of herds are dual purpose and some are also used as working animals.
Buffaloes of Murrah, Mediterranean, Jaffarabadi, and Carabao breeds were introduced to Brazil at the beginning of the 20th century (Cassiano et al. 2004) and recent estimates suggest that there are about 2.8 million buffalos distributed across all Brazilian states (Ramos et al. 2006) out of a global population of more than 174 million individuals (FAO 2006). Perhaps surprisingly, Brazilian buffaloes have had an even higher population growth than cattle during recent years ).
There have been several studies on the breeding systems of Murrah and Mediterranean breeds in Brazil (Ramos et al. 2006, Malhado et al. 2008, Aspilcueta-Borquis et al. 2010. However, equivalent studies are lacking for the Jaffarabadi breed in this and other countries. The objective of this study was to estimate genetic parameters for dual purpose Jaffarabadi buffaloes using Bayesian inference. REML and Bayesian methods have been applied extensively in animal breeding to estimate covariance components and genetic parameters. The Bayes methodology provides a solution for the finite sample size problem, because an exact a posteriori distribution exists for each large or small data set from which inferences can then be drawn. When a large data set is analysed, a priori information tends to be subjugated by the likelihood function in the establishment of the a posteriori distribution. In this case, the parameter estimates are close to those obtained by methods based on likelihood functions. However, this may not be the case when the sample size is limited because the maximum likelihood procedure only produces well-defined properties when the sample size is sufficiently large (Gianola & Fernando 1986). Carneiro Júnior et al. (2007) reported that the Bayesian method is well suited for analysing small populations when extensive historical information is available.

Material and methods
Data on milk yield (n=788), lactation length (n=788), body weight at 205 (n=513) and 365 (n=513) days of age and body weight gain of 205 to 365 days of age (n=513) were collected from buffaloes born between 1990 and 2002 that belonged to two dual purpose Jaffarabadi herds. Variance components and heritability coefficients were calculated for milk yield (MY), lactation length (LL), weight at 205 (W205) and 365 (W365) days of age and weight gain of 205 to 365 days of age (ADG). We also estimated the genetic and phenotypic correlations of MY with LL, and W205 with W365 and ADG.
Bivariate analyses were performed with the Gibbs sampler using the program MTGSAM (Multiple Trait Gibbs Sampling for Animal Models -as described by Van Tassel & Van Vleck [1995]) to obtain estimates of variance and covariance.
The model a for MY and LL adopted, represented in matrix notation, was where y is a vector of observed traits, X is the incidence matrix of fixed effects, β is a vector of fixed effects (CG and CO), Z is the incidence matrix of additive genetic random effects, a is a vector of additive genetic random effects, W is the incidence matrix of permanent environmental random effect, p is a vector of permanent environmental random effects and e is a vector of random error effects. The model for W205, W365, and ADG included additive direct random effects and fixed effects for contemporary groups. Maternal effects were removed from the models, because in the preliminary analysis the estimative was null, possibly due to the small data set.
Four seasons of birth (January to March, April to June, July to September, and October to December) were used for the formation of the contemporary groups (CG). The CG contained animals of the same herd, season and year of birth.
To improve computation efficiency an inverted Wishart distribution was used as an a priori estimation of the covariance component, in which a non-informative prior was assumed for all parameters. Wishart density describes the distribution of sum of squares and sum of products for normally distributed random variables. The random effects were assumed to conform to a multivariate normal distribution. The residual effects were assumed to be normally distributed.
The MTGSAM program uses the Gauss-Seidel iterative method in a mixed model equation to obtain initial values for the fixed and random effects to be used in the Gibbs sampler. The initial numbers were arbitrarily obtained using a single chain with 100 000 iterations, and a burn-in period of a 10 000 samples was used with samples taken each 10 cycles.
The convergence diagnosis was analysed with the method described by Raftery & Lewis (1992) using an algorithm implemented in R software through the package BOA, Bayesian Output Analysis (Smith 2005). The diagnostics proposed by Raftery & Lewis (1992) include convergence to a stationary distribution which finds the required size of chain to accurately estimate the quantile of parameter functions.
The discarded burn-in periods were 1200, 1350, 1500 and 180 for the following respective analyses: MY-LL, W205-W365, W205-ADG and W365-ADG. The diagnosis of convergence by the Raftery & Lewis (1992) method generated sampling intervals of 7 (MY-LL), 8 (W205-W365), 13 (W205-ADG) and 16 (W365-ADG). It should be noted that successive samples are not independent, so it is necessary to discard several iterations between every two saved samples -as the process uses iterated Markov chains, the dependence decreases with the increase in the distance between iterations thereby increasing the independence between saved samples. The greater challenge in using samples with potential autocorrelation between values is that the estimated standard error of the mean will be biased, because the variance of mean is biased (Carlin & Louis 2000). If the specified number of iterations is adequate, the mean values of the a posteriori samples will be valid estimates of the a posteriori distribution of the parameters.
The descriptive statistics (mean, median, mode, and standard deviation) of the a posteriori distribution for each parameter were obtained from effective samples. The highest posterior density (HPD) region or confidence interval provides the interval that includes 95 % of samples and is a measure of reliability. The HPD can also be applied to non-symmetric distributions (Hyndman 1996).
A kernel estimator (available in the software SAS 8; SAS Institute Inc., Cary, NC, USA) was used to generate graphical descriptions of the a posteriori distributions. The general form of a kernel estimator is: The value of l determines the degree of averaging in the estimate of the density function and is called a smoothing parameter. A bandwidth l can be selected for each Kernel estimator by specifying c in the expression l=n-1/5Qc, where Q is the sample interquartile range of the Y variable and c is a specific constant. This formulation makes c independent of the units of Y.

Results and discussion
Average milk yield was 1 620.2±450.9 kg. Ramos et al. (2006) and Aspilcueta-Borquis et al. (2010) reported means of 1 813.5±697.40 kg and 1 650±687.0 kg, respectively for buffaloes of the Murrah breed. In a study conducted in Italy by Rosati & Van Vleck (2002) with dairy buffaloes, the average was even higher (2 286.8 kg). However, lower values have been reported in Mediterranean (1 042.5 kg), Jaffarabadi (1 062.8 kg), Murrah (1 418.4 kg) and crossbred (1 062.8 kg) individuals (Tonhati et al. 2000). The relatively low mean in the current study is probably associated with the lack of specialisation for milk. The mean milk production is similar to that observed by Malhado et al. (2009) Ramos et al. (2006) and Malhado et al. (2009) reported averages of 252 (crossbred) and 256 days (Murrah buffaloes).
Buffaloes appear to be unable to maintain a lactation period that is considered adequate for dairy cattle (approximately 305 days). However, a comprehensive study on different genetic groups of Holstein × Gir breeds in different environments in Brazil observed a wide variation in lactation length, with an average minimum of 148 days in ¼ Holstein-¾ Gir animals raised on an extensive feeding system, and maximum average of 327 days in 7/8 Holstein-1/8Gir animals using an intensive system (Facó et al. 2002). This study demonstrates the importance of systematically assessing the effect of different environments and feeding systems to fully understand the factors influencing lactation traits in each breed and environment.
The mean (0.19), median (0.18) and mode (0.15) of estimates of heritability coefficient for milk yield (MY) were relatively consistent (Table 1), with a slight positive skew (Figure 1). The mode is considered to be the most appropriate measurement for a posteriori distributions, and best reflects the higher frequency values (Wright et al. 2000). It is also important to recognise that other measures of central tendency such as the mean and median can also summarise a posteriori distributions, especially if the densities are approximately symmetrical -under which circumstances such measures are similar.
A study on Murrah buffaloes using Bayesian inference (BI) reported heritability similar and equal to 0.22 (Aspilcueta-Borquis et al. 2010). Similarly, Ramos et al. (2006), Tonhati et al. (2008) and Rodrigues et al. (2010) reported heritabilities of 0.21, 0.19 and 0.25 using REML on the same breed. In contrast, much higher variability (0.39) was reported by Araújo et al. (2008) using BI suggesting that there is considerable genetic variability in this breed with a high potential of genetic gain for milk yield through selection.
The mean, median and mode were similar (0.10) for lactation length. This result suggests that direct selection for this trait is unlikely to result in substantial gains. Studies have reported values comparable to the present ones: 0.10 for Murrah buffaloes (Rodrigues et al. 2010) and 0.08 for crossbred buffaloes (Malhado et al. 2009).
The repeatability (mode) of milk production traits were 0.54 (MY) and 0.12 (LL) (Figure 2). This indicates that selection of animals based on limited information about the first lactation will only result in improved milk yield. Other studies have reported repeatability coefficients for MY of 0.28 (Tonhati et al. 2000), 0.32 (Ramos et al. 2006), 0.51 (Malhado et al. 2009) and 0.33 (Rodrigues et al. 2010). Malhado et al. (2009) also estimated a repeatability of 0.12 for LL suggesting that there may be large oscillations in the duration of lactation.
The estimated correlations between milk yield and lactation length were very high (Table 2) with negatively skewed distribution and higher density in the class of 0.9 to 1 (Figure 3). Malhado et al. (2009) estimated a genetic correlation of 0.89.
A recent study of crossbred dairy cattle reported a genetic correlation of 0.92 between MY and LL (Vercesi Filho et al. 2007). These authors support not removing records on short lactations from the data set for genetic analyses in dairy cattle, as described by Madalena (1988). Removing these records will incorrectly reduce the estimated genetic variability and Table 1 Means, standard deviation, median, mode and higher a posteriori density interval of the genetic parameters for traits milk yield, lactation length, weight at 205 and 365 days of age and weight gain of 205 at 365 days of age  The recommendations of these authors for dairy cattle may be applied to buffaloes in that selection for milk yield could promote substantial changes in the lactation length. The heritability coefficients (mode) for body weight traits were high: 0.43 (W205) and 0.48 (W365) ( Table 1 and Figure 4). The heritability for ADG was 0.32. Malhado et al. (2008), utilizing REML, reported similar heritability values (0.42 for W205, 0.41 for W365) for Mediterranean buffaloes. Coefficients for the three traits were high and showed high additive genetic variability, indicating the potential for genetic gain by selection.
The genetic correlations between the three growth traits were close to 1, and the three a posteriori distributions are almost entirely within the values 0.95 to 1 ( Figure 5). The short HPD regions and small standard deviations reflect the high genetic correlation between traits.
The results of heritability and correlations between traits imply that selection at an early age would be successful for Jaffarabadi buffalos. In the three estimates, the genetic correlations were higher than the phenotypic correlations ( Table 2). The lesser magnitude of the phenotypic correlations may be attributable to a low environmental correlation affecting the traits under consideration. Table 2 Phenotypic and genetic correlations between milk yield with lactation length, weight at 205 with weight at 365, and weight gain of 205 at 365 with W205 and W365. In conclusion, Jaffarabadi buffaloes that originate from dairy herds can also be used for meat production, given satisfactory growth performance. The milk yield reported in this study, despite its consistently low values, can be accounted for due to the dual purpose of the herds. Milk yield and growth traits have clear potential for genetic improvement by direct selection. The selection for weight at an early age would be successful and selection for MY can be performed in the first lactation. The genetic correlation between MY-LL indicates that indirect selection using milk yield is a potentially beneficial strategy to improve LL.