- Research
- Open Access

# Estimation of genetic variance for macro- and micro-environmental sensitivity using double hierarchical generalized linear models

- Han A Mulder
^{1, 2}Email author, - Lars Rönnegård
^{3, 4}, - W Freddy Fikse
^{4}, - Roel F Veerkamp
^{1}and - Erling Strandberg
^{4}

**45**:23

https://doi.org/10.1186/1297-9686-45-23

© Mulder et al.; licensee BioMed Central Ltd. 2013

**Received:**1 August 2012**Accepted:**24 April 2013**Published:**4 July 2013

## Abstract

### Background

Genetic variation for environmental sensitivity indicates that animals are genetically different in their response to environmental factors. Environmental factors are either identifiable (e.g. temperature) and called macro-environmental or unknown and called micro-environmental. The objectives of this study were to develop a statistical method to estimate genetic parameters for macro- and micro-environmental sensitivities simultaneously, to investigate bias and precision of resulting estimates of genetic parameters and to develop and evaluate use of Akaike’s information criterion using h-likelihood to select the best fitting model.

### Methods

We assumed that genetic variation in macro- and micro-environmental sensitivities is expressed as genetic variance in the slope of a linear reaction norm and environmental variance, respectively. A reaction norm model to estimate genetic variance for macro-environmental sensitivity was combined with a structural model for residual variance to estimate genetic variance for micro-environmental sensitivity using a double hierarchical generalized linear model in ASReml. Akaike’s information criterion was constructed as model selection criterion using approximated h-likelihood. Populations of sires with large half-sib offspring groups were simulated to investigate bias and precision of estimated genetic parameters.

### Results

Designs with 100 sires, each with at least 100 offspring, are required to have standard deviations of estimated variances lower than 50% of the true value. When the number of offspring increased, standard deviations of estimates across replicates decreased substantially, especially for genetic variances of macro- and micro-environmental sensitivities. Standard deviations of estimated genetic correlations across replicates were quite large (between 0.1 and 0.4), especially when sires had few offspring. Practically, no bias was observed for estimates of any of the parameters. Using Akaike’s information criterion the true genetic model was selected as the best statistical model in at least 90% of 100 replicates when the number of offspring per sire was 100. Application of the model to lactation milk yield in dairy cattle showed that genetic variance for micro- and macro-environmental sensitivities existed.

### Conclusion

The algorithm and model selection criterion presented here can contribute to better understand genetic control of macro- and micro-environmental sensitivities. Designs or datasets should have at least 100 sires each with 100 offspring.

## Keywords

- Genetic Correlation
- Residual Variance
- Dairy Cattle
- Genetic Parameter
- Sire Model

## Background

The term “genotype by environment (G × E) interaction” refers to the fact that the best genotype in one environment may not be the best genotype in another environment[1] and that genotypes differ in their response to environmental factors, which means that genetic variance for environmental sensitivity or phenotypic plasticity exists[2]. Some environmental factors (e.g., temperature, soil, diet, etc.) are identifiable and can be categorised (e.g., temperate or tropical climate) or quantified (e.g., temperature) and thus are referred to as macro-environmental factors. Other environmental factors are unknown and referred to as micro-environmental factors[1]. Therefore, genetic variance in macro-environmental sensitivity is the genetic variance due to known environmental factors and can be expressed as the genetic variance in the slope of a reaction norm when environments can be quantified on a continuous scale. If environments are categorized, then phenotypes in different environments are considered as separate traits and the genetic covariances between environments are a measure of genetic variation in macro-environmental sensitivity. Genetic variance in micro-environmental sensitivity is the genetic variance due to unknown environmental factors and can be expressed as differences in environmental variance, sometimes called genetic heterogeneity of environmental variance[3].

Numerous studies in the last 70 years have studied G × E interactions or genetics of macro-environmental sensitivity in animal and plant breeding as well as in evolutionary biology[4–6]. Different modeling approaches have been used[4]. In recent years, reaction norm models have been applied in animal breeding to better understand the environmental factors that determine G × E interactions[7–9]. In evolutionary genetics, many experiments on Drosophila and other laboratory species have been carried out to understand the genetics of macro-environmental sensitivity or phenotypic plasticity[6, 10, 11]. Other examples of studies on wild life populations include analyses of great tit[12] and butterfly[13] populations that showed the existence of genetic variation in phenotypic plasticity to temperature change.

The genetics of micro-environmental sensitivity or environmental variance have been studied less extensively than the genetics of macro-environmental sensitivity. In evolutionary genetics, several studies have focused on canalization, i.e. selection for reduced variance[14–16]. Recently, there has been renewed interest on this topic due to the development of methods to estimate genetic variance in environmental variance, e.g. Bayesian methods[17, 18] and double hierarchical generalized linear models (DHGLM) in a REML setting[19]. Hill and Mulder[20] reviewed 14 studies on this subject and concluded that there is empirical evidence for the existence of genetic variance in environmental variance.

Although there is substantial evidence for genetic variance in macro- and micro-environmental sensitivities, very few studies have studied them together or studied their genetic relationship. Jinks and Pooni[21] postulated the concepts of macro- and micro-environmental sensitivity and studied both in *Nicotiana rustica* but did not investigate the relationships between these two types of environmental sensitivity. However, Perkins and Jinks[22] reported that both types of environmental sensitivity were weakly genetically correlated for most traits in *Nicotiana rustica*. In a study on *Populus*, Wu[23] showed that micro-environmental sensitivity was less heritable than macro-environmental sensitivity and that they were weakly genetically correlated. In another paper, also on *Populus*, Wu and O'Malley[24] detected different sets of genes for macro- and micro-environmental sensitivities. In a mapping study of QTL for macro- and micro-environmental sensitivities in barley, Kraakman et al.[25] identified only a few QTL for micro-environmental sensitivity, but did not detect any QTL for macro-environmental sensitivity. Yampolsky and Scheiner[26] observed a weak genetic correlation between developmental noise, similar to micro-environmental sensitivity, and macro-environmental sensitivity in *Daphnia Magna*. In *Drosophila melanogaster*, both positive and negative genetic correlations between macro- and micro-environmental sensitivities have been reported[27, 28]. In farm animals, the genetic correlation between macro- and micro-environmental sensitivities has not been investigated. In addition, to date, no suitable method was available to estimate such correlations in outbred populations and the requirements in terms of number of families and family sizes of the designs necessary to estimate genetic parameters of macro- and micro-environmental sensitivities were unknown.

The objectives of this study were to extend the double hierarchical generalized linear model (DHGLM) of Rönnegård et al.[19] with a reaction norm model to estimate genetic variance for macro- and micro-environmental sensitivities and to investigate bias and precision of estimated variance components in designs resembling dairy cattle populations. In addition, we evaluated Akaike’s information criterion (AIC) using approximated h-likelihood to select the best fitting model and studied situations in which true genetic and best fitting statistical models differed. Finally, we applied the model to lactation milk yield in dairy cattle.

## Methods

In this section, we describe the assumed quantitative genetic model underlying genetic variance for macro- and micro-environmental sensitivities, the statistical model to estimate genetic variance in macro- and micro-environmental sensitivities, and the AIC based on an approximation of the h-likelihood that can be used to select the best fitting model. Finally, we present the simulation used to test the statistical model and describe the evaluated scenarios as well as an application to milk yield in dairy cattle.

### Quantitative genetic model

*P*is the phenotype,

*μ*is the population mean for the phenotype,

*x*is the environmental parameter (continuous or discrete) that cause genotypes to respond differently

*, A*

_{ int }

*, A*

_{ sl }and

*A*

_{ v }are the additive genetic or breeding values for the intercept, for the slope of the reaction norm (= macro-environmental sensitivity), and for the environmental variance (= micro-environmental sensitivity), respectively,${\sigma}_{E}^{2}$ is the environmental variance of the exponential model and

