Genetic components of litter size variability in sheep

Classical selection for increasing prolificacy in sheep leads to a concomitant increase in its variability, even though the objective of the breeder is to maximise the frequency of an intermediate litter size rather than the frequency of high litter sizes. For instance, in the Lacaune sheep breed raised in semi-intensive conditions, ewes lambing twins represent the economic optimum. Data for this breed, obtained from the national recording scheme, were analysed. Variance components were estimated in an infinitesimal model involving genes controlling the mean level as well as its environmental variability. Large heritability was found for the mean prolificacy, but a high potential for increasing the percentage of twins at lambing while reducing the environmental variability of prolificacy is also suspected. Quantification of the response to such a canalising selection was achieved.


INTRODUCTION
Selection for increasing prolificacy in sheep, although leading to a better average litter size in selected lines, also leads to an increase in prolificacy variability. This phenomenon is well known for qualitative traits, where mean and variance are linked. Extreme litters are encountered in prolific ewes (Romanov; Finnish) with five or even more lambs per lambing, which is obviously unacceptable for ewe and lamb viability. Breeders would like to have litter sizes of two exactly -and not on average -or as often as possible. In many situations twins are the most profitable (Benoit, personal communication).
Based on the example of the French Lacaune breed, the aim of this work was to evaluate if sheep can be selected for the objective: "concentrating prolificacy on 2". For that purpose, data consisting of litter size measurements on Lacaune sheep were analysed, using a direct adaptation to ordered categorical data of the quantitative genetic model described by SanCristobal-Gaudy et al. [22] relative to continuous traits. The hypothesis was stated that factors affect the underlying mean and/or the underlying environmental variability. These factors can be environmental, but also genetic. Variance components were estimated, giving the amount of genetic control on the mean and on the environmental variability, in a polygenic context. Prediction of the response to a selection for twins, based on the previous genetic parameter estimates, was derived using Monte Carlo simulation. Finally, this approach was compared with more traditional methods.

Threshold model for polytomous data -Likelihood approach
As Gianola and Foulley [10], Foulley and Gianola [8] or SanCristobal-Gaudy et al. [23] for example, we consider the threshold Wright model, based on an underlying Gaussian random variable. Thresholds transform this continuous variable into a multinomial variable with J ordered categories. Let us define I as cells indexed by i as combinations of levels of explanatory factors. Multinomial data are observed: (N i1 , . . . , N ij , . . . , N iJ ) ∼ M n i+ ; (Π i1 , . . . , Π ij , . . . , Π iJ ) (1) with N ij as the number of counts in cell i for the jth category, and Π ij the probability that an unobservable Gaussian random variable Y i ∼ N (µ i , σ 2 i ) lies between two thresholds τ j−1 and τ j (falls into the j th ordered category). Setting τ 0 = −∞ and τ J = +∞, the following is obtained: where n i+ is the observed number of counts in cell i for all J categories: n i+ = j n ij . The underlying means µ i and variances σ 2 i are linear combinations of parameters to estimate: (3) ln σ 2 i = p i δ, (4) where x i and p i are incidence vectors, β is a vector of location parameters, and δ is a vector of dispersion parameters.

Estimation and hypothesis testing
The estimation procedure can simply be maximum likelihood, implementing for example a Fisher-scoring algorithm, exactly as in [8]. Moreover, the test of H 0 : K δ = 0 vs. H 1 =H 0 , where K is a full-rank matrix, is achieved with the log-likelihood ratio λ = −2(L 1 − L 0 ), where L 0 (resp. L 1 ) is the log-likelihood of model M 0 (resp. M 1 ) corresponding to H 0 (resp. H 1 ). Asymptotically, the statistic λ follows a chi-square distribution under the null hypothesis H 0 , with degrees of freedom equal to the difference in the number of estimated parameters between models M 0 and M 1 .

Bayesian approach
Furthermore, the Bayesian quantitative genetic model developed by SanCristobal-Gaudy et al. [22] is based upon the underlying continuous variable Y as follows: (5) ln (6) where are location parameters, and γ = (δ , v ) are dispersion parameters. The parameters β and δ have flat priors, in order to mimic a mixed model structure, while u and v represent genetic values, with a joint normal prior distribution: where ⊗ denotes the Kronecker product, A is the relationship matrix between the animals present in the analysis, σ 2 u and σ 2 v are additive genetic variances relative to the location and log variance of the trait, respectively, and r is the correlation coefficient between genetic values u and v. Note that the continuous random variable Y is Gaussian conditional on (u, v). Using a now common incorrect terminology, the expressions "fixed effects" and "random effects" will sometimes be used in the following.
Here, focus is on the genetic aspect of the modelling of multinomial data, by the introduction of two (possibly) related groups of polygenes acting on the trait mean and log variance respectively.
Following SanCristobal-Gaudy et al. [22,23], a sire model is written with replacing (5) and (6). Vectors u and v are genetic values of sires, and data are collected on their progeny.

Model fitting
Let us denote N = (N ij ) (i=1,...I)(j=1,...J) as the observation, σ 2 = (σ 2 u , σ 2 v , r) the set of variance component parameters, and ζ = (τ , θ , γ ) the other parameters with τ = (τ j ) j=1,...J as the thresholds. The logarithm L of the joint posterior distribution reads: (10) where q denotes the number of elements in vector u (or v). Estimation of parameters ζ via the maximisation of L with respect to τ, θ, γ presents no theoretical difficulty when variance components are known. A Fisher-scoring algorithm leads to extended mixed-model equations (see Appendix).
When variance components have to be estimated, we chose to base the inference on the mode of the log marginal posterior distribution of variance components σ 2 :σ 2 = Argmax ln p(σ 2 |N), (11) by extension of the usual case (σ 2 v = 0) where the previous equation leads to REML estimates of variance components.
An EM-type algorithm was implemented as in [9,22], using an iterative algorithm where two systems are involved. The first system consists of BLUP-like mixed-model equations, where variance components are replaced by their current estimates. Solutions of these equations give current estimates of ζ. The second system updates the variance component estimates. When r is set to zero, equation (11) reduces to usual REML equations. However, numerical integration is required for multinomial data; details can be found in the Appendix.

Data
Data were collected from Lacaune ewe lambs born over 11 years as the result of inseminations made from 157 sires in 57 flocks. These flocks were a part of a selection scheme implemented in the Lacaune population since 1975 for Litter sizes greater than 5 were grouped into the 5th and last category. The percentages of litters with 1, 2, 3, 4 and 5 or more lambs were 41.1, 47.5, 9.8, 1.5 and 0.1 respectively. The overall prolificacy of these ewes at their first lambing was 1.72.

Homoscedastic models
A usual homoscedastic threshold model is fitted, including the fixed effects YEAR, HERD, SEASON, AGE in an additive way, and a random sire effect (u/2), symbolically written as: (13) on the underlying mean, where u ∼ N 157 (0, σ 2 u A) is the vector of sire genetic values and A is the relationship matrix. Interactions were not taken into account in the model because of non-(or bad) estimability or statistical non-significance. The significance tests for the explanatory factors on the underlying mean are shown in Table I.
The estimation procedure of Gianola and Foulley [10] gave an estimate of heritability equal toĥ 2 u = 0.39.

Heteroscedastic models
The previous additive model for the mean was used throughout the next analyses.
(i) First, factors that have a significant effect on the underlying trait environmental variability were sought. A likelihood ratio test was implemented. The reference model is the homoscedastic model with only fixed effects, including a sire fixed effect (model of the form (8)-(9), without u nor v): The current model for the significance test for, say, the YEAR factor, is for example: Table II gives the results of a forward selection procedure for the model on log variances. It shows that only the sire (considered as a fixed effect) has a significant effect.
(ii) Then a mixed sire model (8) [23]). The fixed effects and breeding value estimates are compared with those obtained with the mixed homoscedastic threshold model. They are close to each other, although the ranking is not exactly the same (not shown).
A plot of estimated breeding values (û,v) ( Fig. 1) allows to apprehend the joint ability of the 157 sires to produce high or low litter size on average and with a high or low variability.
In Table III, two sires with a mean prolificacy of the same order of magnitude are compared. The former has a high dispersion while the latter is canalised. The heteroscedastic model detects these differences and predicts slightly better the probabilities for the five categories. The total number of parameters is higher in the heteroscedastic than in the homoscedastic model, The high estimate of genetic variance (σ 2 v = 0.23) and of heritability (ĥ 2 u = 0.34) can be viewed as a great potential for the population to be canalised toward the phenotypic optimum of two (twins are economically the best), with a reduction of the environmental variability. The next section is a first attempt to quantify the expected response to such a selection, as was done for continuous traits [22].

