Hierarchical Bayesian models for genotype × environment estimates in post-weaning gain of Hereford bovine via reaction norms

It was evaluated statistical models with different assumptions to define the one that best describes the presence of genotype × environment interaction on adjusted post-weaning weight gain (PWG345) of Hereford cattle, through the study of reactions norms to the environment, obtained by random regression using a Bayesian approach. Four reaction norms hierarchical models (RNHM) were used through the INTERGEN program. The RNHMK uses the solutions of contemporary groups previously estimated by the standard animal model (AM) and considers them as environmental level for predicting the reaction norms and the RNHMS, which jointly estimate these two sets of unknowns. For both models, two versions were considered, one with a homogeneous (hm) and another with a heterogeneous (ht) residual variance. Based on the deviance information criterion and Bayes factor, RNHMshm showed the best fit to the data, and by the deviance based on conditional predictive ordinate, the best fit was the RNHMKht, whereas, by all the three criteria used, the worst fit was obtained by using the standard animal model. Heritabilities estimated on RNHM were increasing in the environmental gradients for PWG345, at -60 kg, 0 and +60 kg. The genetic correlation estimated between the level and slope of reaction norms was high, from 0.97 to 0.99, characterizing a scale effect on genotype × environment interaction. The reaction norms hierarchical models are efficient to describe the changes in variance components due to the environment and to describe the presence of genotype × environment interaction on PWG345 trait of Hereford


Introduction
When an organism produces a phenotype that varies as a continuous function of the environment to which it is exposed, this response is called the reaction norms (Woltereck, 1909).From a statistical point of view, with the use of covariance functions (Kirkpatrick & Heckman, 1990;Gomulkiewicz & Kirkpatrick, 1992), the model of reaction norms can be generalized to an infinite number of environments (or traits), enabling the study of animal's reaction norms to gradual changes in the production environment, using a regression analysis of the genotype's performance on the average performance value observed in each environment (De Jong, 1995;Falconer & Mackay, 1996).
The advantage of this model is that the selection response can be predicted not only in the phenotypic expression in each environment, but also in terms of the trait environmental sensitivity (robustness or plasticity) to the changes in the environment (De Jong & Bijma, 2002), being feasible its implementation in genotype × environment interaction studies.Gianola & Fernando (1986) introduced the Bayesian methodology in animal breeding.The foundations of the Bayesian method consist in describing all errors which may exist about a parameter, considering as error measurement the probability of the parameter of taking certain values (Faria et al., 2007).
More recently, the introduction of Monte Carlo methods based in Markov chains (sequences), called MCMC (Markov Chain Monte Carlo), has substantially contributed to enable the implementation of the Bayesian paradigm (Sorensen, 2002).The MCMC methods represent a family of iterative processes used to approximate the generation of multi-trait distribution samples (in Monte Carlo processes with Markov chains properties).The Gibbs sampling is a procedure of numerical integration, used in the estimation of joint and marginal distributions of all the model parameters through the resampling of all fully conditional distributions of the Markov chains (Blasco, 2001).
Thus, in order to identify the occurrence of genotype × environment interaction on the post weaning weight gain of Hereford cattle, a study on environmental reactions norms was carried out, via random regression, using a Bayesian approach.

Material and Methods
The data used in this study derived from information routinely collected in Hereford herds which are participants in the Programa de Melhoramento de Bovinos de Carne (PROMEBO ® , 2008), conducted by the Associação Nacional de Produtores -" Herd Book Collares "-with animals born from 1972 to 2003 and containing 63,796 observations data.
The evaluated trait was post-weaning weight gain standardized to 345 days (PWG 345), which considered the interval from weaning to long-yearling.Initially, preparation, formatting and data description analyses were done using routines in SAS language (SAS, 2000) described in Cardoso (2008).
Contemporary groups were formed including breederherd, year and season of production, sex, management group, diet and weighing date of the animals.These groups were formed to group the animals that had a common environment and to set the environmental gradient by their deviations in relation to the average for PWG345.These averages were standardized and the contemporary groups were grouped into three classes: -1 standard deviation, zero and +1 standard deviation (60 kg).
Extreme records beyond ± 3 standard deviations of their contemporary groups means, of days in post-weaning test and of yearling age were eliminated.Moreover, contemporary groups with less than six animals, sires with less than two calves, animals with repeated numbers, codes of sex (castrated or missing) and missing cow age were also eliminated.
Finally, it was tested the connectedness among contemporary groups based on the total number of genetics ties (minimum 10), through the AMC program (Roso & Schenkel, 2006), which eliminated the disconnected record of 1,765 animals belonging to 148 contemporary groups, resulting in 88,221 animals in the pedigree file and 62,004 PWG345 records.Key details about the database are listed in Table 1.
Five models were studied for the data analysis, namely: the animal model (AM) that ignores the genotype × environment interaction and estimates the animal genetic value and the environmental effects, which later are used as covariates in the linear reaction norms models.The genetic value of animal i in environment j is described as below: in which y ij is the PWG345 record of animal i in j environment; β, a vector of fixed effects (linear and quadratic coefficients for age of cow and calf); x´i, the corresponding incidence vector; X j , environment random effect (contemporary group); a i , additive genetic effect of animal i; and e ij , the residual error.Additionally, two methodologies of analysis were studied for the reaction norms hierarchical models (RNHM): the model proposed by Kolmodin et al. (2002) which uses environmental solutions from AM as covariates in the RNHM, called RNHM k : in which φ = fixed regression coefficient; a i = additive genetic value of the intercept or random level of the animal reaction norm of animal i; b i = random regression coefficient or slope of reaction norm of animal i in environment represented by X ^j; X ^j = prediction of X j obtained in ( 1) and e ij = residual error.Moreover, the proposition of Su et al. (2006) called here RNHM S , in which the estimates of environmental effects are obtained together with the reaction norm, that is X j and b i are jointly estimated in model below: (3) Two different assumptions were adopted for the residual variance of the models: (a) for the AM, RNHM K and RNHM S homoscedastic (RNHM K hm and RNHM S hm): e i ~N(0, σ 2 e ); and (b) for models RNHM K and RNHM S heteroscedastic (RNHM K ht and RNHM S ht): e i ~N(0,σ 2 eij ).Using a Bayesian approach, Su et al. (2006) proposed the RNHM S in which the environmental effects and the reaction norms are estimating in a single analysis, without the need to use results of previous analysis and, therefore, solving the limitation of using estimates of an unknown covariate as if the estimates were true values without considering the uncertainty about these estimates.Although the model is similar to RNHM K , the estimation process is different and simultaneous for contemporary group and norm of reaction.The effects of the contemporary group are the unknown covariates for reaction norms, since their solutions are used as covariate in order to obtain the reaction norms slopes.For this one step analysis, the software is referred to the column of the data file with the contemporary groups identifications and no longer the column containing the contemporary group solutions of the animal model previously studied.One can describe this model for all n records in matrix form as follows: in which y = n order observations vector; β = vector for the fixed effects of order p; ξ = = the vector of environmental order effects n x ; a = = the vector of intercepts with orderq; b = = the vector reaction norms slopes also order q and e (n x1) = the residual vector of order n, while X, E, Z and Ξ = incidence matrices.
The additive genetic variance in the environment X, σ 2 A |X, was obtained by Thus, the heritability was estimated by the ratio of genetic variance to the phenotypic variance (genetic + environmental), as follows: where σ 2 e = residual variance in environment X given by σ 2 e |X = σ 2 e × η ^x on the heteroscedastic model and simply σ 2 e on the homoscedastic model, in which η = variance of heterogeneity parameter in the environmental gradient (X), following the structural model proposed by Cardoso et al. (2005).
The inference was based on Markov Chain Monte Carlo methods (MCMC) of 220,000 cycles, after 22,000 cycles of burn-in, and the samples were stored every ten cycles, for AM, RNHM S ht, RNHM S hm, RNHM K ht and RNHM K hm, using the INTERGEN program (Cardoso, 2008).
Descriptive analyses were developed by means of UNIVARIATE Procedure in the SAS software (2000).With this analysis, posterior means, modes, percentiles (2.5 and 97.5%) and standard deviations were obtained for all parameters from their posterior marginal densities.
In order to verify the best adjustment model it was used the following criteria: the Deviance Information Criterion (DIC) (Spiegelhalter et al., 2002); Based Conditional Predictive Ordinate deviance (DCPO), as described by Gelfand (1996), and deviance based on the Monte Carlo estimator proposed by Newton & Raftery (1994) for the

Results and Discussion
The RNHM models were superior to the AM based on all evaluated criteria, presenting smaller deviances.Based on the deviance obtained by Bayes Factors (DBF) and the Deviance Information Criterion (DIC), the one step RNHM was superior than the two steps versions.However, the Deviance based on Conditional Predictive Ordinate (DCPO) indicated a better fit for MNHR K (two steps) when compared on the assumption of homogeneity (Table 2).
In RNHM S , Deviance Bayes Factors and Deviance Information Criterion criteria indicated that the assumption of homogeneity of variances was better than that of heterogeneity.However, for Conditional Predictive Ordinate criteria the response was reversed.At RNHM K , all three criteria agreed that the best assumptions would be the homogeneity of variances.
The evaluation of gain post-weaning standardized to 345 days by Cardoso et al. (2007a) on Angus cattle, as well as by Corrêa et al. (2007) on Devon cattle, when assessing the same trait, obtained similar results when compared to animal model with RNHM using the DIC.When using the DBF, the DIC and the CPO, Corrêa et al. (2009) found the best fit at RNHM homoscedastic.According to these results, the animal model is the least indicated compared to RNHM whereas homoscedastic RNHM was more suitable.Corrêa et al. (2009) using the diagnosis of Geweke (1992) for studying post-weaning gain in Devon cattle found convergence for the standard animal model (Z: 0.89 and p-value: 0.3752), for the models RNHMhm (Z:1.06 and p-value: 0.2911) and for the model RNHMht (Z: 1.22 and p-value: 0.2218), demonstrating that the Bayesian models indicate convergence to the stationary posterior distribution.
The heritability estimate for PWG345 of 0.10 ± 0.01 in the animal model (Table 3) is somewhat below the results reported in the literature, for example, Bullock et al. (1993) studying growth traits in the same breed with North-Americans herds presented a heritability for post-weaning gains of 0.15.This value is similar to that found by Cardoso et al. (2004), who obtained a heritability estimate of 0.19 ± 0.02 in Brazilian Angus cattle.
The heritability estimates for the standard animal model of this study were very close to those estimated for RNHM in low environments, indicating that the proportion of genetic variance with respect to the environmental factors in the animal model may be underestimated.The heritability estimate value of 0.10 for PWG345 is considered low according to Cardellino & Rovira (1987), indicating that the selection criterion based on the phenotype of the animal will result in small genetic gains.The heritability estimates for different environmental levels obtained in this study (Table 3) are consistent with those of Cardoso et al. (2007a) with estimates of heritability for environments low, medium and high of 0.09 ± 0.01, 0.36 ± 0.01 and 0.54 ± 0.01, respectively, and Cardoso et al. (2007b) with values of 0.18 ± 0.01, 0.29 ± 0.02, and 0.45 ± 0.02, respectively.Such results showed that heritability estimates increased in response to environmental improvement.The heritability estimates for models that include the genotype × environment interaction have been increasing in the environmental gradient and have been higher in more favorable conditions, ranging from 0.07 to 0.64 in the RNHM s hm, from 0.11 to 0.59 in the RNHM k hm, and for the RNHM s ht, between 0.12 and 0.46 and, finally, for the RNHM k ht between 0.08 and 0.23 (Table 3).
Genetic parameters change with the environment gradient, indicating the allocation of greater proportion of phenotypic variation to genetic factors in relation to environmental factors by reaction norms hierarchical models and greater response to selection, especially in better environments.
Genetic correlations between the level and slope of reaction norm were similar and of high magnitude among and within models, indicating that animals of higher average breeding value were the ones which presented a best response to environmental improvement.
Phenotypic variances were similar among the reaction norms models in their respective environmental gradients (low, medium and high); for the animal model it was similar to estimates in the average environments.Confirming the assumption made by Cardoso et al. (2005), the genetic variance ( 2 a s ) in RNHMhm was inflated because of the unadjusted residual variance ( 2 e s ) for environmental gradient in this model.In RNHMht, ocurred this 2 e s adjustment, reducing the 2 a s and therefore the estimate of h 2 .The RNHMht, assuming the existence of residual heterogeneity, better partitioned the phenotypic variance between their genetic and environmental components.
De Mattos et al. (2000), when studying genotype × environment interaction in a multitrait animal model with Hereford cattle herds in three countries, Canada (CA), United States (US) and Uruguay (UY), obtained estimates of direct h 2 for weaning weight of 0.18 and 0.21 (US-CA); 0.21 and 0.22 (US-UY), and 0.17 and 0.19 (UY-CA).Moreover, genetic correlations of 0.86 (US-CA), of 0.90 (US-UY) and 0.88 (CA-UY), evidencing little relevance of genotype × environment interaction in that study.
Using the restricted maximum likelihood (REML) method for covariance components estimation, Alencar et al. (2005) found evidence of genotype × birth period (first and second semester) interaction on weaning (WW) and yearling weights (YW), on average daily gain from weaning to yearling age (ADG) and the performance based on a principal components index (PCI) involving these three traits in a beef cattle herd of Canchim breed.The heritability estimates in the first and second semester for WW, YW, ADG and PCI were 0.41 and 0.40; 0.36 and 0.38; 0.09 and 0.12 and 0.44 and 0.44, respectively.Additionally, the genetic correlation between the first and second semester for these traits were 0.87, 0.97, 0.91 and 0.88; these high values characterized genotype × environment interaction of low magnitude, as described by Falconer (1987).

Table 1 -
Observed means, standard deviations (SD) and ranges for weaning weight, long-yearling weight, post-weaning gain, postweaning average daily gain and post-weaning weight gain standardized 345 days of spectral density at zero frequency obtained by the SPECTRA Procedure of SAS(SAS, 2002), for the first n A and last n B cycles of MCMC chain of length m.Extreme absolute values of Z i , for a two-tailed test, indicate rejection of the convergence test.

Table 2 -
Geweke diagnostic test, deviance based on Bayes factor (DBF), deviance information criterion (DIC) and deviance based on conditional predictive ordinate (DCPO)Values 1 st , 2 nd , 3 rd , 4 th and 5 th indicate increasing order models of best fit.

Table 3 -
Posterior means standard deviations of variances components and heritability at different levels of environmental gradient and genetic correlation between level and slope of reaction norms obtained through different models