*ϵ*is a scaled environmental deviation with variance one. Note that the average environmental variance is$\overline{{\sigma}_{E}^{2}}={\sigma}_{E}^{2}\mathrm{exp}\left(0.5{\sigma}_{{A}_{v}}^{2}\right)$[3]. The additive genetic values

*A*

_{ int }

*, A*

_{ sl }and

*A*

_{ v }have variance

### Statistical model

*y*

_{ i }=

*P*

_{ i }= phenotypic observation of animal

*i*; and a Gamma GLM for the residual variance

*φ*

_{ i }, where$\phantom{\rule{0.25em}{0ex}}{\phi}_{i}\phantom{\rule{0.5em}{0ex}}=\phantom{\rule{0.5em}{0ex}}\mathrm{E}\left(\frac{{\widehat{e}}_{i}^{2}}{1-{h}_{i}}\right),\phantom{\rule{0.5em}{0ex}}{\widehat{e}}_{i}^{2}$ is the squared estimated residual of

*y*

_{ i }and

*h*

_{ i }is the diagonal element of the hat-matrix of

**y**corresponding to

*y*

_{ i }[30]. Felleki et al.[31] extended the model of[19] by allowing estimation of the genetic covariance (${\sigma}_{a,{a}_{v}}$) between the genetic effect in the mean model and the genetic effect in the variance model (Equation (2)) in the absence of a linear reaction norm model and reported on computational details and a theoretical assessment of this algorithm. In brief, equivalent to using a Gamma GLM, linearizing$\frac{{\widehat{e}}_{i}^{2}}{1-{h}_{i}}$ around its current fitted value results in vector

**ψ**and yields the following bivariate linear model:

where **X(X**_{
v
}**)** is the incidence matrix for fixed effects for **y** (**ψ**), **b (b**_{
v
}**)** is the vector with solutions for fixed effects for **y**(**ψ**), **Z(Z**_{
v
}**)** is the incidence matrix for relating observations of **y**(**ψ**) to the additive genetic values **a (a**_{
v
}**)** for phenotype (environmental variance)$\left[\begin{array}{c}\hfill \mathbf{a}\hfill \\ \hfill {\mathbf{a}}_{\mathbf{v}}\hfill \end{array}\right]~\mathrm{N}\left(\mathbf{0},\left[\begin{array}{cc}\hfill {\mathit{\sigma}}_{\mathit{a}}^{\mathit{2}}\hfill & \hfill {\mathit{\sigma}}_{\mathit{a}\mathit{,}{\mathit{a}}_{\mathit{v}}}\hfill \\ \hfill {\mathit{\sigma}}_{\mathit{a}\mathit{,}{\mathit{a}}_{\mathit{v}}}\hfill & \hfill {\mathit{\sigma}}_{{\mathit{a}}_{\mathit{v}}}^{\mathit{2}}\hfill \end{array}\right]\otimes \mathbf{A}\right)$, where **A** is the additive genetic or numerator relationship matrix and${\sigma}_{a}^{2}\phantom{\rule{0.5em}{0ex}}\left({\sigma}_{{a}_{v}}^{2}\right)$ is the genetic variance of **a**(**a**_{
v
}). The residuals of **y** (**e**) and **ψ** (**e**_{
v
}) are assumed to be independent normally distributed$\left[\begin{array}{c}\hfill \mathbf{e}\hfill \\ \hfill {\mathbf{e}}_{\mathbf{v}}\hfill \end{array}\right]~\mathrm{N}\left(\begin{array}{c}\hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill \end{array},\left[\begin{array}{cc}\hfill {\mathbf{W}}^{-1}{\sigma}_{\mathit{\in}}^{2}\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill {\mathbf{W}}_{\mathbf{v}}^{-1}{\sigma}_{{\mathit{\in}}_{v}}^{2}\hfill \end{array}\right]\phantom{\rule{0.2em}{0ex}}\right)$,$\mathbf{W}=\mathrm{diag}{\left(\widehat{\mathbf{\psi}}\right)}^{-1}$ and${\mathbf{W}}_{\mathbf{v}}=\mathrm{diag}\left(\frac{1-\mathbf{h}}{2}\right)$**,**${\sigma}_{\mathit{\epsilon}}^{2}$ and${\sigma}_{{\mathit{\epsilon}}_{v}}^{2}$ are scaling variances, which are expected to be equal to 1, because **W** and **W**_{
v
} already contain the reciprocals of the residual variances per observation. Estimating${\sigma}_{\mathit{\epsilon}}^{2}$ and${\sigma}_{{\mathit{\epsilon}}_{v}}^{2}$ allows more flexibility. An iterative estimation procedure is required to obtain estimates for all parameters because of the dependence of the model for **y** on the results of the model for **ψ** and vice versa. The vector **ψ** and the diagonals of **W** and **W**_{
v
} are updated at each iteration until convergence. The model in Equation (2) can be considered as combining the DHGLM method[19] with the iterative method of Mulder et al.[32].

**y**in Equation (2) was extended with a linear reaction norm model to estimate genetic variance in macro-environmental sensitivity, whereas the residual variance model using

**ψ**estimates the genetic variance in micro-environmental sensitivity or environmental variance. Instead of using an animal model, we used a sire model because implementation of the animal model DHGLM resulted in severely biased estimates of variance components with single observations per animal. As an indication of poor model fit, implementation of the animal model DHGLM gave a lower adjusted profile h-likelihood (see next section “Model selection”) than an animal model without heterogeneity of residual variance. Sire models are commonly used in combination with reaction norm models[7, 8, 33]. In the reaction norm model, we assumed that the environmental parameter

*x*(see Equation (1)) is known without error and does not need to be estimated from the data. The resulting combined macro–micro environmental sensitivity sire model can be formulated as:

where$\phantom{\rule{0.25em}{0ex}}{\mathbf{\psi}}_{\mathbf{s}}$ is a vector with linearized values of transformed squared residuals (see calculation in Appendix), **Z**_{
s
} and **Z**_{
sv
} are respectively the incidence matrices for the sire effects for the intercept of the reaction norm and for the environmental variance, **Z**_{
x
} is the matrix with the environmental parameter *x* as a covariate for the sire effects for the slope of the reaction norm, and **s**_{
int
}, **s**_{
sl
} and **s**_{
v
} are the vectors with the estimated sire effects for intercept, slope, and environmental variance. The sire effects **s**_{
int
}, **s**_{
sl
}, and **s**_{
v
} are assumed trivariate normally distributed$\mathrm{N}\left(\mathbf{0},\frac{1}{4}\mathbf{G}\otimes \mathbf{A}\right)$, assuming that sire (co)variances are a quarter of the additive genetic variance. The residuals of **y (e**_{
s
}**)** and **ψ**_{
s
}(**e**_{
sv
}) are assumed to be independent normally distributed, because cov(*e*, *e*^{2}) = 0, when *e* is normally distributed ($\left[\begin{array}{c}\hfill {\mathbf{e}}_{\mathbf{s}}\hfill \\ \hfill {\mathbf{e}}_{\mathbf{sv}}\hfill \end{array}\right]~\mathrm{N}\left(\begin{array}{c}\hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill \end{array},\left[\begin{array}{cc}\hfill {\mathbf{W}}_{\mathbf{s}}^{-1}{\sigma}_{{\mathit{\in}}_{s}}^{2}\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill {\mathbf{W}}_{\mathbf{sv}}^{-1}{\sigma}_{{\mathit{\in}}_{\mathit{sv}}}^{2}\hfill \end{array}\right]\phantom{\rule{0.2em}{0ex}}\right)$), where${\sigma}_{{\mathit{\in}}_{s}}^{2}$ and${\sigma}_{{\mathit{\in}}_{\mathit{sv}}}^{2}$ are scaling variances for the residual variances in a sire model. Computations of **ψ**_{
s
} and of the diagonals of **W**_{
s
} and **W**_{
sv
} are different due to the use of a sire model compared to an animal model and are explained in the Appendix. The algorithm was implemented by iterating over several runs of ASReml[34]. In each ASReml run, REML-estimates of the variance components were obtained for the current values of **ψ**_{
s
}, **W**_{
s
} and **W**_{
sv
}. The vector **ψ**_{
s
} and the diagonals of **W**_{
s
} and **W**_{
sv
} were updated after each ASReml run. On simulated data, 10 ASReml runs were sufficient to obtain converged parameter estimates. ASReml was run with an option to check whether the variance-covariance matrices were positive definite and were forced to be positive definite[34] if they were not.

### Model selection

*logL*) for the bivariate model in Equation (3) fitted in ASReml, after correcting for the fact that the squared estimated residuals of

**y**are used to compute the response variable

**ψ**

_{ s }(see Appendix for derivation):

*i*, i.e. the

*i*

^{ th }diagonal of

**W**

_{ sv }. Note that Equation (4) can also be used for animal models by replacing elements${w}_{s{v}_{i}}$ and${\sigma}_{{\mathit{\epsilon}}_{\mathit{sv}}}^{2}$ with the corresponding elements of Equation (2). To compare nested and non-nested models, we propose to use APHL in combination with Akaike’s information criterion (AIC):

where t is the number of variance parameters.

### Simulation

*x*was sampled from N(0,1). No fixed effects other than a general mean were simulated (

*μ*= 0) but offspring were randomly allocated to contemporary groups with size 10 to investigate the effect of fixed effects on the precision of genetic parameters. Additive genetic effects were sampled from N(

**0**,

**G**⊗

**A**). Different values for genetic variances and genetic correlations were used for simulation. Table 1 gives an overview of the parameter values used in the simulation, both the default and alternative values. In each scenario, only one parameter was varied at a time, i.e. other parameters were kept at their default values. For each set of parameters, 100 replicates were generated and means and standard deviations of estimates were calculated across replicates.

**Default and alternative parameters values used in Monte Carlo simulation**

Parameter | Default | Alternative values |
---|---|---|

${\sigma}_{{A}_{\mathit{int}}}^{2}$ | 0.3 | 0.1, 0.5 |

${\sigma}_{{A}_{\mathit{sl}}}^{2}$ | 0.05 | 0.025, 0.10 |

${\sigma}_{{A}_{v}}^{2}$ | 0.1 | 0.05, 0.2 |

Genetic correlations (see text below Equation (1)) | 0 | 0.5 |

Number of offspring per sire | 100 | 20, 50, 200 |

Number of sires | 100 | 50, 200 |

Number of replicates | 100 |

### Scenarios in which true genetic and statistical models differ

First, we tested situations in which the true genetic and statistical models contained both types of environmental sensitivities for different sets of genetic parameters and designs. In addition, we investigated the power of AIC to select the correct model given the true genetic model. Three situations were simulated with default parameters: presence of genetic variance in (i) macro- or (ii) micro-environmental sensitivity or in (iii) both. Four statistical models were applied to these simulated situations: a combined macro–micro environmental sensitivity model “macro–micro”, a macro-environmental sensitivity model “macro”, a micro-environmental sensitivity model “micro” and a simple model “simple” with only one additive genetic effect for the phenotype. The “macro–micro” model and the “micro” model were run for 10 ASReml runs to update the weights and the transformed squared residuals, whereas the “macro” model and the “simple” model required only one ASReml run. AIC (Equation 5) was used to select the best fitting model.

**Scenarios with different combinations of true genetic models and statistical models**

Statistical model | |||
---|---|---|---|

True genetic model | Macro ES | Micro ES | Macro and micro ES |

Macro ES | Not addressed | A | B |

Micro ES | C | Not addressed | D |

Macro and micro ES | E | F | G (Default in this study) |

### Application to milk yield in dairy cattle

The “macro–micro”, “macro”, “micro” and “simple” models were applied to 305-day lactation milk yield of 142 565 first parity Swedish Holstein cows. Lactation yields were calculated based on test-day data, as described in[37] using the test-day interval method[38]. Seasons were defined as three-month periods (January-March, April-June, July-September, October-December). Herd-year-season (HYS) classes with less than five cows were excluded. Data on a total of 762 sires with on average 187 daughters were available and 213 of these sires had more than 100 daughters. In the model for phenotype (**y**), HYS was a fixed effect and in the model for residual variance (**ψ**_{
s
}), HYS was a random effect to allow regression to the mean with few observations per HYS and avoid extreme HYS estimates. The model for **y** included a fixed fourth order polynomial for age at calving and a fixed ninth-order polynomial for lactation length, while the model for **ψ**_{
s
} included a third and a sixth order polynomial for age at calving and lactation length, respectively. The orders of polynomials were determined using a Wald test in univariate models for **y** and **ψ**_{
s
} and were kept the same for “macro–micro”, “macro”, “micro” and “simple” models. Herd-year average milk yield standardized to a mean of zero and a variance of one was used as an environmental parameter for macro-environmental sensitivity. The “macro–micro”, “macro”, “micro” and “simple” models were all run for 50 ASReml runs. After 50 runs, the change in variance components was less than 0.2% of the change between consecutive runs and therefore parameters were considered converged.

## Results

### Variance components and effect of fixed effects

**Means and standard deviations of estimated genetic parameters across 100 replicates when genetic correlations are zero**

Fixed effects | True parameters | Estimated parameters (SD) | Np | |||||||
---|---|---|---|---|---|---|---|---|---|---|

${\sigma}_{{A}_{\mathit{int}}}^{2}$ | ${\sigma}_{{A}_{\mathit{sl}}}^{2}$ | ${\sigma}_{{A}_{v}}^{2}$ | ${\sigma}_{{A}_{\mathit{int}}}^{2}$ | ${\sigma}_{{A}_{\mathit{sl}}}^{2}$ | ${\sigma}_{{A}_{v}}^{2}$ | |||||

No | 0.1 | 0.05 | 0.1 | 0.101 | (0.020) | 0.051 | (0.012) | 0.095 | (0.035) | 0 |

No | 0.3 | 0.05 | 0.1 | 0.315 | (0.053) | 0.054 | (0.015) | 0.107 | (0.046) | 0 |

No | 0.5 | 0.05 | 0.1 | 0.507 | (0.066) | 0.052 | (0.013) | 0.115 | (0.061) | 3 |

No | 0.3 | 0.05 | 0.05 | 0.297 | (0.047) | 0.053 | (0.012) | 0.053 | (0.029) | 4 |

No | 0.3 | 0.05 | 0.2 | 0.308 | (0.049) | 0.053 | (0.014) | 0.186 | (0.065) | 0 |

No | 0.3 | 0.025 | 0.1 | 0.298 | (0.047) | 0.026 | (0.008) | 0.097 | (0.041) | 0 |

No | 0.3 | 0.1 | 0.1 | 0.296 | (0.052) | 0.104 | (0.021) | 0.083 | (0.045) | 1 |

Yes | 0.1 | 0.05 | 0.1 | 0.099 | (0.022) | 0.053 | (0.014) | 0.091 | (0.037) | 1 |

Yes | 0.3 | 0.05 | 0.1 | 0.309 | (0.044) | 0.057 | (0.015) | 0.104 | (0.051) | 0 |

Yes | 0.5 | 0.05 | 0.1 | 0.503 | (0.084) | 0.052 | (0.014) | 0.104 | (0.074) | 10 |

### Effect of genetic correlations

**Means and standard deviations across 100 replicates of estimated genetic parameters when genetic correlations are not zero**

True parameters | Estimated parameters (SD) | ||||||||
---|---|---|---|---|---|---|---|---|---|