Objective
One of the general objectives is the minimisation of discrepancies from an optimum of the descendence performances. The simple example of sheep breeders who wish to maximise the proportion of twins, first prompted this work. A single lamb and more than three lambs are economically undesirable. The optimum is then Π 0 = (0, 1, 0, . . . , 0). In the remainder of the text, the focus will be on this particular target. Obviously, generalisations are straightforward without any conceptual addition.

Selection schemes
Simulated selection schemes were run 1 000 times in order to have accurate empirical responses to canalising selection. A fixed number (n p ) of unrelated sires were mated to n unrelated dams each, producing n daughters per sire family. Each daughter had one record (litter size), and the set of n performances in a sire family was used to evaluate this sire. Different indices were compared and are detailed later. For the likelihood-based indices, animals were treated as if they were unrelated. True variance components were used (otherwise mentioned). After sire ranking, n s sires were selected and produce n p males for the next generation. The selection scheme was hence the same as in SanCristobal-Gaudy et al. [22], except that the phenotype was not directly Let us denote by i the sire, j the category, Π ij the probability that father i has daughters with a litter size equal to j for j in the {1, 2, 3, 4, 5} set, n ij the number of daughters of sire i that have a j litter size, I(n i ) the index of sire i with n i = (n i1 , . . . n i5 ), 5 j=1 n ij = n. Two phenotypic selection indices were considered: (16) the empirical estimate of Π i2 , where the index P stands for phenotypic and O denotes on the observed scale; if the discrete trait is treated as continuous, as in [22], the index is: where C stands for continuous (data are considered as such),n i and S 2 i are the empirical mean and variance, respectively, of n i and y 0 = 2.
Data were also generated with σ 2 v = 0.001 and likelihood calculations were performed with σ 2 v = 0.25 and vice versa, to apprehend the impact of using a wrong model on selection efficiency. Moreover, the model was slightly complicated by adding a fixed effect with two levels, say a HERD factor. Each sire i was given at generation t a proportion α it (resp. 1−α it ) of daughters in herd 1 (resp. 2), with α it drawn from a uniform distribution U(0, 1). The following parameterisation was adopted: the two levels had effects equal to a and −a, respectively. The particular value 2a = 1.5 was used in the simulations. It corresponds to a large effect encountered in the analysis of the Lacaune data.
At this point the following question arises: how can one introduce fixed effects in the index of selection when the relation between breeding values and phenotype (or index) is nonlinear? In the traditional linear case, let us denotê µ k +û i the estimated index of animal i in environment k. Evidently, the ranks of these indices do not depend on the environments. This is not the case in the threshold model since the ranks of do depend on environment k. In our particular case, the aim was to select sires giving the maximum of twins whatever the herd. The chosen index was since each sire has a probability of 1/2 of having a daughter in herd 1, by construction. More generally, each likelihood-based index I L * of equations (18), (19), (20), and (21) is replaced by The effect of the herd was not taken into account in the phenotypic indices PO and PC.

