Genetic parameters related to environmental variability of weight traits in a selection experiment for weight gain in mice; signs of correlated canalised response

Data from an experimental mice population selected from 18 generations to increase weight gain were used to estimate the genetic parameters associated with environmental variability. The analysis involved three traits: weight at 21 days, weight at 42 days and weight gain between 21 and 42 days. A dataset of 5273 records for males was studied. Data were analysed using Bayesian procedures by comparing the Deviance Information Criterion (DIC) value of two different models: one assuming homogeneous environmental variances and another assuming them as heterogeneous. The model assuming heterogeneity was better in all cases and also showed higher additive genetic variances and lower common environmental variances. The heterogeneity of residual variance was associated with systematic and additive genetic effects thus making reduction by selection possible. Genetic correlations between the additive genetic effects on mean and environmental variance of the traits analysed were always negative, ranging from -0.19 to -0.38. An increase in the heritability of the traits was found when considering the genetic determination of the environmental variability. A suggested correlated canalised response was found in terms of coefficient of variation but it could be insufficient to compensate for the scale effect associated with an increase of the mean.


INTRODUCTION
The body weight at a given age and the weight gain at a given period of time are important economic traits in animal production. Feed efficiency is indirectly evaluated [18,35] or selected [28] through daily gain. Changes due to selection have been widely reported in mice and estimation of realised heritability for post-weaning gain and realised genetic correlations for postweaning gain and body weight are available [22]. Likewise, several studies relate the weight gain to fertility and prolificacy [11,26]. Moreover, unequal growth of contemporary animals creates a competition that results in differential mortality [2,25].
The models used in animal breeding usually assume homogeneous residual variances. However, there is some evidence of heterogeneity in residual variance for growth in beef cattle [10], backfat thickness in swine [32], and milk yield in dairy cattle [21]. Hill [15] and Garrick and van Vleck [9] studied the consequences of ignoring heterogeneity in residual variance and found it results in a loss of expected selection response.
The modelling of heterogeneity is based on the hypothesis of the existence of a pool of genes controlling the mean of the performance and another pool of genes controlling the homogeneity of the performance when the environment is modified [31]. San Cristobal-Gaudy et al. [29] developed a model to deal with the genetics of variability together with a way of solving it using an algorithm. This model has been applied to estimate genetic parameters of variability in different species and traits: litter size in sheep [30], weight at birth in pigs [1,16,17] and in rabbits [8]. Recently, Sorensen and Waagepetersen [33] described a Bayesian implementation of this model that has been applied to analyse litter size in pigs [33], adult growth in snails [27], litter size in mice [14] and uterine capacity in rabbits [19].
These studies provide statistical evidence for the additive genetic control of environmental variation. The presence of genetic variation at the level of the residual variance suggests the possibility of modifying it by selection. More homogeneous production will allow an easier processing of animal products, with a consequent reduction in costs.
A better biological understanding of the genetics of variability is needed before carrying out any improvement program at the commercial population level. Some selection experiments involving livestock species have been designed and are being carried out [2] to reach this goal, but, in order to reduce the generation interval, selection experiments with laboratory mammals are necessary [14].

281
There is also increasing scientific evidence on the existence of a correlated genetic response for environmental variability in major production traits in mammals. However, the importance of this response on the mean of the traits under selection is still poorly understood. The aim of this paper was to estimate the genetic parameters for environmental variability on weight at 21 days (W21), weight at 42 days (W42) and weight gain between 21 and 42 days (WG) in a selection experiment conducted to improve the weight gain in mice. Even though genetic trends were not an objective of this work, exploratory signs of correlated canalised response were investigated, and the corresponding expected consequences of a combined selection with the objective of increasing mean values and reducing variance, are addressed.

Data
The population of mice used in this study came from a previous project carried out to compare the response of three different selection methods for WG: (A) the classic selection, choosing animals according to their performance and randomly mating selected individuals; (B) weighted selection unbalancing the offspring of each animal according to their genetic superiority; and (C) the minimum coancestry method, as in selection method (A) but designing mating according to the minimum coancestry criterion. The selection experiment was carried out during 18 generations with three replicates per selection method [24]. Within each line and replicate, 32 males were evaluated for weight gain between 21 and 42 days (WG) and those with the largest WG were selected. Eight males were individually selected among these 32 evaluated males. Each selected male was mated with two females and contributed an equal number of offspring (4 |) to the next generation. The females were neither evaluated nor selected. At the end of this process, the whole data set consisted of 5273 records for W21, W42 and WG in males, and 9152 individuals in the whole pedigree file.

