Genetic parameters for litter size in sheep: natural versus hormone-induced oestrus

The litter size in Suffolk and Texel-sheep was analysed using REML and Bayesian methods. Litters born after hormonal induced oestrus and after natural oestrus were treated as different traits in order to estimate the genetic correlation between the traits. Explanatory variables were the age of the ewe at lambing, period of lambing, a year*flock-effect, a permanent environmental effect associated with the ewe, and the additive genetic effect. The heritability estimates for litter size ranged from 0.06 to 0.13 using REML in bi-variate linear models. Transformation of the estimates to the underlying scale resulted in heritability estimates from 0.12 to 0.17. Posterior means of the heritability of litter size in the Bayesian approach with bi-variate threshold models varied from 0.05 to 0.18. REML estimates of the genetic correlations between the two types of litter size ranged from 0.57 to 0.64 in the Suffolk and from 0.75 to 0.81 in the Texel. The posterior means of the genetic correlation (Bayesian analysis) were 0.40 and 0.44 for the Suffolk and 0.56 and 0.75 for the Texel in the sire and animal model respectively. A bivariate threshold model seems appropriate for the genetic evaluation of prolificacy in the breeds concerned.


INTRODUCTION
Litter size (LS) is economically the most important trait in lamb meat production [20] but it also has an important indirect effect on the improvement of other traits. A higher litter size allows more selection pressure on other economically important traits [22]. Because the heritability of LS is usually low, a selection on phenotype will be quite ineffective in improving litter size. The use of estimated breeding values, using BLUP and including information from relatives will substantially accelerate genetic progress. For this reason, the focus of the breeding programme in Belgian meat sheep is on improving the number of lambs born and appropriate genetic parameters are needed for the breeding value estimation procedure. In Belgian pedigree flocks of Suffolk (S) and Texel (T), the practice of hormonal induction of oestrus, followed by natural mating or artificial insemination (AI), is relatively common. In the period under study (1994−2002), 10% (S) and 24% (T) of all litters were born after hormonal treatment. There is no indication that these proportions have changed since. The use of hormones (such as Pregnant Mare Serum Gonadotrophin, PMSG) in sheep is known to cause an additional variation in litter size [3,5]. The effect of hormone administration varied with the level of prolificacy of the breed, the seasonal state of the ewe (anoestrus vs. oestrus) and the dosage. Consequently, the question was raised how litter size after natural oestrus (LSN) and litter size after induced oestrus (LSI) could be combined in a genetic evaluation of natural prolificacy. The average flock size of pedigree sheep breeders in Belgium is small, and discarding the litters after hormonal oestrus induction would lead to a significant loss of information. Moreover, the practice of AI would have been out of the picture because the litters resulting from AI, which is usually done after induced oestrus, would not have been processed in the genetic evaluation for litter size.
Genetic parameters for LSI are scarce in the literature as compared to estimates for LSN. In the Lacaune ewe lambs, the heritability of LSI was 0.05 and 0.06 in 2 data sets and the genetic correlation between the two types of litter size was 0.39 [2]. Other references on the heritability of litter size after induced oestrus and on the genetic correlation with natural oestrus in livestock were not found. With a correlation of less than unity and a difference in variance, it might not be optimal to treat "type of oestrus" as a fixed effect of LS.
Litter size in sheep, defined as the total number of lambs born per lambing is expressed in discrete numbers (1, 2, 3, 4 and 5). In many studies, LS is analysed by a linear model and variance components are obtained by REMLmethods. Arguments are that (1) non-linear models have no big advantages in goodness of fit or predictive ability as compared to linear models and (2) more computing time is required in non-linear models, which might be prohibitive for routine calculations [11]. Different non-linear models for litter size have been proposed, mostly sire models [13,19] but also animal models [17]. The results indicate that non-linear models are able to explain a larger proportion of the variation and increase the accuracy of prediction as compared to linear models. Especially for traits with low heritability and low incidence of some categories, a non-linear model becomes more advantageous compared to a linear model [6,19]. However, the analytical representation of non-linear models becomes difficult. A Bayesian analysis of an ordered, categorical trait is feasible using Markov chain Monte Carlo-methods (MCMC). MCMC methods allow for inferences to be made about the joint or marginal distribution of fixed effects, random effects and variance components. The Gibbs sampler, as an implementation of MCMC, has demonstrated its usefulness for inference on genetic parameters for litter size in sheep [1,28] and on twinning and ovulation rate in cattle [26]. However, simulation studies [12,15] with threshold animal models have shown convergence problems with the Gibbs sampler. This problem is related to the extreme category problem (ECP), which occurs when all observations within a level of an effect are in the same category. It turned out that sire threshold models are more reliable than animal threshold models [12,15].
The aim of this study was to study the genetic parameters of, and between, natural and induced litter size (LSN and LSI) in Suffolk and Texel sheep using different models and different approaches.

Animals
Litter size records from Suffolk and Texel in Belgium were selected from the National Association of Meat Sheep Breeders database (NAMSB). This database holds pedigree information and performance data from 20 breeds of sheep. All information in the system was supplied by the breeders and included the identification of ewes and rams, the date of mating (only for a limited number of litters) and the type of oestrus (natural or hormonally induced). At birth, lambs were individually recorded and the identification number, the date of birth, a code for mortality and the dam identification number were registered. Litter size was defined as the "total number of lambs born per ewe at lambing". Barren ewes were not recorded in the database.
Prior to analysis, data checks were made to discover incomplete or aberrant data. Records with obvious typing errors or with missing information were removed (3% of the records across breeds). The mating date was not generally available for analysis and this made the removal of records of ewes, which returned on induced oestrus impossible. Also, no details were given on the hormonal treatment at oestrus induction. For instance, active components and dosage of hormones were unknown and could not be included in the analysis.
From the complete data sets, only a subset (for each breed) was used in the analysis. Ewes were retained in the data only when their sire had at least two daughters with records and when among his daughters, the records of both types of oestrus were represented. The latter constraint could be met by subsequent records from the same ewe or by records from different ewes. Also, a minimum of three records per year*flock (YF)-effect was imposed for each trait.
This editing resulted in datasets with 1955 and 17 199 ewes descending from 249 and 1767 Suffolk and Texel rams respectively. All known ancestors of the ewes were included in the pedigree, resulting in 4123 Suffolk and 29 877 Texel animals in the analyses.

Models
Preliminary least squares analyses were conducted [21] for the traits separately, to find significant explanatory variables. The records were attributed to 6 classes according to the age of the ewe at lambing (AL), viz. 1 to 5 and 6 years and older. Records were also classified in two periods of lambing (PL) only in the Texel breed. Period 1 was the 30-day period covering the lambing peak, whereas period 2 contained all other litters. The PL-effect is particularly important in breeds with extended mating periods [11]. PL and AL were included as fixed effects in the model. A fixed year*flock-effect (YF) was included in the preliminary least squares analysis, together with a random sire effect.
For variance components estimation with REML and for Gibbs sampling, PL and AL-effects were included as fixed whereas the YF-effect was modelled as the random effect. Handling the YF as a random effect avoided convergence problems of the Gibbs sampler in threshold models [12,15]. Other random factors were the additive genetic effect of the animal in animal models (AM) or the sire-effect in sire models (SM) and a permanent environmental (PE) effect of the ewe. Random effects were assumed independent from each other, but a covariance was fitted between the same effect in the two traits. The PE-effect of the ewe affects the performance in consecutive litters. The significance of the PE-effect was evaluated in the REML analysis by fitting four models to the data; model I included the animal and PE effect for both traits (LSN and LSI). In models II to IV the PE effect for one or both traits was not fitted. The models can be considered as nested models and were compared using the Akaike Information Criterion (AIC). The AIC was computed from the likelihood value by properly accounting for the number of independently adjusted parameters and the rank of the matrix of fixed effects [27]. REML-estimates in the linear animal models were obtained using VCE4.2.5 [10,18] by fitting bivariate AM and treating litter size as a normally distributed trait. Heritability estimates and the variance of the YF and PE-effect were calculated as ratios to the phenotypic variance. Heritability estimates, calculated on the observed scale can be used to approximate the heritability for the underlying variable by the following formula where h 2 U is the heritability of the underlying variable and h 2 Cat the heritability of the observed variable [9]. The number of categories m equals 3, namely litters of singles, twins and triples and this implies 2 thresholds. Each category ( j) has an incidence of p j and a mean value of i j , z j is the height of the standard normal distribution corresponding with the proportion p j . Also the variance proportion of the PE-effect was rescaled in the same way. Genetic correlations and PE correlations estimated in the observed scale are direct estimators of the correlations in the underlying variable [9].
Inference on the genetic parameters for litter size was also made by a Bayesian approach. LSN and LSI were analysed in bi-variate threshold AM and SM, assuming an underlying variable and two thresholds for each trait. Inference on the parameters of interest was based on the marginal posterior distributions. The posterior distributions were obtained by Markov chain Monte Carlo using Gibbs sampling (GS), implemented in MTGSAM [25]. The technique of data augmentation was used to improve convergence [23]. Single chains of 1 000 000 and 500 000 samples were generated for the Suffolk and Texel respectively. Prior distributions for the variances and covariances were inverted Wishart distributions with expected values of 0.10 for the heritability, 0.05 for the PE-effect and YF-effect and 0.50 for the correlations. The residual covariance was set to zero because the observations on LSN and LSI are mutually exclusive. Prior distributions for the thresholds were uniform distributions with starting values calculated from the data [25].
Sequences of variance proportions and correlations were computed from the variances and covariances in the original chains. These sequences were further analysed with GIBANAL [24]. The burn-in period was determined for each parameter separately, as the number of samples before crossing the mean value of the chain for the second time. The burn-in samples were removed, leaving a "stationary phase" of the chain for further processing. Serial correlations were computed between samples and this was repeated at different lags to determine how many samples should be omitted between any two samples (thinning). The more correlated samples are, the more thinning is needed, resulting in smaller final sets. The final set was formed when the serial correlation between the samples was at the most 0.10. Location parameters of the posterior distribution (mean, mode, median and standard deviation) were computed from this subset of independent GS-samples.

RESULTS
Some details on the structure of the datasets are given in Table I. On average, 6.1 and 12.6 litters per year*flock were counted in S and T respectively. Only 2% (S) and 6% (T) of the YF-effects showed no variation in litter size because all litters were singles.
On average, sires had 7.9 (S) and 9.7 (T) daughters in the dataset with an average of 20.6 and 22.5 litters per sire. These numbers corresponded with an average lifetime production of 2.56 (S) and 2.31 (T) litters per ewe. In the Suffolk, the ewes had 1.99 litters after natural oestrus and 0.57 after induced oestrus, whereas in Texel the values were 1.59 and 0.73, respectively. From all active ewes, 27% (S) and 28% (T) had records for both types of LS. Most ewes exclusively had records on LSN (60% in S and 54% in T) and the remaining 13% and 19% of the ewes only had records on LSI.
The average litter size after hormonal treatment was +0.04 (S) and +0.10 (T) lamb/litter higher than litter size after natural oestrus (Tab. I). Also, litter size variance was significantly increased after hormonal induction by 16% in the Suffolk and 27% in the Texel. Hormonal treatment of the ewes resulted in less litters with singles and more litters of three lambs being born, compared to natural breeding circumstances. In the Texel, also more twins were born after oestrus induction.
The preliminary analysis of variance of LSN and LSI resulted in models explaining between 19 and 31% of the total variance. The univariate model for LSI in the Suffolk breed showed the lowest R 2 and a non-significant effect of the sire (Tab. II).
The comparison of models in the REML analysis, using the AIC (Tab. III), led to a different best choice in the two breeds. In the Suffolk, a model without a PE-effect for LSI was most likely. In the Texel, model I, with an additive genetic effect and a permanent environmental effect for both traits, had the lowest AIC.
The heritability estimates for litter size on the observed scale, obtained with REML, ranged from 0.06 to 0.13 depending on the breed and the model used. Somewhat higher heritability estimates were obtained for LSN when the permanent environmental effect was not included in the model. This effect was not observed in LSI. In Suffolk, the heritabilities of LSI and LSN were comparable   REML-estimates of the genetic correlation between litter size traits ranged from 0.57 to 0.64 in the Suffolk and from 0.75 to 0.81 in the Texel. Standard errors in the Suffolk were about 4-times larger than in the Texel. The unconstrained estimate for the PE correlation was 1.26 (S) and 0.72 (T). The result in S can theoretically not be called REML.
The results obtained with Gibbs sampling are presented in Tables IV and V and marginal posterior distributions are depicted in Figures 1 and 2. The burnin periods varied from 5 to 928 samples, so most samples in the chain were kept for further processing. The final sample size for the different components varied between 321 and 20 000. After thinning, serial correlations between samples within each chain were at the most 0.10. The values in the final samples could therefore be considered sufficiently independent.
In the Suffolk, AM and SM yielded similar posterior distributions. Distributions were fairly symmetric with the mean, mode and median close to each other and they appeared close to normal. The posterior mean of the heritability of litter size was 0.10 and 0.08 for LSN and 0.11 and 0.10 for LSI in AM and SM respectively. The PE-effect accounted for about 0.06 to 0.09 of the variance and the YF-effect for 0.05 to 0.06. Table III. REML estimates for the heritability (h 2 ), fraction of the permanent environmental effect (PE 2 ) and genetic (r g ) and environmental correlation (r pe ) for litter size after natural oestrus (LSN) and litter size after induced oestrus (LSI) in the Suffolk and Texel. For the different linear models the Akaike Information Criterion (AIC) is given, standard errors are given between brackets. In the Texel, posterior distributions appeared close to normal. The differences were noted between AM and SM for the location parameters. In the AM, the posterior mean of the heritability of litter size was 0.18 for LSN and 0.11 for LSI. In the SM, corresponding values for the heritability were 0.07 and 0.05. The posterior means for the proportion of variance of the PE-effect were 0.06 and 0.16 for LSN and 0.05 and 0.10 for LSI in AM and SM respectively. For the YF-effects, the means of the posterior distribution ranged from 0.05 to 0.08 and there was no difference between AM and SM. Table IV. Marginal posterior location parameters (mean, mode, median and standard deviation) and final sample characteristics (burn in, sample size (N) and serial correlation (corr)) for the heritability (h 2 ), variance fraction of the permanent environmental effect (PE 2 ), variance fraction of the year*flock effect (YF 2 ) and genetic (r g ), permanent environmental correlation (r pe ), year*flock correlation (r yf ) for litter size after natural oestrus (LSN) and litter size after induced oestrus (LSI) in the Suffolk. The posterior means of the genetic correlation between LSN and LSI, were 0.44 (AM) and 0.40 (SM) in the Suffolk and 0.72 (AM) and 0.56 (SM) in the Texel (Figs. 1 and 2). The mean of the distribution of the correlation between PE-effects was 0.41 (AM) and 0.50 (SM) in the Suffolk and 0.44 (AM) Table V. Marginal posterior location parameters (mean, mode, median and standard deviation) and final sample characteristics (burn in, sample size (N) and serial correlation (corr)) for the heritability (h 2 ), variance fraction of the permanent environmental effect (PE 2 ), variance fraction of the year*flock effect (YF 2 ) and genetic (r g ), permanent environmental correlation (r pe ), year*flock correlation (r yf ) for litter size after natural oestrus (LSN) and litter size after induced oestrus (LSI) in the Texel.

Data issues
Both the quality of field data and the distribution of effects over animals cannot be controlled. Imbalance in the data and/or disconnectedness may cause biased estimates. In order to minimise the risk of bias, data were edited and "sparse cells" were removed from the dataset. From Figures 3 and 4, the cumulative % of sires appearing in 1, 2, 3, to n different YF-effects, can be evaluated. The percentage of sires appearing in only one YF-effect is 10.2% for LSN and 37% for LSI in the Suffolk. In the Texel breed, these percentages are respectively 12.4% and 30%. This indicates that most sires are represented in at least 2 YF-effects. In the Suffolk, one sire had observations on LSN in 44 different YF-effects. In the Texel, the highest count of YF-effects for one sire was 166. The distribution of sires over YF-effects could be interpreted as a minimal evaluation of the genetic connectedness. Also the genetic relationships between ewes and between ewes and rams will contribute to the genetic links.

Modelling issues
In the preliminary least squares analysis, about 20 to 30% of the variation in litter size could be attributed to different explanatory variables. The lowest coefficient of variation was obtained in the Suffolk and for the LSI trait. This set was the most limited in terms of the number of records and this might explain the poorer modelling. There is probably scope for improvement of the model by including other explanatory variables. The preceding litter size (at birth or at weaning), the weaning to mating interval or a variable describing the body condition of the ewe could possibly explain non-genetic variation. However, this information is currently not recorded by the NAMSB and an extension of the current data collection system in the flock book will be necessary. This extension is only advisable if the extra information can reduce the residual variance in the model which would result in a higher heritability and more genetic improvement.
The increase of the mean and the variance of litter size due to hormonal treatment is also noted in most of the French sheep breeds [3] and in Rasa Aragonesa sheep [1]. An adverse effect of hormonal treatment on LS has been found with a reduction of average LS only in prolific breeds, such as the Romanov, Bleu du Maine, the Flemish breed and in Flemish × Texelcrosses [3,5]. However, the effects of oestrus induction on LS found in our data might be slightly biased due to the following reasons: firstly, the ewes were not randomly assigned to oestrus induction. In Suffolk and Texel, the ewes treated with hormones before AI or planned matings were probably "preselected" on general appearance and this could have been a disturbing effect. Secondly, the ewes returning on induced oestrus cannot be distinguished (in our data) from ewes mated after normal, i.e. natural oestrus. The ewes returning on induced oestrus were included as "natural" but an impact of the hormonal treatment on the next fertilising cycle cannot be excluded. We think, however, that the disturbing effects are limited and that the data were sufficiently good for genetic parameter estimation.
The separation of the random PE-effect from the residual term is most critical for LSI because on average, only 1.45 (S) and 1.57 (T) records per ewe for LSI were available (i.e. 1124 records by 774 ewes in S and 12 550 records by 7972 ewes in T). For LSN, 2.28 (S) and 1.95 (T) litters per ewe were available (3885 records by 1701 ewes in S and 27 265 records by 14 009 ewes in T) which makes the separation of the PE-effect from the residual term feasible. With REML-techniques, almost no variance was attributed to the PE-effect, especially for LSI in S. Although S and T do not differ so much in the number of records per trait, per ewe, differences in the amount of genetic covariance between animals and in pedigree depth may also have affected the separation of variance in different components. In S, less daughters per sire were present than in T (on average 7.9 versus 9.7 ewes/sire) and also pedigree depth was lower in S as compared to T (3.5 versus 6 generations of known ancestors on average). The dataset properties mentioned above might explain why the "best" REML-model was different in S and T.
In the initial Markov chains, burn-in periods were of variable length but considerably shorter than in other studies where up to 5000 and 10 000 samples were discarded [1,28]. A flexible choice of the burn-in period is indicated to save samples. Initial serial correlations between samples were high, indicating that samples were strongly autocorrelated. The usage of data augmentation is known to increase auto-correlation in the chain [23] so that more wastage of samples occurs. In general, final sample sizes were smaller in the runs involving animal models compared to sire models. Another difference was seen for posterior distributions involving the PE-effect. In the sire models, the smallest samples were obtained for chains involving one or both PE-effects. This might indicate a higher autocorrelation of the Gibbs samples for the permanent environmental effects than for the genetic or the year*flock effect. In the animal models, chains involving genetic effects were more autocorrelated and resulted in smaller final samples for the heritabilities and genetic correlations.
Simulation studies with threshold sire and animal models [15] indicated convergence problems in the animal models and this was related to the extreme category problem. Convergence problems were not encountered in our analyses, possibly as the result of treating the YF-effect as random [12,15]. Also the proportion of YF-effects with a zero variance was restricted in our data (2% in S and 6% in T) so that the impact of extreme categories was limited.

Genetic parameters
A mean value of 0.07 for h 2 (LSN) in sheep can be computed from REMLestimates published before 1995 [8]. Recent REML-estimates for litter size in sheep ranged from 0.05 to 0.26 [4,7,11,14,16,17]. Our REML heritability estimates for LSN are in agreement with these studies. Heritability estimates for LSI in Lacaune ewe lambs were 0.05 and 0.06 in two datasets with 1328 and 2201 records, respectively. Only the last estimate was significantly different from zero [2] and this value corresponds well with our REML-estimates of the heritability of LSI in S and T.
In Suffolk, rescaled REML-estimates are higher than the posterior mean obtained with GS whereas in Texel there is good agreement between the rescaled REML-estimates and posterior means or modes obtained with GS.
Posterior heritabilities obtained with the Gibbs sampler in AM for Suffolk were close to the posterior mean of 0.08 reported in Rasa Aragonesa sheep [1]. The corresponding Texel results were above this value. GS-results in SM in this study were lower than the posterior heritabilities for litter size (0.26) in Sabi sheep [16] and in Baluchi sheep evaluated at three ages (0.34 to 0.43) [28]. In the latter study, no PE-effect was fitted which is different from our model.
Variance proportions of the PE-effect of litter size are usually small and in the range of 0.01 to 0.03 [4,17]. Also in our REML analyses, very little variance was attributed to the permanent environmental effect whereas with GS, the distribution of the PE-proportions had a mean between 0.05 and 0.16. Also, more variance was attributed to PE-effects in SM as compared to AM and this was seen in both breeds.
The genetic correlation estimates between LSN and LSI obtained with REML (0.57 to 0.81), were above the value of 0.39 reported in Lacaune ewe lambs using a sire model [2]. Genetic correlations obtained with SM and GS corresponded well with the results obtained by Bodin [2].
The differences in genetic correlation estimates found in this study were partly due to the use of an animal versus a sire model. But, breed differences might also exist. In prolific breeds, hormonal induction is known to depress average LS while increasing its variance and this might lead to differences in genetic correlation estimates between breeds. The moderate to high genetic correlation, but different from unity, implies that the two types of litter size should be considered as different traits. Possibly, different sets of genes are in play. Another reason for the genetic correlation being different from one might be that insufficient information is provided by the data.
The REML PE-correlation estimate has a large standard error. From the GS results, it is clear that the PE-correlation was difficult to assess accurately.
Final sample size of the PE-effect was smaller than for the YF or genetic effect, indicating a slower mixing of the chain. Slower mixing of the chain is caused by less information coming from the data. In the case of the covariance between PE-effects, information has to come from subsequent observations of the two traits in the same animal. Only 26% (S) and 28% (T) of the ewes have records for both types of oestrus, which means that the within-animal covariance is difficult to estimate and larger standard errors on the covariances and correlations can be expected.

IMPLICATIONS
The magnitude of the heritabilities (0.08 to 0.18 for LSN and 0.05 to 0.11 for LSI) indicate that genetic improvement of litter size is achievable. Usage of a threshold model will have an advantage over the use of a linear model because the heritabilities in the threshold models were higher than in the linear model. Accounting for the categorical nature of litter size is thus recommended.
A moderately high and positive genetic correlation between natural and induced litter size was found in both breeds and values ranged from 0.40 to 0.81. Selection for higher (natural) litter size can also be based on observations of LSI, but progress will be slower. Although natural oestrus is more "informative" about litter size than induced oestrus, the latter technique is needed to apply AI or planned matings. These reproductive techniques allow the spread of the best genes in the population and provide genetic links in the data. Not discarding the litters after oestrus induction, but including them as a correlated trait, makes the best use of the information. This can be achieved by implementing a bivariate threshold model for litter size in the Belgian Suffolk and Texel.