Results
The six selection indices are compared in terms of mean prolificacy (Fig. 2), phenotypic standard deviation (Fig. 3) with the corresponding genetic progress for v (Fig. 4), and percentage of twins (Fig. 5) during 20 generations of selection, and n = 100 daughters per sire. The shape of the u genetic progress is the same as the shape of the phenotypic mean in Figure 2 (not shown). Similarly, the percentage of quintuplets (not shown) behaves like the phenotypic standard deviation (Fig. 3). More importantly, the equivalence of indices corresponding to the same model, no matter the scale in which it is calculated (Observed or Underlying), is to be mentioned: LhomO behaves like LhomU, and LhetO like LhetU.
The phenotypic variance and the percentage of quintuplets are stabilised by the PO index, while the phenotypic mean tends very slowly towards the optimum. The PC index shows no progress in the mean prolificacy. This can be explained by the fact that the strong effect of the environment is not taken into account; this omission increases the residual variance and hence drastically decreases the heritability. The selection is consequently quite inefficient in moving the mean towards the target. The selection is nevertheless very efficient in decreasing the variance. In contrast the likelihood-based indices show a fast increase in the main criterion, that is the twin percentage and consequently the mean prolificacy. Because of the discrete nature of the data, the strong increase in the mean is accompanied by an increase in phenotypic variance. As soon as the population has reached the optimum on average, the phenotypic variance decreases provided that a heteroscedastic model is used (indices LhetO and LhetU). If not, the variance and the percentage of quintuplets are maintained at a high and constant level. Note that the PC index, also leading to a high genetic progress for v but with a lower mean than the LhetO and LhetU indices, shows a reduction in phenotypic variance.
Since data are discrete, the link between the mean and variance is so strong that the underlying genetic progress in v, which is indeed high for the LhetO and LhetU indices (one genetic standard deviation gain in 10 generations of selection), is not visible on the phenotypic scale until the mean stops increasing. It is however possible to slow down the genetic progress of u in order to privilege the genetic progress of v and its phenotypic expression. This can be achieved by putting different weights in the index, like: Figure 6, the particular values w 1 = 1 and w 2 = 50 were chosen. Compared to the PO index (Fig. 6), the mean evolves faster towards the optimum, while the variance decreases, showing that the weighted index LhetU has the highest performances whatever the point of view (mean or variance evolution).     When a parameter σ 2 v is set to 0.252 in the heteroscedastic model, while its true value is 0, then the selection based on the heteroscedastic indices LhetO or LhetU acts as if the genetic variance σ 2 v was already null, i.e. the indices LhetO or LhetU are quite equivalent to indices LhomO or LhomU in this case. For example, the mean prolificacy is only 3% lower with heteroscedastic than with homoscedastic models, while the phenotypic standard deviation is also 2% lower after three generations of selection. This means that the heteroscedastic approach does not slow down the efficiency of the selection if a higher genetic variance in v is wrongly put in the model.
The previous figures aimed at understanding the global long-term behaviour of some canalising selection indices. In practice, for the particular Lacaune breed, the short-term response to selection is given in Table IV in terms of mean prolificacy, phenotypic standard deviation, underlying genetic progress and percentages of single, twin, triplets, quadruplets and quintuplets or more. In this case, n = 30 progeny per sire is assumed.