Models
Sorensen and Waagepetersen [33] have proposed the use of a Bayesian approach for canalisation analysis to better manage the model defined by San Cristobal-Gaudy et al. [29]. This Bayesian approach has previously been used in a mice population closely related to that analysed here [14]. Under this Bayesian approach two models were fitted: -The homoscedastic model (Model HO) is the classical additive genetic model, which assumes homogeneity of environmental variation: where y i is the performance of animal i, b the vector of unknown parameters for the mixed method-replicate-generation systematic effect with 163 levels (18 generations, 3 selection methods and 3 replicates by method = 18 × 3 × 3, and one level for founder population), u the vector of unknown parameters for the direct animal genetic effect, c the vector of unknown parameters for litter effect with 2649 levels, x i , z i and w i the incidence vectors for fixed effects, animal effect and litter effect respectively and e i the residual. A maternal effect was not explicitly fitted in the model. Ignoring such an effect might increase the genetic variability of the direct genetic effect. However, a previous analysis on performances, fitting together both litter and maternal genetic effects, showed that both effects are confounded and cannot be separated. Thus, maternal influence cannot be considered as ignored in the model, but fitted to a large extent throughout the litter effect. Vectors c and u were assumed to be a priori independent and with a normal distribution, that is: c|σ 2 c ∼ N(0, I c σ 2 c ) and u|A, σ 2 u ∼ N(0, Aσ 2 u ) , where A is the known additive relationship matrix.
-The heteroscedastic model (Model HE [29]) assumes that the environmental variance is heterogeneous and partly under genetic control: where * indicates the parameters associated with environmental variance, b and b * are the vectors associated with the systematic effect, u and u * the vectors associated with the direct genetic effect and c and c * the vectors associated with the litter effect. Incidence vectors x i , z i and w i have been defined in the previous HO Model. It must be noted that c and c * are fitting the litter effect but, as previously mentioned, it is assumed that they are also fitting most of the maternal effect. The genetic effects u and u * are assumed to be Gaussian: where A is the additive genetic relationship matrix, σ 2 u is the additive genetic variance of the trait, and σ 2 u * is the additive genetic variance affecting environmental variance of the trait, ρ is the coefficient of genetic correlation and ⊗ denotes the Kronecker product. The vectors c and c * are also assumed to be independent, with c|σ 2 c ∼ N(0, I c σ 2 c ) and c * |σ 2 c * ∼ N(0, I c σ 2 c * ) where I c is the identity matrix of equal order to the number of females having litters and σ 2 c and σ 2 c * are the litter effect variances affecting, respectively, each trait and its variation. There are several estimations of heritability for the traits under this procedure because residual variance varies among levels of the b effects [14,19,27]. In this case, the phenotypic variance is the variance of the conditional distribution of y i given b and b * , and the heritability parameter h 2 is the usual ratio of additive to phenotypic variance. Under the heteroscedastic model, these parameters are the following: and It has to be pointed out that under Model HE different ratios for h 2 i are obtained for each combination of levels of the systematic effects ((Xb * ) i ). Details can be found in Sorensen and Waagepetersen [33] and Ros et al. [27].
The vectors b and b * were assigned bounded uniform prior distributions. Scaled inverted chi-squared (ν = 4 and S = 0.45) distributions were assigned for variance parameters σ 2 u , σ 2 u * and σ 2 c , σ 2 c * , and a uniform prior bounded between −1 and 1 was assigned for ρ.
The results for each model were computed by averaging the results obtained from two independent Markov chain Monte Carlo (MCMC) samples after running 1 000 000 iterations of the MCMC algorithms described by Sorensen and Waagepetersen [33]. Only one sample of each 50 was saved to avoid the high correlation between consecutive samples. The effective sample size was evaluated using the algorithm of Geyer [13] and Monte Carlo sampling errors were computed using time-series procedures described in Geyer [13], which were always smaller than 0.01. Taking into consideration the Monte Carlo error does not change the conclusions of the paper regarding the posterior means. Convergence was tested using the criterion given in Geweke [12]. For each variance, a scale parameter ("shrink" factor, √ R) was computed, which involves variance between and within chains. The shrink factor can be interpreted as the factor by which the scale of the marginal posterior distribution of each variable would be reduced if the chains were run to infinity. It should be close to 1 to convey convergence. The shrink factor was always between 0.99 and 1.15. In order to study the influence of the prior distribution on the posterior distributions, the models were analysed using different parameters for the inverted chi-squared prior distributions; the S parameter of the scaled inverted chi-squared prior distributions was set equal to 0.1 instead of 0.45. The use of proper priors for the variance components was deliberately chosen in order to avoid improper marginal posterior distributions.
The DIC (Deviance Information Criterion) by Spiegelhalter et al. [34], is a combined measure of model fit and complexity. It has two terms, the first term measures the goodness of fit and the second term introduces a penalty factor for the complexity of the model. Between two models with the same goodness of fit, the DIC chooses the model with the fewest parameters. This was used to test the second model compared with the first one.

RESULTS
Variance components estimated using Model HO are given in Table I for all the traits. Heritability values ranged from 0.09 for WG, to 0.21 for W42. Variance components for the environmental litter component were higher ranging from 0.38 for W42 to 0.52 for W21. The posterior means of variance components, genetic correlations and their highest posterior density at 95% for the three traits under Model HE are given in Table II. These correlations assume that there is a linear association between the additive genetic value affecting the mean and the additive genetic value affecting the environmental variance. Therefore, the boxplot for posterior MCMC realisations under Model HO of averaged squared standardised residuals against groups of additive genetic values ordered according to increasing size were drawn to ensure that they had   an approximate linear trend [27]. In both models HO and HE, the litter component was more important than the additive genetic component, and the highest value was found for trait W21. Under Model HE, genetic correlation between traits and environmental variance, was negative for the three traits (−0.19 to −0.38).
In Table III, Model HO is compared with Model HE, for percentage change of the main variance components, and for differences in the Deviance Information Criterion (DIC). The DIC favours Model HE for all the traits. Under Model HE, genetic additive variance increased for all traits, particularly for the trait with the lowest heritability (WG) which had a 200% increase in value compared to Model HO. These increases in genetic additive variance were accompanied by a decrease in the variance of the litter variance also for all the traits. This change was also more important for WG (−25%). Under Model HE, heritability was estimated for each level of the methodreplicate-generation effect, which was the only fixed effect in the model. Heritabilities estimated in each replicate were averaged within selection method and generation and further plotted in Figure 1 for the three traits analysed. Heritability estimated using Model HE, compared to heritability estimated under Model HO (which is illustrated in Fig. 1 as a dotted horizontal line), reached in general, higher values, particularly for WG. Additional information in Figure 1 shows the (linear) trends of the heritabilities estimated using Model HE; they all decreased with generation regardless of the selection method and trait. Note that since these different estimations of the ratio h 2 i are based on different levels of the fixed effect for the variability, these trends may be understood as non genetic.

