Genetic parameters for canalisation analysis of litter size and litter weight traits at birth in mice

The aim of this research was to explore the genetic parameters associated with environmental variability for litter size (LS), litter weight (LW) and mean individual birth weight (IW) in mice before canalisation. The analyses were conducted on an experimental mice population designed to reduce environmental variability for LS. The analysed database included 1976 records for LW and IW and 4129 records for LS. The total number of individuals included in the analysed pedigree was 3997. Heritabilities estimated for the traits under an initial exploratory approach varied from 0.099 to 0.101 for LS, from 0.112 to 0.148 for LW and from 0.028 to 0.033 for IW. The means of the posterior distribution of the heritability under a Bayesian approach were the following: 0.10 (LS), 0.13 (LW) and 0.03 (IW). In general, the heritabilities estimated under the initial exploratory approach for the environmental variability of the analysed traits were low. Genetic correlations estimated between the trait and its variability reached values of -0.929 (LS), -0.815 (LW) and 0.969 (IW). The results presented here for the first time in mice may suggest a genetic basis for variability of the evaluated traits, thus opening the possibility to be implemented in selection schemes.


INTRODUCTION
The ability of an individual to maintain its performance in the presence of environmental changes is known as phenotypic stability. In contrast, the ability of an individual to adapt its performance in the presence of environmental changes is called phenotypic plasticity. Phenotypic stability and plasticity are therefore two opposing concepts that refer to the same underlying base, the lower or higher adaptability of individuals to changing environmental conditions. Since the commercial success of livestock companies often depends on the homogeneity of animal performance, the selection to reduce environmental variability around an optimum production level, known as canalisation, has been a major focus for research in animal production. Some genes controlling environmental variability of traits were found to be different from those controlling the trait [28] and experimental results have confirmed the existence of plasticity genes [4,9,21]. Moreover, the interest in canalisation has increased due to the recent publication of evidence identifying molecular mechanisms [14] affecting canalisation in Arabidobsis thaliana [20] and D. melanogaster [25].
The modelling of environmental variability 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 variability when the environment is modified [28]. SanCristobal-Gaudy et al. [26] have proposed a model to deal with genetics of variability together with an EM-REML algorithm to estimate all the parameters simultaneously. Sorensen and Waagepetersen [29] extended this method to obtain results under a Bayesian approach. Based on these methods, the analyses of genetic parameters for canalisation were carried out on different traits and species such as litter size in sheep [27], within-litter standard deviation of birth weight in pigs [1,2,[11][12][13], birth weight in rabbits [6] or adult weight in snails [23].
Although the studies above have increased the knowledge on canalisation, a better biological understanding of the genetics of environmental variability is still needed. Hence, well-designed selection experiments with laboratory mammals, as models of livestock mammal species, are a necessary contribution to the literature on this important topic [23].
Litter size in domestic animals is a target trait for canalisation analysis. Successful selection for prolificacy in species such as pigs or rabbits has induced an increase in litter size. As a correlated consequence, mortality rate during lactation has also increased, whilst the ability of the young to survive has decreased [13,15,22]. The reduction in phenotypic variability of birth weight may likely improve survival rate, which is important for both economic and welfare considerations.
The aim of the present study was to estimate the genetic parameters associated with environmental variability for litter size, litter birth weight and individual mean birth weight in mice. For our study, we used an experimental mice population designed to reduce environmental variability in litter size. As a secondary objective, we compared two different approaches to solve the model defined for genetics of stability.

Experimental population
The experimental population analysed here started from a pre-extant mice population originating from a balanced genetic contribution of three inbred mice lines: Balb/c, C57BL and CBA. The three-way crossed population was maintained in panmixia during 20 generations thus ensuring high levels of both genetic and phenotypic variability.
From this panmictic population a total of 43 males were randomly selected to be mated with Balb/c inbred females. From these matings, we obtained a total of 875 females that were further mated with Balb/c inbred males. Litter size (LS), litter birth weight (LW) and mean individual birth weight (IW) were recorded for two consecutive parities. Records for LS and LW were obtained during birth inspections carried out every 24 h. IW was computed as the ratio IW = LW / LS. Note how the mating plan was designed to enable two different evaluations. The first was an evaluation that allowed dealing with within-male litter size variability on a constant female genetic background. The other was a search for the evaluation of litter size variability on half-sib daughters, by mating female F 1 with inbred line males in order to avoid most paternal influences on litter size. The inbreeding coefficients of animals belonging to an inbred line were asymptotically equal to one and were therefore considered as genetically identical animals in the pedigree file. All the available genealogical information for the non-inbred line individuals (including 18 generations back) was used for further analysis. Table I describes the structure of the analysed database. A total of 1976 records were available for LW and IW. The available records for LS included those obtained in the panmictic original population, which totalled 4129. The total number of individuals included in the analysed pedigree was 3997.

The model for canalisation analysis
Genetic parameters were estimated using the model proposed by SanCristobal-Gaudy et al. [26]. where y i j is the jth performance of a particular animal in a particular (animal × environment) combination i, the vectors b and b* contain the effect associated with generation (18 levels) as the only fixed effect considered, and x i , z i and w i are known incidence matrices. 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 p and p* contain permanent environmental effects for the trait and the log-variance, respectively, and are also assumed to be independent, with p|σ 2 p ∼ N(0, I p σ 2 p ) and p*|σ 2 p * ∼ N(0, I p σ 2 p * ) where I p is the identity matrix of equal order to the number of females having litters and σ 2 p and σ 2 p * are the permanent environmental variances affecting each trait and its log-variance respectively.
Two different approaches were used to solve the model. Exploratory estimates were first obtained via a rough solving procedure based on the availability of conventional software without Bayesian features. Then, a Bayesian approach was finally used to obtain final estimates.

First approach
This first approach was based on REML methodology and has been successfully used in rabbits [6]. Following this approach, the p* effect was excluded because conceptually the permanent effect is basically the trait that is being analysed.
The analysis was performed in three steps: 1. Firstly, genetic parameters of the trait were estimated via a homoscedastic model using REML. The matrix notation of the sets of mixed model equations to be solved was y = Xb + Zu + Wp + e, with: where G = Aσ 2 u , P = I p σ 2 p , R = I e σ 2 e , and y is the vector of observations, X the incidence matrix of fixed effects, Z the incidence matrix of animal effect, W the incidence matrix of permanent environmental effect, b the vector of unknown parameters for fixed effect, u the vector of unknown parameters for direct animal genetic effect, p the vector of unknown parameters for permanent environmental effect, e the vector of residuals, I e the identity matrix of equal order to the number of records, I p the identity matrix of equal order to the number of females having litters, A the numerator relationship matrix, σ 2 u the direct genetic variance, σ 2 p the permanent environmental variance, and σ 2 e the error variance.
The next two steps were carried out in order to test a heteroscedastic model. 2. The second step consisted in solving the homoscedastic model, obtaining a new vector of residuals by ê = y − Xb − Zû − Wp, and defining a new vector z by transforming each value in the ê vector by z i = ln(ê 2 i ). These new variables were called VLS, VLW and VIW, as measures of the environmental variability for LS, LW and IW respectively. Their means and variances are shown in Table I. 3. Finally, the new variables were analysed to estimate genetic parameters to complete the heteroscedastic model following the same sets of mixed model equations described above, but this time z instead of y, u* instead of u, e* instead of e and leaving out Wp. The permanent environmental effect, p* was excluded as described above, because it is basically the conceptual trait which is analysed in this step. This last step was carried out under a univariate model, but also under a bivariate model together with the original trait to obtain genetic and phenotypic correlations between the mean and variability of each trait.
A ratio h 2 v = σ 2 u * /(σ 2 u * + σ 2 e * ) can be obtained, in order to compare with previous estimates from the literature.
All analyses were done using the VCE programme [19], and by our own programmes when needed.

Bayesian approach
Sorensen and Waagepetersen [29] proposed the use of a Bayesian approach for canalisation analysis to better manage the model defined by SanCristobal-Gaudy et al. [26]. This Bayesian approach was successfully used in pigs [29] and snails [23].
Under this Bayesian approach two models were fitted: -The homoscedastic model is the classical additive genetic model, which assumes homogeneity of environmental variation, with the distribution of the vector of data y: y|b,u,p,σ 2 e ∼ N(Xb + Zu + Wp, Iσ 2 e ), where σ 2 e is constant. Vectors p and u were assumed to be a priori independently and normal distributed, that is: p|σ 2 p ∼ N(0, I p σ 2 p ) and u|σ 2 p ∼ N(0, Aσ 2 u ) , where A is the known additive relationship matrix. The vector b and the variance components σ 2 p , σ 2 u and σ 2 e , were assigned independent unbound uniform prior distributions. This model was solved using a Gibbs sampling algorithm. -The heteroscedastic model [26] assumes that the environmental variance is heterogeneous and partly under genetic control. The model and the distributions assigned for p, p*, u and u* were described at the beginning of the section describing the model. Scaled inverted chi-squared (ν = 4 and S = 0.45) distributions were assigned for variance parameters σ 2 u , σ 2 u * and σ 2 p , σ 2 p * , and a uniform prior bounded between -1 and 1 was assigned for ρ. In the statistical analysis for this model we used the Bayesian approach described by Sorensen and Waagepetersen [29].
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 algorithm described by Sorensen and Waagepetersen [29]. The mixing reflected in the Monte Carlo error was always smaller than 0.01. The conclusions of the paper regarding the posterior means are not changed by taking into consideration the Monte Carlo error. Convergence was tested using the criterion given in Gelman and Rubin [8]. 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 chain 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 inferences, a model identical to the Bayesian heteroscedastic model was fitted, except that the scale parameter S of the scaled inverted chi-squared prior distributions was set equal to 0.1 instead of 0.45.
The DIC (deviance information criterion) by Spiegelhalter et al. [30] is a combined measure of model fit and complexity. It is composed of 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 chose the model with the fewer parameters. It was used to test the second model compared with the first one.
Note how in this procedure, a residual effect e* does not exist and, as a consequence, an estimation of the ratio h 2 v defined above is unfeasible. There are several estimations of heritability for the traits under this procedure because residual variance varies among levels of the b effects. 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: Var and

Fitted models
Regarding approaches and random effects, five different analyses were fitted to ascertain: (a) the importance of the genetic parameters affecting LS, LW and IW in mice; and (b) the differential characteristics of two different approaches for canalisation analysis. The models solved under the first approach were the following: The three analysed traits were considered as traits of the mother in all the cases.

First approach
The genetic parameters estimated for LS, LW and IW using Model 1 are given in Table II The number of records, means and variances of traits and the artificial variables useful to deal with environmental variability are shown in Table I. The number of records for LW, VLW, IW and VIW traits (1796) is much lower than for VLS and LS (4129) suggesting that all results would be more reliable for this last trait. Moreover, the fact that inspection of births was carried out every 24 h resulted in differences in the interval from the birth time to the weighing time, which could have induced a very small error when registering the litter size, due to cannibalism, but could have greatly influenced the LW and IW traits.
The estimated variance components using Model 2 are given in Table III. The results for environmental variability must be interpreted with caution since they could be affected by the way the estimation is carried out in three steps. Heritabilities related to the trait were consistent with those previously found under a trivariate model whilst those for the variability of the traits were low ranging from 0.004 (LS) to 0.014 (IW).
The results obtained using Model 3 are given in Table IV. Heritabilities for the traits were slightly higher than those obtained using Model 2. Heritabilities for the environmental variability of the traits were also higher, from 0.0136 for LW to 0.0144 for LS. Table II. Heritabilities, genetic (above diagonal) and phenotypic (below diagonal) correlations and permanent environmental rate and correlation between permanent environments, and repeatability (t) for litter size (LS), litter weight (LW) and mean Individual birth weight (IW), with standard errors of parameters in brackets (Model 1).

Bayesian approach
The mean of the posterior distribution of variance components and corresponding ratios estimated using Model 4 are given in Table V. The results were consistent with those obtained from Model 2 and 3 although for IW the 95% highest posterior density intervals for the additive genetic component were extremely large. The means of the posterior distribution of variance components estimated using Model 5 are given in Table VI. Empirical distribution of variance components are shown in Figure 1. The decreases of the corresponding DIC values with those obtained from Model 4 are also shown in Table VI. The heritability values at different levels of generation effect are shown in Figure 2. Using Model 5, the estimated additive genetic variance increased and permanent environmental variance decreased for LS when compared with Model 4. For the same trait, a lower DIC value was obtained as well as a less negative correlation between the genetic effects affecting the mean and those affecting the environmental variability of the trait when compared with those obtained from Models 2 and 3. Regarding the difference in DIC values, Model 5 does not fit better for LW than Model 4. Finally, Model 5 provided the highest estimations of the heritability for IW, decreasing the DIC value of the model. The permanent environmental variances estimated using Model 5 are considerably lower than those estimated using Model 4 especially for the traits with lower DIC values in Model 5 with respect to Model 4.

DISCUSSION
In this paper, we present the first genetic analysis for environmental variability of LS, LW and IW in mice using an experimental population of mice designed to test canalisation. Although the heritabilities estimated here using Model 1 for LW (0.10) do not exactly coincide with the one we previously reported of 0.21 [5], both the genetic (0.96) and the phenotypic (0.93) correlations (Tab. I) estimated between LS and LW were very similar to those previously provided (1.00 and 0.96, respectively [5]). This is consistent with other authors: a genetic correlation of 1.0 was reported for LS and LW in mice by Luxford et al. [17] and for LS and litter weight at weaning by Luxford  and Beilharz [16]. Hubby et al. [12], in pigs, estimated a genetic correlation between litter size at birth and litter size at weaning of 0.98. IW had a much lower heritability than LS and LW and the genetic correlation between the pairs LS-IW and LW-IW was low. Since IW is a ratio trait, the results involving IW must be considered carefully. Notice that correlations among permanent environments were from moderate to high for any pairs of traits (0.36 and 0.94, respectively, in absolute values). Since both, individual phenotypic variability and permanent environmental effect, are defined on the basis of repeated measurements, these strong correlations among permanent environments suggest that a similar behaviour can be expected when estimating components of environmental variability. The permanent environmental component is between three and four times more important than the additive genetic component in IW. Again, since both, individual environmental variability and permanent environmental effect, are defined on the basis of repeated measurements, if a genetic component of environmental variability is present, it seems that this trait can be better managed by the models analysing environmental variability. The heritabilities for environmental variability provided by an initial exploratory approach are the first estimates carried out in mice. Even though the estimated heritability of environmental variability are very low, they are consistent with the value of 0.012 obtained for individual birth weight in rabbits [6] using this same methodology. Furthermore, they cannot be considered as negligible. In fact they were not estimated on the normal scale of a trait and, as a consequence, they may not be interpreted as normal heritability coefficients because transformation applied on the residuals lead to non linearity of the model and, thus, properties of regression are not applicable here [6]. The highest estimated heritability for environmental variability using Model 2 was found for VIW, as expected according to the important permanent environmental estimates for this trait obtained using Model 1.
Model 3 was used to obtain genetic and phenotypic correlations between each trait and its variability, but also to extract further information provided by the correlation between each trait and its variability. The estimated heritabilities for environmental variability with this Model were all roughly 0.014 for all the traits. These values were consistent with the value of 0.012 obtained for individual birth weight in rabbits [6]. The heritability for environmental variability was previously analysed by other authors but working directly on the standard deviation of individual weights. Damgaard et al. [3] found the heritability for within-litter variation in piglet birth weight to be 0.08 ± 0.03 and for weight at three weeks of age as 0.06 ± 0.03; the heritability for weaning litter size variation in pigs reported by Huby et al. [12] was 0.02 ± 0.01.
Genetic correlations between the trait and its variability were high and negative for LS (-1.00) and LW (-0.97) but low and positive (0.26) for IW. Sorensen and Waagepetersen [29] postulated that skewness of residual distribution provided information about the genetic correlation between the traits and their variability, a point also argued by Ros et al. [23]. Figure 3 shows the distribution of the estimated residuals from the second step described in the methodology. The skewness resulted negative for LS and LW, but positive for IW, agreeing with the correlations estimated between the corresponding direct genetic effect and the genetic effect for environmental variability; LS (-1.00), LW (-0.97) and IW (0.26). Genetic correlations concerning environmental variability have been scarcely studied. Huby et al. [12] reported a genetic correlation in pigs between litter size at birth and environmental variability of weight at birth of 0.29, which strongly agrees with the results found here in mice for similar traits. These authors also found a genetic correlation between litter size at weaning and the environmental variabilities of litter size at birth (0.01), litter weight at birth (-0.59) and litter weight at weaning (-0.75) which seems to indicate that a high negative genetic correlation exists between litter size and litter weight within age. Damgaard et al. [3] found a genetic correlation of 0.71 between within-litter Standard deviations of weights in piglets at birth and after 3 weeks. Zhang et al. [31] reported the existence of the whole range of genetic correlations between traits and their variability. Sorensen and Waagepetersen [29] working on litter size in pigs found a strong negative genetic correlation of -0.62, Ros et al. [23] working on adult weight of snails reported a value of about 0.80 for the same parameter, while Rowe et al. [24] analysing 35-day body weights of broiler chickens, calculated the value of a similar parameter to be about -0.10. Gavrilets and Hastings [7] and Hill [10] postulated that high genetic correlations between traits and their variability suppose that, due to pleiotropic effects, most alleles of genes controlling the mean can also act on the variance. These high genetic correlations found in LS and LW suggest the impossibility of canalising without important changes in the mean of the traits. However, it seems possible to select the environmental variance without changing the mean, despite the high correlation between additive genetic values affecting the mean and environmental variation. Ros et al. [23] showed how additive genetic effect for environmental variability can be decomposed into two terms, one of them independent of the additive genetic effect for the mean of the trait. They considered a restricted index I = u* + ku, choosing k so that Cov(u, I) = 0. Then, selection based on I should leave the mean of the trait approximately unchanged.
Model 4, under a Bayesian procedure, was used to validate or reject previous homoscedastic results and provided, as expected, the same estimates for the traits as univariate models under the initial exploratory approach (Tab. V). It also provided a DIC value for the homoscedastic model thus enabling the heteroscedastic model to be checked and compared for best fit. When the heteroscedastic model, Model 5, was fitted, the agreement between the estimation procedures was only the genetic correlation between each trait and its variability. There was a close agreement between LS and LW traits, and the sign of the genetic correlation for IW was in this case positive. Ros et al. [23] concluded that maximum-likelihood inference for this model was highly intractable. However, the initial approach was shown to provide exploratory estimates of the final parameters and could be useful in certain scenarios. Moreover, software utilities using the Bayesian procedure are still under development to consider a large number of situations, for example to work with highly unbalanced data, or to fit multitrait models. Finally, several freely available software packages, for example DFREML [18] or VCE [19], can be used for the canalisation model under the initial exploratory procedure used here. In contrast, the Bayesian methodology still requires a degree of computer skills to use it.
In model 5 under the Bayesian procedure, a parameter, such as heritability, cannot be considered for the environmental variability because a residual component is not considered. Thus, comparisons with the initial exploratory approach are unfeasible. However, an additive genetic component for the environmental variability of each trait was estimated. In this model, the estimates for permanent environmental variance decreased considerably when compared to the homoscedastic model. This might be because the heteroscedastic model captures the genetic variance of the additive genes concerning environmental variability from the permanent environmental component.
A different estimation of the environmental variance was obtained for each generation, which consequently provided different heritability estimates for the traits. Figure 2 shows the evolution of heritabilities through generations as well as the previous estimates under the homoscedastic model. All the heritability values were higher than those from the homoscedastic model suggesting that additive genetic variance of the environmental variability of the LS increased when compared with the homoscedastic model. Furthermore, the differences in the DIC values suggest an improvement in the performance of the model. For LW, with a higher DIC value, the heteroscedastic model was surprisingly worse than the homoscedastic model that had a lower DIC value, even when parameter estimates were consistent. It must be noted that the increase in the DIC value corresponds simultaneously to a lower estimation of the heritability for this trait. Finally, the results for IW were unexpected because the heritability of the trait improved considerably when the heteroscedastic model was fitted. However, conclusions about LW and IW must be taken with caution since models involving environmental variability should have a high data set size but the number of records for these traits is lower than 2000.
Trends in heritability are assumed to depend on the different residual variance by generation, but they may also be originated from differences in the additive genetic variance. A visual inspection of Figure 2 leads to the rejection of the latter idea since the trends of the graph are obviously erratic. Although only the last three generations for LW and IW had recorded data, a curious effect was seen since all three lines seemed to be approximately parallel which should imply that heterogeneity of residual variance had the same behaviour for all the traits involved.
The results presented in this paper and those from the literature could suggest a genetic base for environmental variability of the evaluated traits. Genetic improvement of these traits by artificial selection should be possible since other traits with low heritabilities have been successfully selected. The running selection experiment will hopefully help to confirm the existent genetic component of variability and may then be used to obtain an estimate of realised heritability.