DISCUSSION
The first aim of this work was the analysis of the genetic components of litter size in the Lacaune sheep breed. A liability model was chosen, as is often done for the analysis of polytomous data in animal genetics. A high = 0.34 on the underlying scale) was found for mean prolificacy. This value is greater than estimates generally found in the literature but it was observed before in this particular sheep population by Bodin et al. [1]. Although the structure of the data seems suitable for giving unbiased heritability estimates, according to Engel et al. [5] and Engel and Buist [6], some authors like Matos et al. [15] remark higher heritability estimates with a sire model than with an animal model for litter size. Other estimation procedures could have been chosen such as the quasi-score used by Jaffrezic et al. [12], or MCMC techniques. The only advantages of an EM approach are the certainty of convergence of the algorithm to a local minimum of the function to optimise, and the slight modification of the traditional REML equations. But the need for a MC step in the EM algorithm leads to heavy computations, which may tell in favour of full MCMC techniques.
The infinitesimal model proposed by SanCristobal-Gaudy et al. [22] for continuous traits was extended here to polytomous traits via a continuous underlying variable, allowing the modelling of the environmental variability as is usually done for the mean. The year, herd, season and age have no significant effects on the variability of litter size in the Lacaune population, but the sire factor has an important influence. The inclusion of the relationship matrix allows the interpretation of the sire variance σ 2 v of the log residual variances in the underlying scale as an additive genetic variance. The estimate of this parameter was found equal toσ 2 v = 0.23; it corresponds to a maximum value of the ratio of sire variances on the underlying scale equal to σ 2 Max /σ 2 min = exp(v Max − v min ) ≈ exp(6σ v ) ≈ 18, which is pretty high. At present, this value, however, has no comparison in the literature.
The second aim of this work was the prediction of the response to a selection for homogenising litter size around the target of two lambs per lambing. This problem is already complicated in standard situations, due to nonlinearity. An immediate extension of the work of Im and Gianola [11] shows that the parent-offspring regression is nonlinear for polytomous data with more than two categories. Some of the heritability estimates proposed by Magnussen and Kremer [13] cannot be extended to multiple-category data. Analytical expressions for the selection response of a binary trait given by Foulley [7] are unfortunately not feasible when a multiplicative model is set on the underlying environmental variance. The simulations performed in the previous section were imposed by these analytical complications.
Quantitatively, canalising selection is less efficient here than for continuous traits, due to the relationship between phenotypic mean and variance for discrete traits. The Lacaune situation is particularly difficult since one aspect of the objective is the increase of mean prolificacy, whose consequence (the increase of phenotypic variance) has an opposite action on the other aspect of the objective (reduction of the environmental variance). Despite a high genetic progress on the underlying environmental variance, only a small part of this is reproduced on the observed scale.
In fact, the model assumes a constant genetic variance in the mean value of the underlying variable Y and fixed threshold values that define a limit to the possible reduction in phenotypic variance, corresponding to the case in which Var(Y) = σ 2 u . At the limit, the expected proportions of litter sizes are equal to 0.12, 0.76, 0.11, 0.003 and 10 −5 , in increasing order. No reduction in the genetic variance was envisaged for this theoretical limit. More flexible models, derived from a physiological analysis (as in the work of Mariana et al. [14]), or involving the effects of QTLs or major genes on mean prolificacy, might probably be required to make such mid-and long-term predictions of the response to canalising selection more realistic.
Qualitatively, the analysed indices can be ranked on the basis of their related selection responses. In every case, the indices based on a heteroscedastic model (LhetO and LhetU) gave the best results for this criterion. A gain in the selection of categorical traits based on a threshold model over a linear model was already pointed out by Meuwissen et al. [17]. Moreover, the omission of an environmental factor with large effect, like the HERD in the simulations, has disastrous consequences on the selection, stressed by the nonlinearity between breeding values and index. Long-term figures were given in order to understand the global dynamics of certain canalising selections. So far, the selection objective had been the increase of twin proportion for the next generation.
In practice however, short-or mid-term figures are interesting for breeders. Then, generation-dependent weights in the selection indices can be envisaged, generalising the use of weights as in index (25): w 1,t (µ +û i /2 − y 0 ) 2 + w 2,t 3σ 2 u + exp(η +v i /2 + 3σ 2 v /8) (26) or j=1,J c j,tΠj,t for generation t, these weights should be chosen optimally to maximise a selection objective over T generations: To be fully comprehensive, the quantity Π j,t in equation 27 must be calculated over all the possible levels of environment k as in (23): where p k,t is the incidence of level k in the whole population. Economic studies will estimate weights c 0, j,t (Benoit, personal communication). One must note that the Lacaune population analysed in this paper has been selected for increasing the mean litter size. The observed high heterogeneity in sire variances may be due to the presence of polygenes controlling the residual variance (sensitivity to the environment), as was done in this paper. Heteroscedasticity may also be due to a major gene controlling the mean and segregating in the population, with the progeny of homozygote sires being less variable than heterozygotes. A canalising selection will favour homozygotes by reducing the variability, and pertaining polygenes will move the population mean to the optimum. The existence of such a major gene is currently being tested by Bodin et al. [3]. However, the genetics of reproduction traits is difficult (see for example Bodin et al. [2]), and no tool is currently available for fully understanding the genetic determinism of litter size variability.
The Fisher-scoring algorithm requires the information matrix, which can be obtained from the Hessian matrix and the fact that (equation (1)) Elements of the gradient of L are equal to: where Ω − denotes a generalised inverse of Ω, with and The results presented in [8] are a special case of these equations with ξ i = 1 and r = 0. We present hereafter the elements of the inverse of the Fisher information matrix: