|
|
||||||||

* Animal Breeding and Genetics Group, Wageningen Institute of Animal Sciences, Wageningen University, PO Box 338, 6700 AH Wageningen, The Netherlands
CR-Delta, PO Box 454, 6800 AL Arnhem, The Netherlands
Corresponding author:
M. P. L. Calus; e-mail:
m.p.l.calus{at}id.dlo.nl. Current address: Institute for Animal Science and Health, ID-Lelystad, PO Box 65, 8200 AB Lelystad, The Netherlands.
| ABSTRACT |
|---|
|
|
|---|
Key Words: genotype x environment interaction protein yield reaction norm model dairy cattle
Abbreviation key: G x E = genotype by environment, HYS = herd-year-season, MSE = mean squares of errors
| INTRODUCTION |
|---|
|
|
|---|
In this study, three different sire models were applied to quantify G x E interaction. Given phenotype P, genotype G, environment E, and genotype x environment interaction G*E, the three models can be represented as:
![]() | 1. |
![]() | 2. |
![]() | 3. |
where the underlined terms represent the sires breeding values. These are, respectively, a model with a sire x herd-year-season (HYS) interaction, a model with HYS subclasses divided in different groups based on an environmental parameter and the reaction norm model. In this study, model 1, the interaction model, included the sire by HYS interaction as a variance component and so calculated EBV excluding G x E interaction. Model 2 calculated EBV per group of environments; the G x E interaction becomes part of the breeding value per group. Model 3 calculated EBV as a linear function, the reaction norm, of an environmental parameter. This model is only recently introduced for use in dairy cattle breeding value estimation (Strandberg et al., 2000). A reaction norm is the environmental sensitivity of a genotype in different environments (Falconer, 1990). The EBV of model 3 are built up of two parts: an environment-independent part, the level, and an environment-dependent part, the slope. In fact, the level is the random intercept of the reaction norm, and the slope is the random linear coefficient of the random regression on the environmental parameter. This model was applied by using a random regression model (Strandberg et al., 2000) rather than a repeatability model (Ravagnolo and Misztal, 2000).
The results of all models were compared to results of a standard model that did not account for G x E interaction.
The objective of this study was to compare the different models for their ability to account for G x E interaction.
| MATERIALS AND METHODS |
|---|
|
|
|---|
The pedigree information, from 1980 onwards, contained 11,656 sires in total. Sires (n = 2150) with at least five daughters with a lactation record in the data were identified. EBV were calculated for all those sires, based on 102,899 records for protein yield. For the sire x HYS interaction model, 65,711 sire x HYS subclasses were formed.
Models
To show the influence of the heterogeneity of variances, some of the models were also applied to the data after correction for heterogeneity of variances. From here, uncorrected and corrected will refer, respectively, to the situations without and with correction on the data for heterogeneity of variances. For the correction, a simple method was used; both within HYS subclass variances and HYS subclass means were standardized for all animals in the data. The within variances for all HYS subclasses were calculated as:
![]() |
where yij is the phenotypic value of animal j in HYS subclass i, yi. is the mean of HYS subclass i, and ni is the number of phenotypic records in HYS subclass i.
The weighted average of all within-HYS standard deviations was arbitrarily used as base standard deviation
b. Also, the weighted average of all HYS averages was arbitrarily used as base average yb. The preadjusted record was (Lofgren et al., 1985; Ibáñez et al., 1996):
![]() |
This adjustment does not account for fixed effects. Because only first-lactation animals are involved, the within-herd variances were comparable to each other.
The variance components of all models were estimated by REML (Gilmour et al., 1999).
The standard model was:
![]() | [Model 0] |
where yijkl is the record of daughter l of sire k, hys
i is a fixed effect of herd-year-season subclass i, agej is a fixed effect of age class j, Gk is a random additive genetic effect of sire k, and eijkl is a residual effect. This model was applied for reasons of comparison and to estimate environmental effects, corrected for genetic influences, for all HYS subclasses. These hys
i effects were used as environmental parameters in the reaction norm model.
The interaction model was:
![]() | [Model 1] |
where Iik is a random effect of G x E interaction of sirek and HYS subclass i.
In matrix algebra models 0 and 1 are:
![]() | [0] |
![]() | [1] |
where
is[b', g',s'] for model 0, and [b', g', s',t'] for model 1. And, for model 0:
|
|
Also, for model 1:
|
|
where
2s is the variance of additive sire effects,
2t is the variance of interaction effects, and
2e is the variance of residual effects.
The heritabilities were calculated as
[0] and
[1].
In this model, the interaction variance as part of the total phenotypic variance explains the amount of G x E interaction.
Model 2 is model 0 applied independently to four different subgroups of the data. The division of the data is based on the estimated HYS effect in model 0 without correction for heterogeneity of variances. Covariances between effects in different groups were not taken into account.
In model 2, calculating the correlation between the sets of EBV shows the G x E interaction. The correlation depends on the reliability of the EBV. Reliability is defined as:
|
|
where P (prediction error variance) is equal to the squared standard error of the EBV and
. And approximately, r2= ne/(ne + k) Wilmink et al. (1983), where ne is the effective number of daughters, and
|
|
The effective number of daughters can be calculated as ne = n*m/(n + m), where n is the number of daughters of a sire and m is the number of daughters of other sires in the same HYS subclasses (Van Vleck, 1987). In this study, sires with a minimum number of 75 daughters in the protein data were included in the calculations of the correlations between different sets of EBV. This made a total of 107 sires with 75,683 daughters with records.
Correlations between the EBV were transformed to genetic correlations (Blanchard et al., 1982),
![]() |
where,
g is the estimator of the genetic correlation, The third model is
![]() | [Model 3a] |
where yijkl is the record of daughter l of sire j, HYSi reflects the HYS subclass i, agej is a fixed effect of age class j, Lk is the random level of the additive genetic effect of sire k, hys
i is the estimated effect of HYS subclass i as estimated by model 0, Sk is the value of the random slope of the additive genetic effect of sire k and eijkl is a residual effect.
The same model is applied in a situation in which, instead of taking HYS subclasses into account as a fixed effect, a fixed regression is performed on estimated HYS effects hys
i. So, hys
i is a fixed covariate of HYS subgroup i and b is the coefficient of a fixed regression of y on hys
i:
![]() | [Model 3b] |
For both models, the estimated HYS effects were scaled between –1.0 and +1.0. This scaling makes internal calculations in the used ASREML-package easier (Gilmour et al., 1999; M. H. Pool, pers. comm.). This scaling was done with the following formula:
![]() |
where
=
max–
min and,
max is the estimated HYS-effect of HYS subclass i as estimated by model 0,
maxis the maximum estimated HYS-effect in the data, and
minis the minimum estimated HYS-effect in the data. For both models:
![]() |
with Gk is the additive genetic effect of sire k.
In matrix algebra model 3a and 3b are:
![]() | [3a] |
![]() | [3b] |
where
is[b', g',sL',sS'] for model 3a, and [ß',g',sL',sS'] for model 3b. As far as the variances are concerned, models 3a and 3b both have two variants. The first variant ignores the covariance between the variance of the level and the slope. These are referred to as respectively model 3a1 and 3b1. The second variant includes the covariance between level and slope. Those are referred to, respectively, as model 3a2 and 3b2. So,
![]() |
![]() |
|
|
![]() |
|
|
where
s2 is the variance of the genetic level of sires,
sS2 is the variance of the slope of the genetic effect of sires, and
sL,S is the covariance of the level and the slope of the genetic effect of sires. The obtained variance components of models 3a1 and 3b1 were used as prior values for respectively models 3a2 and 3b2.
Sire variances were calculated as
in models 3a1 and 3b1, and
in models 3a2 and 3b2. Heritabilities were calculated as
. (hys
)2 is the weighted average of the squared estimated HYS effects and hys
is the weighted average of the estimated HYS effects.
In this model, the G x E interaction is explained by the variance of the slope and the covariance between the level and the slope. All models are summarized in Table 1
.
|
|
|
with
i=v'i
i is the sum of the estimates of the relevant fixed and random effects according to equations 0| RESULTS |
|---|
|
|
|---|
|
In Table 3
, the variance components and heritabilities of model 2 for the uncorrected and the corrected data are shown. The minimum and average estimated HYS effects are given. The heterogeneity of sire variance is clearly visible. The residual variance of the corrected data does not show heterogeneity, due to the correction. The sire variance in general increased more than the residual variance with increasing level of estimated HYS effect. This resulted in higher heritabilities in groups with higher estimated HYS effect.
|
|
|
For 3b2, the covariance component was added to this: var(s) = var(sL) + (2*–0.0427)*cov(sL,sS) + (0.0517)*var(sS).
–0.0427 is the weighted average scaled estimated HYS-effect and 0.0517 is the weighted average of the squared scaled estimated HYS-effects.
The variance of the slope of model 3a1 for corrected protein was very small and not significant. Therefore, models 3b1 and 3b2 were only applied to the uncorrected data. The results of model 3b1 were used as prior values for model 3b2. The prior value of the correlation between level and slope was set to 0.80. This value was derived from Strandberg et al. (2000).
Comparing models 3b1 and 3b2 shows a difference between variances of level and slope. It should be noted that the multiplication factor of the covariance term is negative. Therefore, the sire variance of model 3b2 is smaller than the variance of the level (var(sL)).
A comparison between models 3a1 and 3b1 shows the influence of either including HYS as subclasses or as estimated HYS effect is needed. A decrease of both sire and residual variance is observed if HYS is an estimated HYS effect. However, the ratios between sire and residual variance are almost equal in these situations.
To show the distinctiveness of the reaction norm model, EBV of model 3b2 were calculated for all four groups from model 2. The Spearmans rank correlation coefficients between the groups are shown in Table 6
. The correlations are very high between the different groups. The average estimated HYS effect of groups 1 and 4 are –0.29 and 0.28, respectively, in the range of –1 to +1. This means that very few HYS subclasses are located in the extreme ends. Spearmans rank correlation between EBV in the two extreme ends, e.g., for scaled estimated HYS effects of, respectively, –1 and +1, was 0.42.
|
|
The correlation coefficients between the EBV of corrected protein yield were almost equal to those for uncorrected protein yield.
Mean Squares of Errors
Table 8
shows the MSE of the different models for both uncorrected and corrected protein yield. For uncorrected protein yield, model 1 clearly has the lowest MSE. After that, model 2 is second best. MSE of model 2 were calculated as an average of the groups. For uncorrected protein the MSE increases if the estimated HYS effect increases. Model 0 seems to be the worst model, as far as MSE is concerned. The reaction norm models are only slightly better. Correcting the data resulted in lower MSE for all models. Model 1 had the lowest MSE, but differences between models were small.
|
| DISCUSSION |
|---|
|
|
|---|
The interaction variance in model 1 was 2.5% of the phenotypic variance. The EBV were, however, close to those of the standard model at the cost of a higher calculation time. Ten Napel and Van der Werf (1992) found comparable results with an interaction variance of 3% of the total phenotypic variance for protein yield.
Model 1 did show the lowest MSE, but the benefit of this was not identified in correlation coefficients or estimates of variance components. The decrease of MSE implies that the EBV are more reliable, but still they are very close to those of model 0.
Judging from the results of model 2, it seems that there is more G x E interaction than model 1 detects. Probably, the sire x HYS interaction accounts only partly for G x E interaction. Therefore, it might be worthwhile to consider other possible definitions for G x E interaction then only sire by HYS. A possibility is dividing herds in different groups of production systems (e.g., in levels of intensity of production) and defining a sire by production system interaction.
The results of model 2 clearly show the heterogeneity of variances in the data. The model with the division of the data in groups showed correlations between 0.74 and 0.86 for EBV of different groups for uncorrected protein. Based on literature, correlations between 0.85 and unity were expected.
The records used in model 2 are divided over four equally sized groups and, therefore, the data are different from the other models. The fact that the correlations of EBV between groups are smaller than correlations with models 0, 1, and 3 might have originated because this model only uses a subset of the data that the other models use. The correlations might also be influenced by the method used to calculate them. Application of a multiple-trait model that uses all data together was tried but was not feasible. The results of such a model might be more comparable to those of the other models.
The reaction norm model has one major difference with model 0: the slope. The reaction norm model divides the standard breeding value in a level and a slope. It was expected that the slope would cause reranking of sires over environments. However, the level had such a great impact, that the slope had very little influence on the total breeding value. In fact, the reaction norm model was not able to detect any significant variance of the slope for corrected protein. Therefore, the EBV of corrected protein only contained the level. In both cases, the MSE and the correlations showed almost no differences with model 0.
Scaling the estimated HYS effects between –1 and +1 makes that the average scaled estimated HYS effect is close to 0. The level is then approximately equal to the breeding value for average HYS effect and the EBV of model 0.
In Figure 1
it is shown that for protein yield for eight randomly chosen sires, only two change rank within the range of scaled estimated HYS effects from –0.25 to 0.25. In this range, 75% of the records are located. However, outside this range reranking does occur, but very few animals and herds are located in both extreme ends. The reranking was only identified by the correlation of EBV of the two extreme situations. This might indicate that the model is able to calculate environment specific EBV, but that the environmental parameter was not distinctive enough for the majority of the herds.
|
A few remarks should be made according to the models used. All models applied were sire models. In other research (Schaeffer et al., 2001), it has been shown that application of alternative animal models gives higher reranking for cows than for bulls. Therefore, even if EBV of sires hardly change, the application of an alternative animal model might be considered.
The results of the model with the division in groups might indicate that the relation between environment and EBV is not linear, as was assumed in the applied reaction norm model. In that case, it is logical that calculating EBV over smaller intervals of an environmental parameter is better. EBV are not forced to be on a certain line as in the reaction norm model but are estimated based on the smaller amount of data within the interval. This indicates that a model with a nonlinear reaction norm should be considered. Another difference between models 2 and 3 is the way the residual is calculated. Model 2 calculates a specific residual for each group of environments. Model 3 calculates only one residual for all environments. This might explain the fact that the slope is not specific enough. One way to solve this problem is to include a permanent environmental effect as a function of the environmental effect as done by Ravagnolo and Misztal (2000).
| CONCLUSIONS |
|---|
|
|
|---|
Further research is needed to explore the characteristics and abilities of the reaction norm model further.
Received for publication January 31, 2002. Accepted for publication March 10, 2002.
| REFERENCES |
|---|
|
|
|---|
This article has been cited by other articles:
![]() |
H. Hammami, B. Rekik, H. Soyeurt, C. Bastin, J. Stoll, and N. Gengler Genotype x Environment Interaction for Milk Yield in Holsteins Using Luxembourg and Tunisian Populations J Dairy Sci, September 1, 2008; 91(9): 3661 - 3671. [Abstract] [Full Text] [PDF] |
||||
![]() |
M. Lillehammer, M. E. Goddard, H. Nilsen, E. Sehested, H. G. Olsen, S. Lien, and T. H. E. Meuwissen Quantitative Trait Locus-by-Environment Interaction for Milk Yield Traits on Bos taurus Autosome 6 Genetics, July 1, 2008; 179(3): 1539 - 1546. [Abstract] [Full Text] [PDF] |
||||
![]() |
M. M. Shariati, G. Su, P. Madsen, and D. Sorensen Analysis of Milk Production Traits in Early Lactation Using a Reaction Norm Model with Unknown Covariates J Dairy Sci, December 1, 2007; 90(12): 5759 - 5766. [Abstract] [Full Text] [PDF] |
||||
![]() |
A. G. Fahey, M. M. Schutz, D. L. Lofgren, A. P. Schinckel, and T. S. Stewart Genotype by Environment Interaction for Production Traits While Accounting for Heteroscedasticity J Dairy Sci, August 1, 2007; 90(8): 3889 - 3899. [Abstract] [Full Text] [PDF] |
||||
![]() |
M. Lillehammer, M. Arnyasi, S. Lien, H. G. Olsen, E. Sehested, J. Odegard, and T. H. E. Meuwissen A Genome Scan for Quantitative Trait Locus by Environment Interactions for Production Traits J Dairy Sci, July 1, 2007; 90(7): 3482 - 3489. [Abstract] [Full Text] [PDF] |
||||
![]() |
J. R. Bryant, N. Lopez-Villalobos, J. E. Pryce, C. W. Holmes, D. L. Johnson, and D. J. Garrick Environmental Sensitivity in New Zealand Dairy Cattle J Dairy Sci, March 1, 2007; 90(3): 1538 - 1547. [Abstract] [Full Text] [PDF] |
||||
![]() |
W. J. Nauta, R. F. Veerkamp, E. W. Brascamp, and H. Bovenhuis Genotype by environment interaction for milk production traits between organic and conventional dairy cattle production in The Netherlands. J Dairy Sci, July 1, 2006; 89(7): 2729 - 2737. [Abstract] [Full Text] [PDF] |
||||
![]() |
G. Su, P. Madsen, M. S. Lund, D. Sorensen, I. R. Korsgaard, and J. Jensen Bayesian analysis of the linear reaction norm model with unknown covariates J Anim Sci, July 1, 2006; 84(7): 1651 - 1657. [Abstract] [Full Text] [PDF] |
||||
![]() |
J. J. Windig, M. P. L. Calus, B. Beerda, and R. F. Veerkamp Genetic Correlations Between Milk Production and Health and Fertility Depending on Herd Environment J Dairy Sci, May 1, 2006; 89(5): 1765 - 1775. [Abstract] [Full Text] [PDF] |
||||
![]() |
S. Konig, G. Dietl, I. Raeder, and H. H. Swalve Genetic Relationships for Dairy Performance Between Large-Scale and Small-Scale Farm Conditions J Dairy Sci, November 1, 2005; 88(11): 4087 - 4096. [Abstract] [Full Text] [PDF] |
||||
![]() |
J. J. Windig, M. P. L. Calus, and R. F. Veerkamp Influence of Herd Environment on Health and Fertility and Their Relationship with Milk Production J Dairy Sci, January 1, 2005; 88(1): 335 - 347. [Abstract] [Full Text] [PDF] |
||||
![]() |
B. J. Hayes, M. Carrick, P. Bowman, and M. E. Goddard Genotype x Environment Interaction for Milk Production of Daughters of Australian Dairy Sires from Test-Day Records J Dairy Sci, November 1, 2003; 86(11): 3736 - 3744. [Abstract] [Full Text] [PDF] |
||||
![]() |
M. P. L. Calus and R. F. Veerkamp Estimation of Environmental Sensitivity of Genetic Merit for Milk Production Traits Using a Random Regression Model J Dairy Sci, November 1, 2003; 86(11): 3756 - 3764. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |