|
|
||||||||
,
,1
* Department of Genetics and Biotechnology, Faculty of Agricultural Sciences, University of arhus, PO Box 50, DK-8830 Tjele, Denmark
Department of Large Animal Sciences, Faculty of Life Sciences, University of Copenhagen, Ridebanevej 12, DK-1870 Frederiksberg C, Denmark
Department of Animal Science, Ferdowsi University of Mashhad, 91775-116 Mashhad, Iran
1 Corresponding author: mohammad.shariati{at}agrsci.dk
| ABSTRACT |
|---|
|
|
|---|
Key Words: reaction norm unknown covariate residual variance heterogeneity herd-year environment
| INTRODUCTION |
|---|
|
|
|---|
The reaction norm model is an interesting alternative to the multiple trait approach because it can accommodate a large number of environmental levels with few parameters. It assumes that the response variable is linearly related to the covariate representing the environmental variable. The slope associated with a given individual is a measure of its environmental sensitivity, and the amount of variation in slope in the population indicates the degree of GxE that exists for the trait under consideration. Linear reaction norms have usually been used to detect GxE because they are simple to interpret (Kolmodin et al., 2002; Calus and Veerkamp, 2003; Oseni et al., 2004).
In the reaction norm model the value of the covariate is typically unknown, and an approximation used in the literature consists of replacing the covariate by the average phenotypic performance of the individuals in a given environment [e.g., herd-year (Kolmodin et al., 2002; Fikse et al., 2003; Hayes et al., 2003)]. Other suggestions are to estimate herd(-year) effects in a standard additive genetic model and to use these as proxies for the unknown covariates in the reaction norm model (Calus et al., 2002; Oseni et al., 2004), or using the reaction norm model, to replace the covariates by estimates obtained in an iterative manner (Calus et al., 2004). These methods have a number of shortcomings, especially when the trait of interest is under selection (Calus et al., 2004; Su et al., 2006). Su et al. (2006) developed a Bayesian reaction norm model for inferring genetic parameters and environmental covariates simultaneously. The advantage of this approach over the approximate methods was shown in a simulation study (Su et al., 2006).
Several studies have reported the presence of GxE in milk production traits using approximate methods (Kolmodin et al., 2002; Calus and Veerkamp, 2003; Hayes et al., 2003). The aim of this study is to apply the methodology proposed by Su et al. (2006) to investigate the magnitude of GxE for dairy production traits (milk, protein, and fat yield) in early lactation using a random regression model where covariates (herd-years) are treated as unknown and inferred jointly with the remaining parameters. In this study the model was also extended to account for heterogeneity of within herd residual variance.
| MATERIALS AND METHODS |
|---|
|
|
|---|
|
Statistical Model
Sampling Model.
The linear reaction norm model used to analyze the data was
![]() |
where yijklmp is the first test-day record for milk, protein, or fat from the ith animal in the jth herd-year, kth DIM, lth age at calving, mth calving year, and pth residual variance group (herd); fj is the random effect of the jth herd-year (3,775 classes); dmk is the fixed effect of the kth DIM (41 classes); agel is the fixed effect of the lth age at calving in months (22 classes); cym is the fixed effect of the mth year of calving (14 classes); het is the breed heterozygosity; r is the regression on breed heterozygosity; prop is the proportion of Holstein Friesian genes; s is the regression on the proportion of Holstein genes; a0i is the additive genetic level of the ith animal; afi is the regression on the unknown herd-year effect specific to animal i; and eijklmp is the residual effect. Herd-year effects were treated as random effects because the model with fixed environmental effects is not identifiable. It has been shown that treating herd-year effects as random can increase the accuracy of selection (Visscher and Goddard, 1993) and the predictive ability of the model (Babot et al., 2003). Because herd-year effects were treated as random, calving year was included in the model as a fixed effect to account for the effect of environmental time trend (Babot et al., 2003).
When there is an association between sire effects and herd effects, treating herd-year effects as random can lead to biased estimates of breeding values (Henderson, 1973). Visscher and Goddard (1993) found a substantial bias under negative association of sires and herd environments for all contemporary group sizes. Inclusion of the relationship among animals using an animal model can decrease this possible bias (Ugarte et al., 1992). Ugarte et al. (1992), in a simulation study with preferential use of sires across herds, found that there was no bias in predicted breeding values in scenarios where sizes of contemporary groups were larger than 12. This was under positive association between sires and herds, which is expected to be the case for production traits.
The sampling distribution of the observations was
![]() | [1] |
where b is a vector containing effects of DIM, age at calving, calving year, regressions on breed heterozygosity and proportion of Holstein Friesian genes; f is the vector of herd-year effects; a0 and af are the vectors of genetic levels and slopes, respectively; X, E, and Z are known incidence matrices; and Zf is the unknown incidence matrix constructed with elements using f to associate the genetic slopes to the records. The (co)variance structure of the model is
![]() | [2] |
where
is the (co)variance matrix of genetic level and slope, A is the additive genetic relationship matrix of all animals in the pedigree, I is the identity matrix, and
2f is the variance of herd-year effects. Heterogeneity of residual variation (Short et al., 1990; Ibanez et al., 1999) was taken into account by postulating a residual variance peculiar to each of the 299 herds in the data set. Therefore D = diag (
2ep), P = 1, 2, ..., 299, is a diagonal matrix with herd specific residual variances on the diagonal. Heterogeneity at the level of herd-years could not be ascertained due to numerical instability of the algorithm, presumably due to the large number of residual variance parameters in relation to the number of observations per herd-year.
Prior Distributions.
The vector b was assumed to have an improper prior distribution. The prior distributions of genetic levels and slopes as well as herd effects were assumed to be normal with mean equal to zero and with a (co)variance matrix defined in [2]. Scaled inverted
2 distributions were considered for the variance component of herd effects (
2f) and for the elements of D, and finally, an inverse Wishart distribution was specified for G.
Joint Posterior Distribution.
The joint posterior distribution of all parameters was
![]() | [3] |
where
= {b, a0,af} is the vector of location parameters.
In this model the breeding value of individual i in environment j was given by
![]() |
The genetic variance in environment fj was defined as
![]() | [4] |
The heritability in herd-year fj was defined as
![]() | [5] |
where
2a(j) is defined in [4] and
2ep is the residual variance in herd p.
The genetic covariance between environments fi and fj is defined as
![]() |
and the genetic correlation between environments fi and fj is
![]() | [6] |
If the slope is the same for all individuals,
2af =
a0 af = 0, and raij = 1, indicating absence of GxE. Herein-after, for each trait estimated herd-year effects f have been expressed in units of their standard deviation, yielding standardized herd-year effects f* = f/
2f. As a result, genetic slope variances and covariances between level, slope, and heritabilities have been transformed accordingly, as shown by Su et al. (2006).
Model Comparison
The linear reaction norm model with unknown covariates with herd specific residual variances (RNUCH) was compared with a reaction norm model with homogeneous residual variance (RNUCS) and to an animal model with herd specific residual variances but without random regressions (ANH). The deviance information criterion (DIC) was used as the criterion of comparison. Let
represent parameters of a model. The DIC is computed as
![]() |
where
(
) is the posterior mean of the deviance that is –2logp(y|
). The effective number of parameters PD(
) is determined by PD(
) =
(
) – D(
), where
is the posterior mean of
. Smaller quantities of DIC indicate a better fit, after penalization for model complexity (Spiegelhalter et al., 2002).
Implementation of the Gibbs Sampler
In the Gibbs sampler the parameters were updated from their fully conditional posterior distributions that can be derived from the joint posterior distribution in [3] (Sorensen and Gianola, 2002). In the Bayesian setting of the linear reaction norm model with unknown covariates, all the fully conditional distributions are of standard form and are given in Su et al. (2006).
Post Gibbs analysis of MCMC output from previous analyses of the present data with slightly different models indicated that an adequate strategy in terms of Monte Carlo standard errors was to execute the sampler using a burn in of length 100,000 followed by a chain of length 200,000 for milk and fat yield. Rate of convergence and effective chain length was smaller for protein yield; therefore, 200,000 iterations of the Gibbs sampler were considered as burn-in followed by a chain of length 300,000. Monte Carlo standard errors and effective chain sizes of the MCMC output of RNUCH have been displayed in Table 2
.
|
| RESULTS AND DISCUSSION |
|---|
|
|
|---|
|
Residual Variances
Histograms of Monte Carlo estimates of posterior means of residual (within herd) variances for milk, protein, and fat yield are shown in Figure 1a
. The averages across herds were 16.60, 0.016, and 0.040 for milk, protein, and fat yield, respectively. The empirical correlation between posterior means of residual variances of milk and protein was 0.93. The corresponding correlation coefficient between posterior means of residual variances of milk and fat yield was 0.71, and between protein and fat yield was 0.76. This indicates that environmental variability within herd is more similar for milk yield and protein yield, but less similar for milk yield and fat yield and for protein yield and fat yield. Using RNUCS, the estimated residual variances for milk, protein, and fat yield were 16.63, 0.016, and 0.038, respectively. It shows that by postulating a single residual variance in the reaction norm model, the residual variance in the mean environment is estimated.
|
Herd-year Effects
Posterior means of variances of herd-year effects for milk, protein, and fat yield were 4.13, 0.0046, and 0.0107, respectively. Figure 1b
displays the distribution of estimated herd-year effects. Herd effects, defined as the average of the estimated herd-year effects over years, are shown in Figure 1c
. The correlation between posterior means of herd effects and of herd-specific residual variances was smaller than 0.10 for all traits studied. This suggests that heterogeneity of residual variances cannot be removed by a simple transformation of the data (Hayes et al., 2003).
In the present study, only data from herds with records from at least 30 first lactation cows per year were included. If small herds had lower management levels, the exclusion of small herds would result in somewhat of a decrease in the variation of environmental effects. This would not have an impact on the estimate of variance of slope because the slope is independent of the range of environmental effects if GxE has a linear reaction norm over all environments. However, the genetic correlation between extremely poor and extremely good environments could be slightly overestimated, and the extent of reranking could be somewhat underestimated because of a possible exclusion of the worst environments.
Genetic Variance Components and Correlations
Table 4
shows Monte Carlo estimates of posterior means of various genetic parameters, including the elements of the covariance matrix G (
2a0,
a0af and
2af), the additive genetic variance (defined in equation [4]) in environments 1 standard deviation above and below the mean environmental effect (
2a(+) and
2a(–), respectively), and the genetic correlation (defined in equation [6]) between these environments (ra(–,+)).The posterior distribution of the variance of the slope
2afwas shifted away from zero for all traits under study, supporting the presence of GxE. The ratio of interaction (slope) variance to the genetic level variance was largest for fat production and smallest for milk yield, indicating that the amount of GxE is highest in the former. Fat yield showed the lowest genetic correlation between environments, followed by protein and milk yield. The genetic correlation indicates the extent of reranking of genotypes across environments (Falconer, 1990). The magnitude of the genetic correlation in the present study indicates that there cannot be major reranking of genotypes in the different environments (Robertson, 1959). The rank correlations between breeding values at 5th and 95th quantiles of the estimated environmental effects were respectively equal to 0.91, 0.90, and 0.76 for milk, protein, and fat yield. The correlations for milk and protein yield are high, and little reranking across extreme environments is expected. In the case of fat yield, high genetic variability in early lactation along with a high magnitude of GxE variance translates into variable response of genotypes to different environments and lower rank correlation across extremely different environments.
|
|
2a0,
a0af,and
2af, respectively. The corresponding ratios of milk to fat yield were (x1,000) 0.33, 0.16, and 0.11. The constant ratios of milk yield to protein yield suggest that the two traits have a similar reaction norm.
The maximum genetic variability of fat yield occurs during the first days of lactation (Swalve, 1995; de Roos et al., 2004). In the current study, the genetic level variance for fat yield was 4 times larger than the variance of genetic level for protein yield (Table 4
). This variance determines the additive genetic variance in the mean environment (at f = 0). Swalve (1995) detected the same result with first test-day records of first lactation Friesian cows from northern Germany, where the environmental conditions are similar to Denmark. In his study, additive genetic variances were slightly lower than the values in the present study.
If the aim of a breeding program is to obtain individuals that perform well across the whole range of the environmental gradient, then selection could be based on the (posterior mean of the) breeding values in the mean environment (at f = 0). The Pearson correlation coefficients of breeding values in the mean environment and the breeding values from an additive genetic model without genetic effects for slope were 0.99, 0.99, and 0.97, for milk, protein, and fat yield, respectively. The analogous Spearmans rank correlation coefficients were 0.98, 0.98, and 0.96. The high correlations show that a simple additive genetic model is reasonable when the purpose is to rank animals for overall performance across a range of environments.
Heritabilities
Trajectories of heritabilities (defined in equation [5]) for milk, protein, and fat yield over environmental effects are displayed in Figure 2
. The mean heritability shows the trajectory of heritability in a herd with average residual variance that is approximately the mean of heritabilities in each environment. The plots in Figure 2
show that both heritability and variation of heritability between herd-years increased with increasing herd-year effects. The increase in variation of heritability between herd-years was largest for fat and protein yield. These results are in agreement with the study on protein yield by Calus et al. (2002) and on milk production traits by Hayes et al. (2003), who observed higher heritabilities in the "best" environments. On the other hand, Calus and Veerkamp (2003) did not find heterogeneity in heritability for milk production traits across environmental values defined as average protein production, although heterogeneous genetic variances were detected.
| CONCLUSIONS |
|---|
|
|
|---|
| ACKNOWLEDGEMENTS |
|---|
|
|
|---|
Received for publication January 23, 2007. Accepted for publication September 4, 2007.
| REFERENCES |
|---|
|
|
|---|
This article has been cited by other articles:
![]() |
G. Su, P. Madsen, and M. S. Lund Reaction norm model with unknown environmental covariate to analyze heterosis by environment interaction J Dairy Sci, May 1, 2009; 92(5): 2204 - 2213. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |