Estimation of Genetic Parameters for Litter Size and Growth Traits in a Prolific Line of Tunisian Barbarine Sheep

Simple Summary: The Barbarine breed is an authoctonous sheep breed that originated in Tunisia and was raised in low-input production systems. This study was conducted to estimate the genetic parameters of several ewe productivity traits in a prolific line of this Barbarine breed selected for increasing litter size. Estimates of heritabilities were low in general. Genetic correlations between litter size and growth traits were negative, indicating that genetic selection for litter size will result in negative genetic responses to growth traits. A selection index including both types of traits could be proposed to correct for this antagonism. The results from this study can be used as a basis for the genetic breeding program of this prolific line. Abstract: This study aimed to assess genetic parameters for ewe productivity in a Tunisian Barbarine sheep line. The traits studied were litter size (LS), birth weight (BW), weight at 90 days (W90), and average daily gain between 10 and 30 days (ADG13). A total of 3804 growth and 2726 lambing records were used. Bivariate linear and threshold animal models were fitted and analyzed using the Gibbs sampling methodology. Heritabilities for LS obtained with univariate threshold, bivariate linear, and threshold models were around 0.15, higher than the estimate obtained by a univariate linear model (0.09 ± 0.03). Direct heritability for growth traits remained consistent across models, except for W90 in the bivariate linear threshold model. Maternal heritability for growth traits was higher than direct heritability, ranging from 0.07 to 0.15, except for BW. The covariances between the direct and maternal effects of growth traits were slightly negative. Repeatability oscillated between 0.16 and 0.62. Direct genetic correlations between LS and the other traits were negative, varying from − 0.18 (LS-BW) to − 0.83 (LS-W90). Our results suggest that the threshold model may be the most appropriate for LS. A selection index including LS and growth traits may be proposed for routine genetic evaluation in this population.


Introduction
The Barbarine fat-tailed sheep is the largest breed in Tunisia (64% of the overall sheep population).It can be found in all environments, but mostly in low-input production systems.The breed is well adapted to prevailing agroecological environments, with ecotypes that enable farmers to access specific production niches [1].
Several breeding programs have been implemented to improve Tunisian sheep meat production, concerning traits as diverse as growth rate, meat quality, environmental adaptation, and prolificacy.In conventional Barbarine flocks, growth and carcass traits are the main breeding objectives [2,3].Most of the sheep breeds in Tunisia are lowly prolific, with an average prolificacy of about 1. 1 [4].Thus, a prolificacy-based selection program has been implemented since 1979 in the research center of Oueslatia, Kairouan, of the Tunisian National Institute for Agronomic Research (INRAT), creating a line of prolific ewes of the fat-tailed Barbarine sheep.This selection process resulted in a prolificacy of 1.6 [4].
The productivity of the Barbarine breed remains low due to the large heterogeneity of environmental conditions [5].To solve this problem, reliable estimations of genetic parameters for each trait are necessary to reach an optimum rate of genetic progress.However, recent parameter estimates of productivity are not available in this line.Moreover, genetic parameters have only been estimated in a few studies using the REML (restricted maximum likelihood) methodology [4][5][6][7].Bivariate animal models that combine litter size and growth traits, or threshold models, have never been used for this breed.
Litter size is a discrete trait and is not distributed normally.Wright [8] considered that a linear variable underlies categorical traits, and thresholds define which category is observed.However, several studies showed that the estimates of heritability or breeding values from linear and threshold models are highly correlated [9][10][11].Furthermore, some studies have reported that threshold models give biased estimates of variance components [12,13].
Then, the objective of this study was to estimate genetic parameters of productivity traits in a prolific line of Barbarine sheep using, for the first time, Gibbs sampling applied to univariate and bivariate linear and threshold animal models.

Data
Records of productivity and pedigree information were collected from a Barbarine prolific line raised at the INRAT research station in Ouesslatia, located in the center-east of Tunisia.A total of 3804 lamb weights and 2726 lambing records collected from 1989 to 2016 from 985 breeding ewes were used in this study (an average of 3.86 lambs born and 2.76 parities per ewe).Reproduction was managed within sire families, and one ram and its replacement were assigned to each family.For mating, females born in a given family were systematically moved to another family to avoid the increase in inbreeding.Replacement animals were selected based on litter size records.Mating was practiced in the spring, and consequently, the lambing season spanned from September to November.Lambs were weighed at birth and seven other times at intervals of approximately 21 days.Weaning was at approximately 3 months of age.The rearing system is based mainly on natural pasture; however, concentrate supplements were given during the breeding and late gestation periods.The total number of animals in the pedigree was 4443 (987 ewes and 285 rams).

Traits Analysed
Traits included in this study were litter size (LS) and weight and growth traits: birth weight (BW), weight at 90 days (W90), and average daily gain between 10 and 30 days (ADG13).The litter size was recorded on the day of lambing as the total number of lambs born.Birth weight was defined as the live weight of lamb measured at the latest 24 h after birth.Weights at 10, 30, and 90 days and average daily gain between 10 and 30 days were estimated by interpolation between weight dates, bracketing the standard age specified with the formula of Warwick and Cartwright [14].Values out of the range defined by the mean ± 2 standard deviations were excluded.The characteristics of the data used in the analyses are presented in Table 1.

Statistical Analysis
A preliminary analysis using the general linear model implemented in the 'GLM2' R package [15] was applied to test the significance of non-genetic effects considered to influence ewe productivity and to be included in the genetic model.The statistical models included the year of lambing, the month of lambing, the age of the dam, and a combination of sex and type of birth.
Univariate and bivariate repeatability animal models were fitted to estimate the variance components and genetic parameters.Each bivariate model included litter size and one of the remaining traits.Genetic correlations between direct and maternal genetic effects for the same trait were taken into account.
The following bivariate animal model was fitted for all traits: where A i is the effect of year of lambing (28 levels for LS and BW, 26 levels for W90, 27 levels for ADG13), M j is the effect of the month of lambing (3 levels, from September to November), ST k is the effect of combined sex-type of birth with four levels, Ag l is the effect of the age of the dam with six levels (Ag l , l = 2, 3, . .., 7 and more years), a m and m n are the effects of direct additive and maternal additive genetic effects, respectively, p n is the permanent environmental effect of the dam and e i is the residual of the model.The model for litter size had neither the sex type of birth nor the maternal additive genetic effects.For growth traits, the repeatability of ewe performance was computed as the total maternal effect (t m ) using the formula of Notter [16].

Bivariate Linear-Linear Model
All traits were modeled as continuous variables and assumed to be sampled from the following bivariate normal distribution: where b LS and b G were random vectors, including the effects of A, M, ST, and age; a LS and a G were vectors of individual additive genetic effects; m G was a vector of maternal additive genetic effect; p LS and p G were vectors of permanent environmental effects.X, Z, and W were known incidence matrices; R was the residual (co)variance matrix.

Bivariate Threshold-Linear Model
LS was modeled as a categorical trait, as described by Sorensen et al. [17].We assumed that LS was related to an unobservable continuous variable (liability), u LS being the liability vector.Thus, the joint distribution of LS and the other trait records was assumed to be conditionally normally distributed, as in the case of the bivariate linear-linear model, and y LS was replaced by u LS (LS reliabilities).The response of LS (y LS ) was modeled with the following distribution as in Casellas et al. [18]: where I(.) was an indicator function with the argument as defined within parentheses and t 1 , t 2 and t 3 are thresholds defining the three categories of response.

Univariate Linear and Threshold Models
Litter size was analyzed separately, with linear [19] and threshold [17] univariate models to compare with the results from the bivariate models.To compare linear and threshold models for litter size for goodness-of-fit, the MSE (mean squared error) statistic was computed according to Casellas et al. [18].Birth weight, weight at 90 days, and average daily gain between 10 and 30 days were analyzed under linear univariate models.
The marginal posterior distributions of all unknowns were estimated by Bayesian inference using Gibbs sampling [19].The TM software developed by Legarra et al. [20] was used for all Gibbs sampling procedures.After some exploratory analyses, chains of 1,000,000 samples were used, with a burn-in period of 200,000.One sample of each 100 was saved to avoid high correlations between consecutive samples.Convergence was tested using the Z criterion of Geweke [19].

Results
The descriptive statistics for litter size and growth traits for the prolific Barbarine line are presented in Table 1.Average litter size in the prolific Barbarine line was higher than the one observed in conventional flocks of the same breed (1.13 ± 0.34 in Djemali et al. [21]; 1.08 ± 0.23 in Lassoued et al. [22]).For birth weight and growth traits, values were lower in this prolific line than the ones reported in conventional flocks [5,21,23].
For all traits, coefficients of phenotype variation were high, varying from 19.9% for BW to 40.2% for ADG13.In the case of litter size, CV was estimated at 40% for this line and 25% for a conventional flock raised at the INRAT research station (result not published).
The estimates of variance components and direct-maternal genetic correlations obtained with univariate and bivariate models are summarized in Tables 2 and 3, respectively.The estimates of additive genetic and permanent environmental variances for LS in the univariate threshold analysis were higher (0.06 and 0.04, respectively, Table 2) than those in the linear analysis (0.03 and 0.025, respectively, Table 2), but the estimates of residual variances were similar (Table 2) in both models.In bivariate models, the additive genetic variance was similar to the value estimated in the univariate threshold model (Table 3).The permanent environmental variance increased (0.18, Table 3), and the residual variance was reduced (0.14, Table 3) in the bivariate threshold model.
The estimated parameters for BW were similar in all models (Tables 2 and 3).For W90 and ADG13, few differences were observed between models (Tables 2 and 3).In all models, the estimated residual variances showed an incremental increase according to the age of the animal.The maternal components were higher than the direct effects for W90 and ADG13.However, for BW, direct effects were more important than maternal ones.The estimates of genetic correlations between direct and maternal effects were high and negative (−0.57 for W90 in the bivariate linear model to −0.86 for BW in the univariate model), except for W90 in the univariate and bivariate threshold models (−0.45 and −0.32, respectively).
Features of the estimated marginal posterior distributions of heritability and repeatability for the traits studied with univariate and bivariate analysis are summarized in Tables 4 and 5, respectively.The mean and median were similar for all the traits, showing that, in all cases, the marginal posterior distributions were symmetric.
The heritability estimate of LS obtained with the univariate threshold model was 0.14, with a probability of 81% of being at least 0.10 and greater than the estimate obtained with the linear model (0.08; Table 4).Bivariate linear-linear and linear threshold models resulted in estimates of h 2 similar to those of the univariate threshold model, but with a greater k value (0.10) and the probabilities of h 2 being greater than 0.10 (95% and 97%, respectively, see Table 5).The estimates of repeatability were quite different between the models.The highest estimate of repeatability resulted when the threshold model was used in bivariate analysis, and its respective high-density interval included high values (Table 5).Estimates of direct and maternal heritability for growth traits were generally similar in all models, except for W90, which showed direct heritability close to zero in the bivariate threshold model in comparison with the other models (Tables 4 and 5).The direct heritability estimate for BW was 0.18, with a probability of 95% being at least 0.10 (k value; Table 4).Both W90 and ADG13 showed low direct heritability estimates (0.03 and 0.04, respectively, Table 4) relative to the estimate at birth.Estimates of maternal heritability for growth traits were greater than direct heritability, except for BW, but both remained relatively small.Values increased with age and varied from 0.06 for BW to 0.15 for W90 (Tables 4 and 5).The repeatability estimates of ewe performance were moderate for all traits in all models except for ADG13, which showed low repeatability in the univariate analysis (0.17, Table 4).Estimates ranged from 0.20 for W90 (Table 4) to 0.29 for BW (Table 5).
Features of the marginal posterior distributions of phenotypic correlations between traits are presented in Table 6.Estimates of phenotypic correlations were negative (P = 1, Table 6) in all models, and their high-density intervals contained only negative values.Estimates were slightly lower than the genetic correlations and ranged between −0.34 for LS-ADG13 and −0.60 for LS-BW.Table 6.Phenotypic correlations between the traits analyzed: litter size (LS), birth weight (BW), weight at 90 days (W90), and average daily gain between 10 and 30 days (ADG13) were obtained with bivariate analysis in the prolific Barbarine line.Features of the marginal posterior distributions of genetic correlations are shown in Table 7.All estimates of correlations between the LS and growth traits were moderately negative with a very low probability of being positive (P value; Table 7), and the highdensity intervals included both positive and negative values.Estimates ranged from −0.18 for LS-BW to −0.83 for LS-W90.

Trait
Table 7. Genetic correlations between the traits analyzed: litter size (LS), birth weight (BW), weight at 90 days (W90), and average daily gain between 10 and 30 days (ADG13) were obtained with bivariate analysis in the prolific Barbarine line.Finally, the standardized average trend in EBVs for the analyzed traits over time is shown in Figure 1.The estimates of slopes were: 0.0031 for LS, 0.0006 for BW, 0.0003 for W90, and 0.0004 for ADG13 in standard deviations of traits.A favorable positive increase in LS was found over time (0.053 lambs in 28 years).For BW, W90, and ADG13, the tendency did not present a clear pattern, and a slight increase was observed (0.0092 kg, 0.0325 kg, and 0.609 g for BW, W90, and ADG13, respectively, in 28 years).

Discussion
Phenotypic values showed that birth weight and growth traits are lower in the prolific Barbarine line, in contrast with the conventional one.These findings show that selection based on litter size produces an increase in mean prolificacy and a decrease in individual lamb weights.The high CV observed for LS can be explained by the existence of multiple births in the prolific line and could have a positive impact on annual genetic gain in litter size.Safari et al. [24] reported CV from 6 ± 1 to 19.8 ± 1.5 for growth traits and 34.1 ± 0.6 for litter size in different breeds (wool, dual-purpose, and meat).
Estimates of genetic parameters for LS in the prolific Barbarine line were reported in a few former studies [4,7].Estimates were 0.09, 0.13 for h 2 , and 0.22, 0.62 for r by Ben Gara [7] and Bedhiaf et al. [4], respectively.Genetic parameters in these studies were estimated using REML methodology with single trait linear models, and standard errors were not available.In the prolific Morrocan D'man sheep, obtained h 2 and r were 0.09 ± 0.04 and 0.19 through the REML procedure [25].
Several studies comparing linear and threshold models with real and simulated data have demonstrated the suitability of the threshold model for studying categorical traits [9,18,26].Nevertheless, other studies showed that the threshold model gave similar results as the linear one in terms of predictive ability and goodness-of-fit [11,27].
The estimated residual variances increased according to the age of the animal in all models.This pattern might be explained by the fact that, at a later age, the lamb is less dependent on its mother and more exposed to environmental conditions.Estimates of ge-

Discussion
Phenotypic values showed that birth weight and growth traits are lower in the prolific Barbarine line, in contrast with the conventional one.These findings show that selection based on litter size produces an increase in mean prolificacy and a decrease in individual lamb weights.The high CV observed for LS can be explained by the existence of multiple births in the prolific line and could have a positive impact on annual genetic gain in litter size.Safari et al. [24] reported CV from 6 ± 1 to 19.8 ± 1.5 for growth traits and 34.1 ± 0.6 for litter size in different breeds (wool, dual-purpose, and meat).
Estimates of genetic parameters for LS in the prolific Barbarine line were reported in a few former studies [4,7].Estimates were 0.09, 0.13 for h 2 , and 0.22, 0.62 for r by Ben Gara [7] and Bedhiaf et al. [4], respectively.Genetic parameters in these studies were estimated using REML methodology with single trait linear models, and standard errors were not available.In the prolific Morrocan D'man sheep, obtained h 2 and r were 0.09 ± 0.04 and 0.19 through the REML procedure [25].
Several studies comparing linear and threshold models with real and simulated data have demonstrated the suitability of the threshold model for studying categorical traits [9,18,26].Nevertheless, other studies showed that the threshold model gave similar results as the linear one in terms of predictive ability and goodness-of-fit [11,27].
The estimated residual variances increased according to the age of the animal in all models.This pattern might be explained by the fact that, at a later age, the lamb is less dependent on its mother and more exposed to environmental conditions.Estimates of genetic covariances between direct and maternal effects were all negative, indicating antagonism between these effects, and this should be considered when formulating breeding goals.
Previous studies on this Barbarine population carried out with the joint data of the prolific line and the conventional one reported σ am moderately positive for BW, W90, and close to zero for ADG13 [5,28].These values were estimated using the REML methodology and presented without standard errors.
A lot of variations in the values of the covariances between direct and maternal genetic effects on growth traits have been found in the literature.In a review by Safari et al. [24] on results in several sheep breeds (wool, dual purpose, and meat), most reported genetic correlations between direct and maternal effects were negative.Robinson [29] suggested that this negative correlation between direct and maternal genetic effects could be reduced by including a sire-by-year interaction in the genetic model.
In this study, direct heritability declined and the maternal one increased with increasing age, in contrast to the tendency reported by Safari et al. [24].The values reported by those authors were moderate and ranged from 0.15 ± 0.02 (BW) to 0.41 ± 0.02 (adult weight) for h 2 d and from 0.24 ± 0.03 to 0.04 ± 0.01 for the same traits for h 2 m .In the prolific Barbarine line, Ben Gara [7], using the REML procedure, obtained greater estimates of direct heritability for BW (0.25), W90 (0.24), and ADG13 (0.15), but those values were reported without SE and maternal heritability was not estimated.
Some studies reported greater direct heritability estimates for BW than maternal heritability in Brazilian Santa Ines sheep [30,31], similar to or higher in several (wool, dual-purpose, and meat) sheep breeds [24].Estimates of greater magnitude for maternal than direct heritability have also been reported in other studies in a joint analysis with prolific and conventional Barbarine lines and Santa Ines sheep [5,32].
Linear and threshold analyses showed similar results, except for the case of the genetic correlation between LS and W90 in the linear threshold model, which was different but remained negative (−0.83; Table 7).This difference could be explained by the lack of convergence detected by the Geweke test for this estimate.
Negative correlations estimated in this study indicated that selection based on increasing litter size would have a negative influence on individual weight traits.However, several studies in this prolific Barbarine line [7,33] and the prolific D'man sheep [34] concluded that litter size is highly and positively correlated with total ewe litter weights at various lambs' age values.Estimates of genetic correlations were 0.89 for LS at birth and total litter weight at birth and 0.70 for LS at weaning and litter weight at weaning [33].Consequently, despite the negative impact of the selection of this prolific line due to the negative correlations between litter size and individual weights, this could be compensated by the positive correlations of litter size to total litter weights.
To this date, there are no available estimates of phenotypic and genetic correlations between litter size and growth traits for the prolific and conventional Barbarine flocks.Estimates of genetic and phenotypic correlations between LS and growth traits in the literature are abundant and cover a wide range of values, including the ones obtained here ( [35] in Merino sheep; [36] in Danish populations of Nordic sheep breeds; [37] in Polypay sheep).The differences in estimates could be due to different causes: data structure, methodology of genetic evaluation, production system, farm management, etc.
In sheep, farm profitability depends on the occurrence of multiple litters and lamb growth.The weight of a lamb is influenced by its genetic potential and its mother's ability during the growth period.The selection of replacement females in the prolific line of Barbarine sheep is currently based on their litter size, ignoring their maternal abilities.The advantage of the prolific Barbarine line in terms of maternal ability motivates its inclusion in the breeding goal.This objective can be achieved by fitting an index for female traits, including litter size and one of the growth traits.

