Bayesian inference in a quantitative genetic study of growth traits in Nelore cattle ( Bos indicus )

The objective of this study was to estimate (co)variance components and genetic parameters for weights (W) and scrotal circumferences (SC) at 365 and 450 days of age, of Nelore (Bos indicus) cattle, using Bayesian inference in single and multiple-trait animal models. The fitted linear models included, besides the animal and residual random effects, the contemporary group (herd-year-season-sex-management group) and age-of-dam as “fixed effects”. The analyses were carried out using a Gibbs sampler implemented through the MTGSAM program. The heritabilities (in parentheses) obtained fitting single-trait models were W365 (0.49), W450 (0.52), SC365 (0.68) and SC450 (0.66). Estimates of means, modes and medians for genetic parameters obtained from marginal posterior distributions were similar for all traits. The W365 and SC365 can be considered as suitable traits to be included as selection criteria in genetic improvement programs, not only because of their relatively high heritabilities but also due to the fact that they can be measured when the animals are relatively young compared to the corresponding traits W450 and SC450. The Bayesian approach appears to be an appropriate alternative for estimating genetic parameters, and has the advantage over point estimation methods of allowing inferences on marginal posterior distributions.


Introduction
As markets become more and more competitive and globalized, efficiency becomes a basic condition for any sector.In Brazil, the beef cattle industry is one of the most important segments of the national agribusiness, with an annual turnover of more than three billion dollars.
With reference to beef cattle, it is important to emphasize that zebu cattle (Bos indicus) today represents about 80% of the Brazilian herd.Additionally, as observed by Eler et al. (1996), more than 70% of zebu cattle are of the Nelore breed because, at least in part, of their adaptability to tropical conditions and their high fertility.Nelore cattle exhibit great differences in performance according to region and cattle raising practices, and these differences are not only related to environmental effects but also due to the genetic make up of the population.These variations can be considered as indicative of the potential for genetic improvement through selection.Thus, Nelore breed efficiency improvement becomes a basic condition for increasing livestock productivity in Brazil.
The Nelore Breed Genetic Improvement Program (Programa de Melhoramento Genético da Raça Nelore, PMGRN), see www.ancp.org.br, has as one of its major goals the genetic improvement of the growth, carcass and fertility traits of Nelore cattle.This program was developed and conducted by the National Association of Breeders and Researchers, and hosted at the Department of Genetics of the College of Medicine of the University of São Paulo, Ribeirão Preto Campus, Ribeirão Preto, São Paulo state, Brazil.To achieve the genetic improvement objectives, the PMGRN relies on sound data collection and processing schemes, including the application of classical genetic evaluation procedures.Most genetic evaluations in animal populations are based on mixed (mainly linear) model methodology and best linear unbiased predictor (BLUP) theory (Henderson, 1984), which require good estimates of (co)variance components among the traits for which individual animals are evaluated.These components need to be correctly estimated, so that predicted expected progeny differences (EPD) can serve as an effective instrument for genetic improvement purposes.Additionally, accurate estimates of (co)variance components is important to animal breeding due to the fact that prediction error variances for predicted breeding values increase as estimated values deviate from true values (Henderson, 1975).
Due to their importance, not only to genetics, but also to other areas of knowledge, several statistical methods of (co)variance estimates have been proposed throughout the years.More recently, the introduction of computer intensive statistical methods and the increasing adoption of Bayesian approaches and techniques, generally referred to as Markov Chain Monte Carlo (MCMC) methods, has been favored.In this respect, Gibbs sampling (GS) is a Monte Carlo numerical integration method that allows inferences to be made about joint or marginal distributions, even if appropriate densities cannot be explicitly formed (Geman and Geman, 1984).The use of GS allows analysis of Bayesian posterior distributions that had been computationally intractable due to the numerical integration required to obtain those distributions.The method presents properties of interest for animal improvement and extensive statistical resources capable of supplying more and more accurate estimates.
The objective of the study described in this paper was to apply the Gibbs sampling method for (co)variance components and the estimation of genetic parameters to growth and scrotal circumference traits using field data from Nelore cattle.

Material and Method
We used data and pedigree information from 69,025 Nelore cattle (Bos indicus) from distinct herds participating in the PMGRN.In general, the cattle were raised on pasture in various regions of Brazil, where the predominant climate ranges from subtropical hot humid to tropical hot humid, with two distinct seasons, one dry and the other rainy.The traits studied were weights (W) and scrotal circumferences (SC) at 365 (W 365 and SC 365 ) and 450 (W 450 and SC 450 ) days of age.The contemporary groups consisted of the herdyear-season-sex-management group.Birth seasons were divided into four periods: January to March, April to June, July to September and October to December.Age-of-dam effects were grouped into six classes: = 3, 4, 5, 6, 7-10 and = 11 years of age at calving.
The original data were edited using the Statistical Analysis System software (SAS, 1996), with preliminary analyses for general file consistency.The variance compo-nents necessary for the estimation of genetic parameters were obtained by the Gibbs sampling method, with the application of the "Multiple Trait using Gibbs Sampler under Animal Model, or MTGSAM program ( Van Tassel and Van Vleck, 1996).
The analyses were conducted by fitting single and two-trait animal models.The general mixed linear model, in matrix notation, was: where y is a vector of observations, X and Z are known incidence matrices relating β and u, respectively, β is a fixed effects vector, u is a random vector of additive genetic effects and e is a random vector of residuals.
From these definitions, it is considered that: where: σ a 2 and σ e 2 are additive genetic and residual variances, respectively; A is the numerator relationship matrix among animals and I is the appropriate identity matrix.Distributional assumptions were u ~N(0, Aσ a 2 ) and e ~N(0, A σ e 2 ), respectively.In the two-trait model, Σ = ⊗ G A and R I R = ⊗ + , where G is the matrix of additive genetic effects for the two-trait; and R + is matrix of residuals for the traits measured in each animal.
Prior distributions for the (co)variances were estimated from an inverted Wishart (IW) distribution for animal genetic and residual effects of traits ( Van Tassel and Van Vleck, 1996).If a random variable, W, is an IW random variable, then the probability density function of W is: As indicated above, several parameters are needed to model the density functions of the random variables.The parameter v is the degree of freedom corresponding to the Wishart random variable and the degree of belief for the prior distribution.The matrix V describes the (co)variance structure of the variable W, and m is the dimension of V.The mean of W is V -1 /v*, where V -1 = v*V 0 , V 0 is the (co)variance matrix specified from prior information.
The (co)variance component means are calculated as the mean of the expected values in the whole chain, after the initial sample burn-in.The estimates of the posterior means for heritabilities and correlations are calculated as the means of the functions using the (co)variance components in each Gibbs sampling cycle, after the burn-in period.For implementation of GS we considered a 200,000 total cycle sampler chain, a 20,000 cycle burn-in and a 100 cycle thinning interval to obtain 1,800 samples of (co)variance components and genetic parameters.The GS was implemented using a long chain scheme (Jensen et al., 1994;Sewalem and Johansson, 2000) and after burn-in one sample was collected from every 100 samples to avoid possible correlation between consecutive samples.The Gibbsit program of Raftery and Lewis (1994) was used to evaluate burn-in for all (co)variance components and also to calculate the thinning interval to determine the frequency of retaining sampled values so that those samples are uncorrelated for each parameter.The initial values used reflected no previous knowledge of the parameters.