${\rho}_{{A}_{\mathit{int}},{A}_{v}}$ | ${\rho}_{{A}_{\mathit{int}},{A}_{\mathit{sl}}}$ | ${\rho}_{{A}_{\mathit{sl}},{A}_{v}}$ | ${\sigma}_{{A}_{\mathit{int}}}^{2}$ | ${\sigma}_{{A}_{\mathit{sl}}}^{2}$ | ${\sigma}_{{A}_{v}}^{2}$ | ${\rho}_{{A}_{\mathit{int}},{A}_{v}}$ | ${\rho}_{{A}_{\mathit{int}},{A}_{\mathit{sl}}}$ | ${\rho}_{{A}_{\mathit{sl}},{A}_{v}}$ | Np |

0 | 0 | 0 | 0.315 | 0.054 | 0.107 | −0.005 | 0.004 | 0.012 | 0 |

(0.053) | (0.015) | (0.046) | (0.165) | (0.155) | (0.249) | ||||

0.5 | 0 | 0 | 0.303 | 0.051 | 0.099 | 0.554 | −0.028 | −0.007 | 1 |

(0.047) | (0.012) | (0.048) | (0.155) | (0.136) | (0.203) | ||||

0 | 0.5 | 0 | 0.303 | 0.053 | 0.094 | −0.010 | 0.508 | −0.014 | 1 |

(0.056) | (0.012) | (0.045) | (0.200) | (0.132) | (0.239) | ||||

0 | 0 | 0.5 | 0.293 | 0.052 | 0.092 | 0.014 | 0.021 | 0.558 | 2 |

(0.045) | (0.013) | (0.034) | (0.185) | (0.147) | (0.208) | ||||

0.5 | 0.5 | 0.5 | 0.301 | 0.053 | 0.089 | 0.537 | 0.517 | 0.530 | 6 |

(0.051) | (0.013) | (0.036) | (0.171) | (0.138) | (0.192) |

### Effect of different designs

**Means and standard deviations across 100 replicates of estimated genetic parameters for different designs**

Estimated parameters (SD) | ||||||||
---|---|---|---|---|---|---|---|---|

NS | NO | ${\sigma}_{{A}_{\mathit{int}}}^{2}$ | ${\sigma}_{{A}_{\mathit{sl}}}^{2}$ | ${\sigma}_{{A}_{v}}^{2}$ | ${\rho}_{{A}_{\mathit{int}},{A}_{v}}$ | ${\rho}_{{A}_{\mathit{int}},{A}_{\mathit{sl}}}$ | ${\rho}_{{A}_{\mathit{sl}},{A}_{v}}$ | Np |

100 | 20 | 0.301 | 0.063 | 0.120 | −0.008 | −0.004 | 0.083 | 51 |

(0.057) | (0.038) | (0.113) | (0.304) | (0.216) | (0.416) | |||

100 | 50 | 0.309 | 0.056 | 0.099 | −0.016 | −0.013 | 0.025 | 8 |

(0.051) | (0.020) | (0.060) | (0.279) | (0.165) | (0.376) | |||

100 | 100 | 0.315 | 0.054 | 0.107 | −0.005 | 0.004 | 0.012 | 0 |

(0.053) | (0.015) | (0.046) | (0.165) | (0.155) | (0.249) | |||

100 | 200 | 0.301 | 0.053 | 0.093 | 0.000 | 0.007 | −0.013 | 0 |

(0.046) | (0.009) | (0.024) | (0.135) | (0.116) | (0.168) | |||

50 | 100 | 0.312 | 0.053 | 0.107 | 0.015 | 0.002 | 0.036 | 5 |

(0.079) | (0.018) | (0.064) | (0.295) | (0.199) | (0.320) | |||

200 | 100 | 0.301 | 0.053 | 0.099 | 0.009 | −0.002 | −0.009 | 0 |

(0.033) | (0.009) | (0.028) | (0.142) | (0.108) | (0.163) |

### Model selection

**The best model selected in 100 replicates according to Akaike’s information criterion and effect of the number of offspring for different true genetic models**

NO | True genetic model | True parameters | Nb of times model selected based on AIC | ||||
---|---|---|---|---|---|---|---|

${\sigma}_{{A}_{v}}^{2}$ | ${\sigma}_{{A}_{\mathit{sl}}}^{2}$ | Macro–micro | Macro | Micro | Simple | ||

100 | Macro–micro | 0.1 | 0.05 | 99 | 1 | 0 | 0 |

Macro | 0 | 0.05 | 5 | 95 | 0 | 0 | |

Micro | 0.1 | 0 | 3 | 1 | 94 | 2 | |

Simple | 0 | 0 | 1 | 4 | 5 | 90 | |

50 | Macro–micro | 0.1 | 0.05 | 63 | 27 | 8 | 2 |

Macro | 0 | 0.05 | 7 | 87 | 1 | 5 | |

Micro | 0.1 | 0 | 4 | 1 | 71 | 24 | |

Simple | 0 | 0 | 1 | 1 | 3 | 95 |

### Scenarios in which true genetic and statistical models differ

**Means and standard deviations across 100 replicates of estimated genetic parameters when true and statistical models differ**

Scenario | True model | True parameters | Statistical model | Estimated parameters (SD) | ||||
---|---|---|---|---|---|---|---|---|

${\sigma}_{{A}_{\mathit{sl}}}^{2}$ | ${\sigma}_{{A}_{v}}^{2}$ | ${\sigma}_{{A}_{\mathit{int}}}^{2}$ | ${\sigma}_{{A}_{\mathit{sl}}}^{2}$ | ${\sigma}_{{A}_{v}}^{2}$ | Np | |||

A | Macro | 0.025 | 0 | Micro | 0.303 | 0.010 | 51 | |

(0.045) | (0.020) | |||||||

Macro | 0.05 | 0 | Micro | 0.303 | 0.010 | 48 | ||

(0.044) | (0.013) | |||||||

Macro | 0.1 | 0 | Micro | 0.308 | 0.018 | 37 | ||

(0.050) | (0.019) | |||||||

B | Macro | 0.025 | 0 | Macro–micro | 0.296 | 0.027 | 0.012 | 60 |

(0.046) | (0.009) | (0.011) | ||||||

Macro | 0.05 | 0 | Macro–micro | 0.298 | 0.052 | 0.013 | 51 | |

(0.048) | (0.011) | (0.013) | ||||||

Macro | 0.1 | 0 | Macro–micro | 0.307 | 0.091 | 0.012 | 54 | |

(0.051) | (0.020) | (0.013) | ||||||

C | Micro | 0 | 0.05 | Macro | 0.296 | 0.002 | 0 | |

(0.041) | (0.003) | |||||||

Micro | 0 | 0.1 | Macro | 0.292 | 0.003 | 0 | ||

(0.039) | (0.004) | |||||||

Micro | 0 | 0.2 | Macro | 0.289 | 0.002 | 0 | ||

(0.043) | (0.003) | |||||||

D | Micro | 0 | 0.05 | Macro–micro | 0.303 | 0.003 | 0.053 | 71 |

(0.050) | (0.003) | (0.024) | ||||||

Micro | 0 | 0.1 | Macro–micro | 0.300 | 0.003 | 0.084 | 60 | |

(0.047) | (0.003) | (0.034) | ||||||

Micro | 0 | 0.2 | Macro–micro | 0.304 | 0.002 | 0.143 | 66 | |

(0.049) | (0.003) | (0.058) | ||||||

E | Macro–micro | 0.05 | 0.05 | Macro | 0.297 | 0.055 | 0 | |

(0.047) | (0.014) | |||||||

Macro–micro | 0.05 | 0.1 | Macro | 0.298 | 0.053 | 0 | ||

(0.052) | (0.012) | |||||||