Conclusions
This study showed that the threshold model could be the most appropriate for estimating the genetic parameters of litter size.Higher values of the estimates of maternal heritability with respect to the direct ones in growth traits indicated that selecting for maternal ability could have a positive impact on improving weaning weight.However, this possibility could not be directly exploited because of the existence of antagonism between litter size and growth traits in this prolific line.A selection index including both traits should be designed.

Figure 1 .
Figure 1.Trend in average standardized EBVs over time for a litter of size (LS), birth weight (BW), weight at 90 days (W90), and average daily gain between 10 and 30 days (ADG13) in the prolific Barbarine line.

Figure 1 .
Figure 1.Trend in average standardized EBVs over time for a litter of size (LS), birth weight (BW), weight at 90 days (W90), and average daily gain between 10 and 30 days (ADG13) in the prolific Barbarine line.

Table 1 .
Descriptive statistics for litter size (LS), birth weight (BW), weight at 90 days (W90), and average daily gain between 10 and 30 days (ADG13) in the prolific Barbarine line.
n: number of records; sd: standard deviation; CV: coefficient of variation.

Table 2 .
Estimations of variance components (SD within parenthesis) and direct-maternal genetic correlations (Mcse within parenthesis) of the studied traits were obtained with univariate models in the prolific Barbarine line.

Table 3 .
Estimations of variance components (SD within parenthesis) and direct-maternal genetic correlations (Mcse within parenthesis) of the studied traits were obtained with bivariate linear and threshold models in the prolific Barbarine line.

Table 4 .
Heritability (h 2 ) and repeatability (r) of litter size (LS), birth weight (BW), weight at 90 days (W90), and average daily gain between 10 and 30 days (ADG13) were obtained with univariate analysis in the prolific Barbarine line.

Table 5 .
Heritability (h 2 ) and repeatability (r) of litter size (LS), birth weight (BW), weight at 90 days (W90), and average daily gain between 10 and 30 days (ADG13) were obtained with bivariate analysis in the prolific Barbarine line.