DISCUSSION
Heritability estimated using Model HO (Tab. I) for W21 (0.15) was lower than that for W42 (0.21), which was in agreement with the results reported by Eisen and Prasetyo [5] and, in a seminal paper, by Falconer [6]. Estimated heritabilities for these two traits were also in close agreement with those reported by Fernández et al. [7] for litter weight using DFREML [23] on the same mice population as analysed here. Heritability estimated using Model HO for WG (0.09) was clearly lower than that for the other traits, but it was in close agreement with the one found for another trait such as litter weight in a similar population [14]. The litter component was much more important than the additive genetic component, ranging from 0.38 for W42 to 0.52 for W21, and substantially higher than that of 0.14 reported by Fernández et al. [7] and by Gutiérrez et al. [14]. According to Gutiérrez et al. [14], traits with a strong second random component, are expected to benefit from the use of models considering a decomposition of the environmental variability (Model HE). This was especially true here for WG, which was the trait with the lowest additive genetic component estimated under Model HO.
The results from Model HE (Tab. II) show an important increase in the additive genetic variance when compared to Model HO (28%, 46% and 200% of the original values, respectively for W21, W42 and WG). These increases were accompanied by a much less important decrease of the variance of the litter component (3%, 6% and 25% of the original values under Model HO, respectively for W21, W42 and WG). This might confirm that Model HE captures the genetic variance of the additive genes concerning phenotypic variability from the permanent environmental component [14]. Additionally, differences between variance components of a given trait substantially decreased when Model HO and Model HE were compared, especially for WG which is the trait with the lowest additive genetic component estimated under Model HO (Tab. III). Furthermore, Model HE had a better fit than Model HO for all the traits when using DIC value to compare between them (Tab. III).
Gutiérrez et al. [14] observed parallel lines in the evolution of the heritabilities estimated over three generations under panmixia in mice for litter size, litter weight and mean individual weight, thus showing that the residual variance equally increased or decreased for the three traits from one generation to the next. Some similar behaviour could be argued from Figure 1 for the traits analysed, but it is difficult to draw any conclusions from this. To carry out such an analysis, heritabilities for the three replicates were averaged within generation, selection method and trait. Then, we computed all the 9 × 9 correlations between the increases from one generation to the next one in the ratios h 2 i . These correlations ranged from a minimum of 0.08 to a maximum of 0.85, which were always positive. This seems to confirm that the changes in the residual variability tend to have the same explanation for all the traits. On the contrary, the observed trend for the ratio h 2 i estimated across generations (Fig. 1), had a negative slope regardless of the traits analysed. This was especially true for W42, which had a high genetic correlation with the selection criterion (WG), but also the highest heritability, and the highest correlation between mean and variance. It is important to remember that this ratio must not be interpreted as heritability in the classical way, and it is only the part of the additive genetic variability in the total variability, which cannot be expressed without environmental references. Moreover, each h 2 i assumes different residual variances depending on the estimated level of the fixed effect b * , but it assumes the same additive genetic variance, which is in fact that estimated for the founder population. Thus, this is a non genetic trend, and there is no easy explanation for these trends. Apart from drift or response variability, other possible unknown causes could be influencing their trend.
The negative correlation found between the estimated posterior means of additive values affecting mean and variance (Tab. II) were consistent with the results reported by Garreau et al. [8] for body weight at birth and its variability in rabbits. However, Ibáñez-Escriche et al. [20] found no correlation for slaughter weight at 175 days in pigs while Gutiérrez et al. [14] found extreme positive and negative correlations depending on the trait, and Damgaard et al. [4] and Huby et al. [17] found positive genetic correlations between mean and variability for weight in pigs. Moreover, Zhang et al. [36] found in the literature a wide range of values for correlations between mean and variability. The negative correlation between mean and environmental variance (−0.19) for WG, means that we should expect the environmental variance for this trait to decrease in an experiment conducted to increase the mean of the trait. Moreover, the genetic correlations estimated between WG and the other two traits analysed here, were high, 0.68 at W21 and 0.94 at W42 [24], and the genetic correlation between mean and environmental variance for these traits was higher than that for WG (−0.31 for W21 and −0.38 for W42). Given the negative genetic correlations between mean and environmental variances found here, it is expected that environmental variance will decrease throughout the generations as a consequence of a correlated decrease in the environmental variability. However, these changes in the environmental variance depend on the functional relationship between mean and variance [27]. In our model we postulate a linear, stochastic relationship between mean and log-variance, and an incorrect choice of functional relationship could give the wrong results for the genetic correlation between mean and log environmental variance [27].
In order to visually check the observed evolution of the variability across generations, we plotted the trend in phenotypic mean of WG, the phenotypic variance of WG computed directly from the data, and the coefficient of variation for WG (Fig. 2). Mean WG seems to have increased with generation as a consequence of the selection process. However, the phenotypic variance did not decrease as would be expected as a correlated response. On the contrary, the trends corresponding to the coefficient of variation seem to have been negative and showed that a correlated canalised response may have effectively been achieved. However, this correlated canalised response seems to be insufficient to affect the sign of the trends in phenotypic variances over generations, even though there could be several reasons acting on the phenotypic variability, for example the Bulmer effect [3] at the beginning of the experiment, the inbreeding, the drift or the response variability. Since the mean of the trait increased across generations and the coefficient of variation decreased, while the phenotypic variance did not decrease, a scale effect seems to have acted on the variance, somehow compensating for the correlated canalised response. However, these trends are only exploratory signs and should be confirmed by estimating genetic trends from BLUP values of animals.
The present analysis confirms the existence of additive genetic control of the environmental variability, which has been widely reported by other authors [1,4,8,14,16,17,20,27] for different species and traits. If the environmental variability is computed for the upper and the lower limits of the 95% highest posterior density interval of u * , a ratio of 4 to 6 is found between the corresponding environmental variability depending on the trait. Additive genetic component controlling environmental variability was estimated on the residual variance. Even though, canalised correlated response is likely to decrease the environmental variance, other factors without genetic control can induce an increase in this variance thus preventing a phenotypic response of a reduction of the variability of the selected trait. San Cristobal-Gaudy et al. [29] suggested the possibility of defining a selection index combining breeding values for the mean and environmental variance of a given trait in order to optimise a selection program. Our results suggest that before implementing such a selection index, further studies are required to understand the functional relationship between mean and variance and their influence on the expected correlated response. Thus, it would be desirable to explore other models [36] with different relationships between mean and variance. Another possibility would be to validate the models by comparing their expected response to selection with a selection experiment for variability.