Results and Discussion
The results obtained for the means, variation coefficients and growth traits minimum and maximum in this study are described in Table 1.The fact that single-trait and two-trait analyses were used did not interfere with these preliminary results.This is because descriptive statistics is a tool for discovering and evaluating the quality and consistency of the information, so that (co)variance component estimation analyses can be started with confidence in the data file.
To obtain the (co)variance component and genetic parameter estimates, the analyses were carried out using the scheme (Jensen et al., 1994;Sewalem and Johansson, 2000), with a view for evaluating the posterior distribution densities obtained from Gibbs chain samples.The results of variance component and heritability estimates for growth traits, obtained through GS in single-trait analyses are presented in Table 2.
The high heritabilities estimates obtained for all the traits considered in this study indicated the existence of considerable genetic potential for gains obtainable through selection (Eler et al., 1996).Indeed, the emphasis given to selection for growth in genetic improvement programs has most often provided substantial responses.The heritability estimates for weights and scrotal circumferences at 365 and 450 days were similar.The magnitude of the heritability estimates for all traits can be considered sufficient to promote substantial genetic gains in the Nelore breed.In this study, the components of variance and heritability estimates obtained for growth traits using single-trait analyses were consistent with analogous values reported in the literature (Gressler et al., 2000;Mercadante et al., 2000;Pereira et al., 2000;Everling et al., 2001).
Bayesian analyses allowed the estimation of the marginal posterior densities of the genetic parameters, and with these densities it was possible to visualize the errors in the estimates of these parameters, making them more accurate (Gianola et al., 1994).Table 2 presents standard deviations, Monte Carlo errors and, for the point estimates, the respective means, modes and medians of the variance and heritability components of growth traits in single-trait analyses.It can be observed that, except for SC 450 , the mean, median and mode estimates for the variance components were similar, since for heritability values the measures of central tendency had identical values as expected in posterior distributions classified as normal distributions.For the SC 450 value the mode of the heritability estimate was greater than the mean or median, although this difference was small.In relation to the Monte Carlo error (Table 2) the standard deviations were very small for all variance components and genetic parameters, which indicated that the size of the Gibbs Bayesian inference for growth traits in Nelore cattle 547  chain was sufficient to estimate reliable posterior means.The Monte Carlo error is the error in parameter estimation due to the number of samples used from the Gibbs chain ( Van Tassel and Van Vleck, 1996), and is inversely proportional to the length of the Gibbs chain.
As shown in Figure 1, when the Gibbs sampling is correctly implemented the marginal posterior distributions appear stable and tend toward normality (Blasco et al., 1998;Garcia Cortés et al., 1998;Sewalem and Johansson, 2000).These distributions are an important aspect of the Bayesian analysis, which provides an alternative of great flexibility in genetic analyses, due to the possibility of generating more accurate estimates arising from marginal posterior distributions not possible with the usual Restricted 548 Faria et al.
Figure 1 -Histograms of posterior densities estimates of additive genetic variance (σ a 2 ), residual variance (σ e 2 ), and heritability (h 2 ) for weight at 365 (upper) and 450 days of age (upper middle), scrotal circumference at 365 (lower middle) and 450 days of age (lower) from a single-trait analyses using Gibbs sampling.
Maximum Likelihood (REML) methodology.Comparing the two methodologies and using Nelore field data, Magnabosco et al. (2000) reported that the Gibbs sampling estimates were greater than those obtained by REML, but with similar results to those of other previous studies on Nelore cattle (Mercadante et al., 2000;Everling et al., 2001).Resende (2000) reported that Bayesian inference estimates were more precise than REML estimates, and have the advantage of producing standard deviations and confidence intervals for (co)variance components and genetic parameters as well as the marginal posterior densities.According to Blasco (2001), both Bayesian and Frequentist schools of inference are well established and it is not necessary to justify why one or the other school is preferred.
Multi-trait analyses can explore inter-trait correlations while improving genetic prediction (Henderson, 1984).In our study we carried out two-trait analyses among growth traits, the (co)variance components and genetic parameter estimates obtained being presented in Table 3.The estimates of heritability for the weights at 365 and 450 days of age in two-trait analyses were similar to those obtained with single-trait analyses (Tables 2 and 3).With respect to scrotal circumferences, the two-trait heritability values were smaller than the single-trait estimates.The mean, mode and median values of (co)variance component estimates were similar, and genetic parameters were practically identical (Table 3).This may be confirmed from the marginal posterior distributions shown in Figure 2, in which the measures of central tendency did not differ.The marginal posterior density distributions, classified as Gaussian, or normal, distributions were also unaltered.The standard deviation of the Monte Carlo error was very small for all variance components and genetic parameters, indicating that the length of the Gibbs chain was sufficient to estimate reliable posterior means (Table 3).
Figure 2 also shows that the 365 day and 450 day heritability estimates were similar, with a genetic correlation of high magnitude.These results indicate that selection for growth may be made at earlier ages, as with 365 day weight.However, Sarmento et al. (2003) have argued that the further away from weaning then the more efficient selection has to be, because the performance would have a greater direct effect with decreasing influence of maternal effects.It is worth noting that, due to the need to identify earlier maturing and faster growing animals, this factor Bayesian inference for growth traits in Nelore cattle 549 : additive genetic variance of trait 1; σ a 12 : additive genetic covariance between traits 1 and 2; σ a 2 2 : additive genetic variance of trait 2; σ e 1 2 : residual variance of trait 1; σ e 12 : residual covariance between traits 1 and 2; σ e 2 2 : residual variance of trait 2; h 1 2 :heritability of trait 1; h 2 2 : heritability of trait 2; r g 12 : genetic correlation between traits 1 and 2; SD: standard deviation; MCE SD : standard deviation of the Monte Carlo error.
would not be enough to exclude one year weight as a selection criterion.
The genetic correlations between growth and scrotal circumference traits obtained in this study were of medium to high magnitude.Correlated responses in weights when selecting for scrotal circumference were considered by Everling et al. (2001).Additionally, the genetic correlation obtained for 365 day weight and scrotal circumference was greater than that found at 450 days of age, a result that also favors selection at 365 days.Similar results were reported by Gressler et al. (2000).As presented in Table 3, the genetic correlation estimate between the 365 and 450 days scrotal circumferences was 0.87, indicating that selection on one trait would be reflected directly on the progress of the other.Because of this and the fact that both traits presented high heritabilities, selection for 365 day scrotal cir-  cumference should be adopted (Gressler et al., 2000;Mercadante et al., 2000;Pereira et al., 2000).
In conclusion, the (co)variance components and genetic parameters estimates obtained in this study for growth traits were, in general, similar to most estimates found in the literature for Nelore cattle.For the 365 and 450 day scrotal circumferences, the results showed higher values than those generally reported in the literature, usually employing Restricted Maximum Likelihood methodology.Bayesian inference showed to be an appropriate alternative for estimating (co)variance components and genetic parameters.It provided reasonable estimates and greater flexibility than the usual likelihood estimates, mainly because of the inferences obtained by employing the posterior marginal distributions.The weight and scrotal circumference at 365 days of age may be considered advantageous as criteria for genetic selection compared to the corresponding traits measured at latter ages, not only because of the similar heritability estimates observed for the traits but because this methodology allows identification of superior genotypes in a shorter time.

Table 1 -
Descriptive statistics for all variables included in the analyses.
365 = weight at 365 days of age; W 450 = weight at 450 days of age; SC 365 = scrotal circumference at 365 days of age; SC 450 = scrotal circumference at 450 days of age.

Table 2 -
Descriptive statistics for variance components and heritabilities estimates for weight at 365 (W 365 ) and 450 (W 450 ) and scrotal circumference at 365 (SC 365 ) and 450 (SC 450 ) days of age, obtained from single-trait analyses using Gibbs sampling.MCE SD : standard deviation of the Monte Carlo error.

Table 3 -
Descriptive statistics for variance components and heritabilities estimates for weight at 365 days of age (W 365 ) and 450 (W 450 ), scrotal circumference at 365 (SC 365 ) and 450 (SC 450 ) days of age based on two-trait analyses using Gibbs sampling.