Macro–micro | 0.05 | 0.2 | Macro | 0.292 | 0.053 | 0 | ||

(0.049) | (0.011) | |||||||

Macro–micro | 0.025 | 0.1 | Macro | 0.306 | 0.026 | 0 | ||

(0.050) | (0.008) | |||||||

Macro–micro | 0.1 | 0.1 | Macro | 0.298 | 0.106 | 0 | ||

(0.045) | (0.022) | |||||||

F | Macro–micro | 0.05 | 0.05 | Micro | 0.299 | 0.050 | 3 | |

(0.051) | (0.031) | |||||||

Macro–micro | 0.05 | 0.1 | Micro | 0.297 | 0.095 | 1 | ||

(0.047) | (0.042) | |||||||

Macro–micro | 0.05 | 0.2 | Micro | 0.304 | 0.191 | 0 | ||

(0.045) | (0.070) | |||||||

Macro–micro | 0.025 | 0.1 | Micro | 0.299 | 0.105 | 0 | ||

(0.049) | (0.043) | |||||||

Macro–micro | 0.1 | 0.1 | Micro | 0.298 | 0.103 | 1 | ||

(0.045) | (0.045) |

### Application to milk yield in dairy cattle

**Estimated genetic parameters for macro- and micro-environmental sensitivity of milk yield in dairy cattle**

Parameter | Macro–micro | Macro | Micro | Simple | ||||
---|---|---|---|---|---|---|---|---|

Estimate | SE | Estimate | SE | Estimate | SE | Estimate | SE | |

${\sigma}_{{A}_{\mathit{int}}}^{2}$ | 420 800 | 27960 | 420 400 | 28 004 | 416 800 | 27 696 | 416 000 | 27 692 |

${\sigma}_{{A}_{\mathit{sl}}}^{2}$ | 11 096 | 2288 | 11 116 | 2320 | ||||

${\sigma}_{{A}_{v}}^{2}$ | 0.043 | 0.008 | 0.042 | 0.008 | ||||

${\rho}_{{A}_{\mathit{int}},{A}_{\mathit{sl}}}$ | 0.808 | 0.062 | 0.812 | 0.063 | ||||

${\rho}_{{A}_{\mathit{int}},{A}_{v}}$ | 0.627 | 0.073 | 0.608 | 0.0751 | ||||

${\rho}_{{A}_{\mathit{sl}},{A}_{v}}$ | 0.765 | 0.098 | ||||||

APHL | 193 704 | 194 179 | 193 692 | 202 832 | ||||

AIC | 193 722 | 194 191 | 193 704 | 202 840 |

## Discussion

### Model and design

In this study, we developed a model to estimate genetic variances for macro- and micro-environmental sensitivities. The model is an extension of the DHGLM as presented by Rönnegård et al.[19]. Here, we combined a linear reaction norm model to estimate genetic variance for macro-environmental sensitivity with the DHGLM to estimate genetic variance for micro-environmental sensitivity. The animal model in Equation (2) was adapted to a sire model because the animal model produced highly biased estimated variance components because of the high dependence of the estimated breeding values and residuals on the variance ratio used in the mixed model equations, which would differ for each animal. Felleki et al.[31] also reported the presence of bias in variance components with few repeated observations per animal but the bias decreased as the number of repeated observations per animal increased. Furthermore, an animal model with heterogeneous residual variance (DHGLM implementation) gave a poorer adjusted profile h-likelihood than an animal model without homogenous residual variance, which indicates that the former did not produce a better fit than the latter in a scenario that included both macro and micro-environmental sensitivities. Therefore, we decided to use a sire model implementation, because it is more robust than the animal model implementation with DHGLM when animals only have a single observation. Furthermore, sire models are commonly used for reaction norm models, because of their substantially lower computational burden compared to animal models and genetic information about environmental sensitivity typically comes from paternal half-sibs that perform in different environments. As far as we know, this is the first time that a model to estimate genetic variance in macro- and micro-environmental sensitivities suitable for outbred animal populations is presented.

Monte Carlo simulation was used to evaluate bias and precision of estimated genetic parameters. Genetic parameters were unbiased in most situations. The precision was not very high, especially of estimates for${\sigma}_{{A}_{v}}^{2}$, as indicated by the high standard deviation of estimates across replicates, particularly in designs with small families. Designs with at least 100 sire families, each with at least 100 offspring, are required to have sufficient precision. Presence of fixed effects, such as contemporary group effects, would increase the required number of sire families and the number of offspring per sire family. These results are in agreement with standard error and power calculations reported by Hill[39], Mulder[40] and Hill and Mulder[20] with respect to estimation of genetic variance for micro-environmental sensitivity or for environmental variance. For instance, Hill and Mulder[20] derived that the optimal family size to estimate genetic variance for environmental variance with a family design is approximately 2/*γ*^{2}, where *γ*^{2} is the square of the coefficient of variation of the within-family variance. Thus, the optimal family size for half-sibs is 137 when${\sigma}_{{A}_{v}}^{2}=0.10$. The only study providing Monte Carlo results for the DHGLM to estimate genetic variance for micro-environmental sensitivity is by Rönnegård et al.[19]. They considered a design with clones, which is more powerful and leads to lower standard deviations of estimates than those in our study.

With respect to macro-environmental sensitivity, the magnitude of the standard deviations of the${\sigma}_{{A}_{\mathit{sl}}}^{2}$ estimates across replicates was similar to that reported by Calus et al.[36], but larger than those reported by Lillehammer et al.[41], which is explained by the fact that the latter authors simulated more sire families, i.e., 1000 sire families with 100 offspring each. In general, the number of required offspring per family is lower for macro-environmental sensitivity than for micro-environmental sensitivity, as indicated by the lower standard deviations of estimates of genetic variance for the former. This is in agreement with a previous study by Mulder[40], which showed that the power to detect G x E interactions between two environments is greater than the power to detect genetic variance for environmental variance or micro-environmental sensitivity. In the present study, we assumed that the environmental parameter *x* used for the reaction norm was known without error. This will be true in some cases, e.g. when using temperature or rainfall or other herd characteristics[42]. In other cases, an estimated herd mean is used as the environmental parameter, which is estimated from the data[7, 8, 33]. Calus et al.[36] showed that genetic variance in macro-environmental sensitivity (${\sigma}_{{A}_{\mathit{sl}}}^{2}$) was severely underestimated when the environmental parameter was estimated from the data. This may have also led to the genetic variance in macro-environmental sensitivity to be underestimated in our application to milk yield. Su et al.[43] reported that a Bayesian approach that estimates simultaneously the herd mean and the reaction norm parameters leads to unbiased estimates of${\sigma}_{{A}_{\mathit{sl}}}^{2}$.

Lillehammer et al.[41] showed that sire models gave upward biased estimates of${\sigma}_{{A}_{\mathit{sl}}}^{2}$ when heterogeneity of residual variance was ignored, because the unexplained genetic variance in the reaction norm parameters becomes part of the residual variance when using a sire model. In our study, the bias in estimates of${\sigma}_{{A}_{\mathit{sl}}}^{2}$ was smaller than in Lillehammer et al.[41]. Lillehammer et al.[41] proposed including a dummy animal effect in the model to account for the residual three-quarters of the genetic variance that is not accounted for by the sire effect. This solution was also tested in our model but gave severely biased variance components because of the high dependency of estimated breeding values and residuals on variance ratios that were used in the mixed model equations, which differ by animal when considering heterogeneity of residual variance.

The algorithm developed in this study allowed estimating genetic correlations between the different genetic effects. Standard deviations of estimated genetic correlations were large, especially those of the genetic correlation between macro- and micro-environmental sensitivities (${\rho}_{{A}_{\mathit{sl}},{A}_{v}}$). This large magnitude of the standard deviation of estimates of the genetic correlation was expected considering that the genetic correlation between macro- and micro-environmental sensitivities is mainly based on paternal half-sib information and that both traits have low heritability. Using the equation in Robertson[44] and assuming a heritability of 0.05 for both traits, the standard error is approximately 0.26 when the true genetic correlation between the traits is zero, which is close to the value found here, i.e. 0.25 (Table 4). To increase the precision of estimates of the genetic correlation between macro- and micro-environmental sensitivities, designs with a larger number of families and larger family sizes are required. The application to milk yield data in dairy cattle shows that it should be possible to estimate genetic correlations with standard errors between 0.06 and 0.10, since in most countries datasets with at least 100 bulls each with 100 daughters are easily obtained. Full-sib families or clones would also reduce standard errors of estimates of genetic correlations in comparison to half-sibs.

Here, we showed that the adjusted profile h-likelihood (APHL) can be approximated from REML-output and used in combination with AIC to provide an efficient model selection tool. In addition, we showed that biases in genetic parameters were relatively small when statistical and true models differed. Both results are re-assuring that these models can discriminate between macro- and micro-environmental sensitivities and that the true model of environmental sensitivity can be elucidated using AIC. We also found that the Bayesian information criterion (BIC) was too conservative and favoured the simpler model too often (results not shown). AIC has the advantage that it can be used independent of the order of the fitted models, whereas the likelihood ratio test requires a hierarchical structure such as in a forward selection scheme[45].

### Improving biological understanding of environmental sensitivity

The proposed model can contribute to better understand the genetic architecture of environmental sensitivity, e.g. whether macro- and micro-environmental sensitivities are genetically related. Most studies in plants and laboratory species indicate that macro- and micro-environmental sensitivities are weakly correlated[22–24, 26–28]. This seems to indicate that selection on one type of environmental sensitivity will hardly affect the other. The first application of the model on milk yield data in dairy cattle, revealed a high genetic correlation between macro- and micro-environmental sensitivities (0.76), suggesting that selection on one type of environmental sensitivity will also affect the other in the same direction. Generally, little is known about these relationships in livestock. Knowledge about these genetic correlations could be used to optimize selection strategies for environmental sensitivity.

In this study, we assumed a linear reaction norm model, but reaction norms can also be non-linear[8]. The model presented here can easily be extended to higher-order polynomials. Furthermore, the genetic basis of micro-environmental sensitivity may not be the same along an environmental gradient and the model for residual variance or micro-environmental sensitivity in Equation (3) could be extended to contain a reaction norm with a known environmental gradient. For instance, in stressful environments, there might be more genetic variance for micro-environmental sensitivity than in less stressful environments. Furthermore, G × E interactions often exist between categorical environments and are often analysed with character state or multivariate models[2]. Character state models do not explicitly estimate breeding values for macro-environmental sensitivity, but these breeding values could be back-calculated by using covariance functions when the environmental parameter responsible for G × E interactions is identified because reaction norm models and character state models are interchangeable[2]. In the case of two environments, a reaction norm model with a dummy environmental parameter with values 0 and 1 would yield results that are identical to a bivariate character state model. Multivariate versions of the DHGLM or reaction norm models with dummy environmental variables could be used to simultaneously investigate macro- and micro-environmental sensitivities when environments are discrete.

### Application to breeding

Taking macro- and micro-environmental sensitivities into consideration is highly relevant in animal breeding. Due to the high level of globalisation in animal breeding programs, it is necessary to breed animals that can perform well in a wide range of environments. Therefore, it may be important to select animals that have limited environmental sensitivity, especially for environments with a higher risk of environmental disturbances. Reduction in environmental sensitivity increases the predictability of performance and reduces risk for farmers[46]. Furthermore, reduction in micro-environmental sensitivity will increase the uniformity of animal products[47], which is a general goal. In plant breeding, application of a model for macro- and micro-environmental sensitivities is also highly relevant since G × E interactions are generally very strong and uniformity of crops is very important. Recent papers by Ordas et al.[48], Makumburage and Stapleton[49] and Kliebenstein[50] show that there is an interest for increased uniformity in plants. Economic values could be derived for micro-environmental sensitivity[47]. The economic value of macro-environmental sensitivity can be determined as a function of the importance of environments along the environmental gradient. Based on Mulder et al.[3, 47], progeny testing schemes are more efficient than sib testing schemes to reduce micro-environmental sensitivity since it behaves as a trait with a small heritability. Genomic selection could be an alternative selection strategy with sufficient accuracy and shorter generation intervals.

## Conclusions

In this study, a model was developed to estimate genetic parameters of macro- and micro-environmental sensitivities, combining a reaction norm model with a double hierarchical generalized linear model within a REML framework. Simulations showed that the genetic parameters obtained were mostly unbiased, but designs with at least 100 sires, each with 100 half-sib offspring, were required to estimate genetic parameters with sufficient precision. Using AIC, the true genetic model was selected as the best statistical model in at least 90% of replicates when the number of offspring per sire was 100. Application of the model to milk yield data in dairy cattle showed that both types of environmental sensitivity existed. Our model and AIC based on h-likelihood can be used to increase our understanding of the genetic control of environmental sensitivity in livestock populations but more research is needed to test the model in a wider range of situations.

## Appendix

### DHGLM algorithm for a sire model

The original DHGLM algorithm of Rönnegård et al.[19] was developed for an animal model. Here we describe the estimation algorithm for the sire model used in the current paper, including a few adjustments of the algorithm in Rönnegård et al.[19] to correct for the fact that the residual variance in a sire model (without permanent environmental effects, e.g. with animals with a single observation) contains three quarters of the additive genetic variance in addition to the environmental variance. The adjustments in the algorithm are as follows: adjustment of the squared residuals (**y**_{v}) accounting for the fact that the residual variance in a sire model includes three quarters of the additive genetic variance, use of average residual variance${\widehat{\sigma}}_{{e}_{s}}^{2}$ to calculate **ψ**_{
s
} (instead of predicted individual values), and computations of the diagonal weight matrices **W**_{
s
} and **W**_{
sv
}. These adjustments resulted in a computationally robust algorithm with small or no bias (as presented in Results).

**ψ**

_{ s }for the bivariate sire model in Equation (3), first we calculate${y}_{{v}_{i}}$ as:

*y*

_{v}in[19], except for the multiplication by the ratio of the average estimated residual variance in a sire model ($\overline{{\widehat{\sigma}}_{{e}_{s}}^{2}}={\sigma}_{{\mathit{\in}}_{s}}^{2}/\overline{{w}_{s}}$, where$\overline{{w}_{s}}=\mathrm{tr}\left({\mathbf{W}}_{\mathbf{s}}\right)/n$and

*n*is the total number of records) and the average residual variance in an animal model, which is calculated as$\overline{{\widehat{\sigma}}_{{e}_{a}}^{2}}=\overline{{\widehat{\sigma}}_{{e}_{s}}^{2}}-\frac{3}{4}{\widehat{\sigma}}_{{A}_{\mathit{int}}}^{2}$, ($\phantom{\rule{0.25em}{0ex}}{\widehat{\sigma}}_{{A}_{\mathit{int}}}^{2}=4\phantom{\rule{0.25em}{0ex}}{\widehat{\sigma}}_{{s}_{\mathit{int}}}^{2}$). Because we use a log link function,${y}_{{v}_{i}}$ is linearized as:

**W**

_{s}are the reciprocals of${\widehat{\mathbf{\psi}}}_{\mathbf{s}}$ (i.e.${W}_{s}{}_{{}_{i}}=1/{\widehat{\mathit{\psi}}}_{{s}_{i}}$), which is the vector of predicted residual variances for each observation based on the previous iteration and is calculated as:

**W**

_{ sv }are the reciprocals of the residual variance of

**ψ**

_{ s }and were calculated as:

since the true residuals are assumed normal and the squared true residuals are therefore Gamma distributed with a variance of$2{\left({\widehat{\sigma}}_{{e}_{s}}^{2}\right)}^{2}$.

- 1.
Run model on

**y**in Equation (3) with homogeneous residual variance. - 2.
Calculate

**ψ**_{ s },**W**_{ s },**W**_{ sv }, where ${\mathbf{W}}_{\mathbf{s}}=\mathrm{diag}\left({\widehat{\sigma}}_{e}^{2}\right)$ in iteration 1. - 3.
Run bivariate model in Equation (3).

- 4.
Update

**ψ**_{ s },**W**_{ s },**W**_{ sv } - 5.
Iterate steps 3 till 4 until convergence

**Approximation of the adjusted profile h-likelihood**

**H**is the Hessian of the h-likelihood and τ is the vector of all fixed and random effects both in the mean and variance parts of the model. For the model in Equation (3), minus two times the h-likelihood (− 2

*h*) is:

*l*is log-likelihood,${\widehat{\sigma}}_{{e}_{i}}^{2}$ (${\widehat{\sigma}}_{{e}_{i}}^{2}={\widehat{\sigma}}_{{\mathit{\in}}_{s}}^{2}/{w}_{{s}_{i}}$, with${w}_{{s}_{i}}$ being the

*i*

^{ th }diagonal of

**W**

_{ s }) is the residual variance of observation

*i*,${\widehat{e}}_{i}^{2}$ is the estimated squared residual of observation

*i*,$\tilde{\mathbf{G}}$ is the covariance matrix of all random effects ($\tilde{\mathbf{G}}=\frac{1}{4}\mathbf{G}\otimes \mathbf{A}$),

**s**is a vector of all random sire effects ($\tilde{\mathbf{s}}\text{'}=\left[\begin{array}{ccc}\hfill {\mathbf{s}}_{\mathbf{int}}\hfill & \hfill {\mathbf{s}}_{\mathbf{sl}}\hfill & \hfill {\mathbf{s}}_{\mathbf{v}}\hfill \end{array}\right]$) and${w}_{{s}_{i}}$ is the weight for the mean model for observation

*i*. The minus two log REML likelihood (

*logL*) from the bivariate model in Equation (3) is:

*n*residual variances come from the first part of the bivariate model (

**y**) and the next

*k*residual variances come from the second part of the bivariate model (

**ψ**

_{ s }) (

*k*=

*n*= number of records), and

**C**is the Hessian of the bivariate model (i.e. left-hand-side of the mixed model equations). Because log (det(

**C**)) is a reasonable approximation of log (det (

**H**))[31], we can approximate

*APHL*as given in Equation (4) in the main text:

Thus, the REML likelihood for the bivariate model is corrected for the fact that the squared residuals are used as “observations” in the bivariate model. Note that equation (14) can also be used for animal models by replacing the elements${w}_{s{v}_{i}}$ and${\sigma}_{{\mathit{\epsilon}}_{\mathit{sv}}}^{2}$ with the corresponding elements of Equation (2) or for models with more or fewer random effects in the mean and variance model.

## Declarations

### Acknowledgements

We would like to thank Jack Windig, Majbritt Felleki, Jessica Franzén, Jorge Urioste and Susan Wijga for useful discussions about this research. This project was financed by the RobustMilk project, which is financially supported by the European Commission under the Seventh Research Framework Programme, Grant Agreement KBBE-211708. The content of this paper is the sole responsibility of the authors, and it does not necessarily represent the views of the Commission or its services.

## Authors’ Affiliations

## References

- Falconer DS, Mackay TFC: Introduction to quantitative genetics. 1996, Essex: Pearson Education Limited, 4Google Scholar
- De Jong G, Bijma P: Selection and phenotypic plasticity in evolutionary biology and animal breeding. Livest Prod Sci. 2002, 78: 195-214. 10.1016/S0301-6226(02)00096-9.View ArticleGoogle Scholar
- Mulder HA, Bijma P, Hill WG: Prediction of breeding values and selection responses with genetic heterogeneity of environmental variance. Genetics. 2007, 175: 1895-1910. 10.1534/genetics.106.063743.PubMed CentralView ArticlePubMedGoogle Scholar
- Meyer K: Factor-analytic models for genotype x environment type problems and structured covariance matrices. Genet Sel Evol. 2009, 41: 21-10.1186/1297-9686-41-21.PubMed CentralView ArticlePubMedGoogle Scholar
- Piepho HP, Mohring J, Melchinger AE, Buchse A: BLUP for phenotypic selection in plant breeding and variety testing. Euphytica. 2008, 161: 209-228. 10.1007/s10681-007-9449-8.View ArticleGoogle Scholar
- Scheiner SM: Selection experiments and the study of phenotypic plasticity. J Evol Biol. 2002, 15: 889-898. 10.1046/j.1420-9101.2002.00468.x.View ArticleGoogle Scholar
- Kolmodin R, Strandberg E, Madsen P, Jensen J, Jorjani H: Genotype by environment interaction in Nordic dairy cattle studied using reaction norms. Acta Agric Scand A. 2002, 52: 11-24.Google Scholar
- Calus MPL, Veerkamp RF: Estimation of environmental sensitivity of genetic merit for milk production traits using a random regression model. J Dairy Sci. 2003, 86: 3756-3764. 10.3168/jds.S0022-0302(03)73982-4.View ArticlePubMedGoogle Scholar
- Bergsma R, Hermesch S: Exploring breeding opportunities for reduced thermal sensitivity of feed intake in the lactating sow. J Anim Sci. 2012, 90: 85-98. 10.2527/jas.2011-4021.View ArticlePubMedGoogle Scholar
- Bradshaw AD: Evolutionary significance of phenotypic plasticity in plants. Adv Genet. 1965, 13: 115-155.View ArticleGoogle Scholar
- Scheiner SM: Genetics and evolution of phenotypic plasticity. Annu Rev Ecol Syst. 1993, 24: 35-68. 10.1146/annurev.es.24.110193.000343.View ArticleGoogle Scholar
- Nussey DH, Postma E, Gienapp P, Visser ME: Selection on heritable phenotypic plasticity in a wild bird population. Science. 2005, 310: 304-306. 10.1126/science.1117004.View ArticlePubMedGoogle Scholar
- Windig JJ: Reaction norms and the genetic basis of phenotypic plasticity in the wing pattern of the butterfly Bicyclus anynana. J Evol Biol. 1994, 7: 665-695. 10.1046/j.1420-9101.1994.7060665.x.View ArticleGoogle Scholar
- Waddington CH: Canalization of development and the inheritance of acquired characters. Nature. 1942, 150: 563-565. 10.1038/150563a0.View ArticleGoogle Scholar
- Flatt T: The evolutionary genetics of canalization. Q Rev Biol. 2005, 80: 287-316. 10.1086/432265.View ArticlePubMedGoogle Scholar
- Debat V, David P: Mapping phenotypes: canalization, plasticity and developmental stability. Trends Ecol Evol. 2001, 16: 555-561. 10.1016/S0169-5347(01)02266-2.View ArticleGoogle Scholar
- Sorensen D, Waagepetersen R: Normal linear models with genetically structured residual variance heterogeneity: a case study. Genet Res. 2003, 82: 207-222. 10.1017/S0016672303006426.View ArticlePubMedGoogle Scholar
- Ibanez-Escriche N, Sorensen D, Waagepetersen R, Blasco A: Selection for environmental variation: a statistical analysis and power calculations to detect reponse. Genetics. 2008, 180: 2209-2226. 10.1534/genetics.108.091678.PubMed CentralView ArticlePubMedGoogle Scholar
- Rönnegård L, Felleki M, Fikse F, Mulder HA, Strandberg E: Genetic heterogeneity of residual variance: estimation of variance components using double hierarchical generalized linear models. Genet Sel Evol. 2010, 42: 8-10.1186/1297-9686-42-8.PubMed CentralView ArticlePubMedGoogle Scholar
- Hill WG, Mulder HA: Genetic analysis of environmental variation. Genet Res. 2010, 92: 381-395. 10.1017/S0016672310000546.View ArticleGoogle Scholar
- Jinks JL, Pooni HS: The genetic basis of environmental sensitivity. Proceedings of the Second International Conference on Quantitative Genetics: 31 May - 5 June 1987; Raleigh. Edited by: Weir BS, Eisen EJ, Goodman MM, Namkoong G. 1988, North Carolina State University, 505-522.Google Scholar
- Perkins JM, Jinks JL: The assessment and specificity of environmental and genotype-environmental components of variability. Heredity. 1973, 30: 111-126. 10.1038/hdy.1973.16.View ArticlePubMedGoogle Scholar
- Wu RL: Genetic control of macro- and micro-environmental sensitivities in Populus. Theor Appl Genet. 1997, 94: 104-114. 10.1007/s001220050388.View ArticlePubMedGoogle Scholar
- Wu RL, O'Malley DM: Nonlinear genotypic response to macro- and micro-environments. Theor Appl Genet. 1998, 96: 669-675. 10.1007/s001220050787.View ArticleGoogle Scholar
- Kraakman ATW, Niks RE, Van den Berg PMMM, Stam P, Van Eeuwijk FA: Linkage disequilibrium mapping of yield and yield stability in modern spring barley cultivars. Genetics. 2004, 168: 435-446. 10.1534/genetics.104.026831.PubMed CentralView ArticlePubMedGoogle Scholar
- Yampolsky LY, Scheiner SM: Developmental noise, phenotypic plasticity, and allozyme heterozygosity in Daphnia. Evolution. 1994, 48: 1715-1722. 10.2307/2410259.View ArticleGoogle Scholar
- Weber SL, Scheiner SM: The genetics of phenotypic plasticity. IV. Chromosomal localization. J Evol Biol. 1992, 5: 109-120. 10.1046/j.1420-9101.1992.5010109.x.View ArticleGoogle Scholar
- Scheiner SM, Caplan RL, Lyman RF: The genetics of phenotypic plasticity. III. Genetic correlations and fluctuating asymmetries. J Evol Biol. 1991, 4: 51-68. 10.1046/j.1420-9101.1991.4010051.x.View ArticleGoogle Scholar
- SanCristobal-Gaudy M, Elsen JM, Bodin L, Chevalet C: Prediction of the response to a selection for canalisation of a continuous trait in animal breeding. Genet Sel Evol. 1998, 30: 423-451. 10.1186/1297-9686-30-5-423.PubMed CentralView ArticleGoogle Scholar
- Hoaglin DC, Welsh RE: The hat matrix in regression and ANOVA. Am Stat. 1978, 32: 17-22.Google Scholar
- Felleki M, Lee D, Lee Y, Gilmour AR, Ronnegard L: Estimation of breeding values for mean and dispersion, their variance and correlation using double hierarchical generalized linear models. Genet Res. 2012, 94: 307-317. 10.1017/S0016672312000766.View ArticleGoogle Scholar
- Mulder HA, Hill WG, Vereijken A, Veerkamp RF: Estimation of genetic variation in residual variance in female and male broiler chickens. Animal. 2009, 3: 1673-1680. 10.1017/S1751731109990668.View ArticlePubMedGoogle Scholar
- Windig JJ, Mulder HA, Bohthe-Wilhelmus DI, Veerkamp RF: Simultaneous estimation of genotype by environment interaction accounting for discrete and continuous environmental descriptors in Irish dairy cattle. J Dairy Sci. 2011, 94: 3137-3147. 10.3168/jds.2010-3725.View ArticlePubMedGoogle Scholar
- Gilmour AR, Gogel BJ, Cullis BR, Thompson R: ASReml User Guide Release 2.0. 2006, VSN International Ltd: Hemel HempsteadGoogle Scholar
- Lee Y, Nelder JA: Double hierarchical generalized linear models. J Roy Stat Soc C. 2006, 55: 139-167. 10.1111/j.1467-9876.2006.00538.x.View ArticleGoogle Scholar
- Calus MPL, Bijma P, Veerkamp RF: Effects of data structure on the estimation of covariance functions to describe genotype by environment interactions in a reaction norm model. Genet Sel Evol. 2004, 36: 489-507. 10.1186/1297-9686-36-5-489.PubMed CentralView ArticlePubMedGoogle Scholar
- Rönnegård L, Felleki M, Fikse WF, Mulder HA, Strandberg E: Variance component and breeding value estimation for genetic heterogeneity of residual variance in Swedish Holstein dairy cattle. J Dairy Sci. 2013, 96: 2627-2636. 10.3168/jds.2012-6198.View ArticlePubMedGoogle Scholar
- ICAR: International agreement of recording practices. 2012,http://www.icar.org/pages/recording_guidelines.htm,Google Scholar
- Hill WG: Heterogeneity of genetic and environmental variance of quantitative traits. J Ind Soc Agric Stat. 2004, 57: 49-63.Google Scholar
- Mulder HA: Methods to optimize livestock breeding programs with genotype by environment interaction and genetic heterogeneity of environmental variance. PhD Thesis. 2007, Wageningen University, Animal Breeding and Genetics GroupGoogle Scholar
- Lillehammer M, Odegard J, Meuwissen THE: Reducing the bias of estimates genotype by environment interactions in random regression sire models. Genet Sel Evol. 2009, 41: 30-10.1186/1297-9686-41-30.PubMed CentralView ArticlePubMedGoogle Scholar
- Windig JJ, Calus MPL, Veerkamp RF: Influence of herd environment on health and fertility and their relationship with milk production. J Dairy Sci. 2005, 88: 335-347. 10.3168/jds.S0022-0302(05)72693-X.View ArticlePubMedGoogle Scholar
- Su G, Madsen P, Lund MS, Sorensen D, Korsgaard IR, Jensen J: Bayesian analysis of the linear reaction norm model with unknown covariates. J Anim Sci. 2006, 84: 1651-1657. 10.2527/jas.2005-517.View ArticlePubMedGoogle Scholar
- Robertson A: The sampling variance of the genetic correlation coefficient. Biometrics. 1959, 15: 469-485. 10.2307/2527750.View ArticleGoogle Scholar
- Burnham KP, Anderson DR: Model selection and multi-model inference: a practical information-theoretic approach. 2002, New York: SpringerGoogle Scholar
- Eskridge KM, Johnson BE: Expected utility maximization and selection of stable plant cultivars. Theor Appl Genet. 1991, 81: 825-832.View ArticlePubMedGoogle Scholar
- Mulder HA, Bijma P, Hill WG: Selection for uniformity in livestock by exploiting genetic heterogeneity of residual variance. Genet Sel Evol. 2008, 40: 37-59.PubMed CentralPubMedGoogle Scholar
- Ordas B, Malvar RA, Hill WG: Genetic variation and quantitative trait loci associated with developmental stability and the environmental correlation between traits in maize. Genet Res. 2008, 90: 385-395. 10.1017/S0016672308009762.View ArticleGoogle Scholar
- Makumburage GB, Stapleton AE: Phenotype uniformity in combined-stress environments has a different genetic architecture than in single-stress treatments. Front Plant Sci. 2011, 2: 12-PubMed CentralView ArticlePubMedGoogle Scholar
- Kliebenstein DJ: The quantitative genetics of phenotypic error or uniformity. Front Genet. 2011, 2: 59-PubMed CentralView ArticlePubMedGoogle Scholar
- Lee Y, Nelder JA, Pawitan Y: Generalized linear models with random effects: unified analysis via H-likelihood. 2006, Chapman & Hall/CRCView ArticleGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.