Longitudinal analysis of direct and indirect effects on average daily gain in rabbits using a structured antedependence model

Indirect genetic effects (IGE) are important components of various traits in several species. Although the intensity of social interactions between partners likely vary over time, very few genetic studies have investigated how IGE vary over time for traits under selection in livestock species. To overcome this issue, our aim was: (1) to analyze longitudinal records of average daily gain (ADG) in rabbits subjected to a 5-week period of feed restriction using a structured antedependence (SAD) model that includes IGE and (2) to evaluate, by simulation, the response to selection when IGE are present and genetic evaluation is based on a SAD model that includes IGE or not. The direct genetic variance for ADG (g/d) increased from week 1 to 3 [from 8.03 to 13.47 (g/d)2] and then decreased [6.20 (g/d)2 at week 5], while the indirect genetic variance decreased from week 1 to 4 [from 0.43 to 0.22 (g/d)2]. The correlation between the direct genetic effects of different weeks was moderate to high (ranging from 0.46 to 0.86) and tended to decrease with time interval between measurements. The same trend was observed for IGE for weeks 2 to 5 (correlations ranging from 0.62 to 0.91). Estimates of the correlation between IGE of week 1 and IGE of the other weeks did not follow the same pattern and correlations were lower. Estimates of correlations between direct and indirect effects were negative at all times. After seven generations of simulated selection, the increase in ADG from selection on EBV from a SAD model that included IGE was higher (~ 30%) than when those effects were omitted. Indirect genetic effects are larger just after mixing animals at weaning than later in the fattening period, probably because of the establishment of social hierarchy that is generally observed at that time. Accounting for IGE in the selection criterion maximizes genetic progress.


Background
The genes that an individual carries and the environment in which it lives determine its phenotype. In various situations, during the pre-weaning period or in the case of group housing for instance, other individuals are part of the environment of the focus animal. Thus, the phenotype of an individual is influenced by its interactions with other individuals of the group. Such interactions may be positive (e.g. maternal effects or sterile helpers in social insects [1]) or negative (e.g. competition for limited resources or aggressive social behavior [2][3][4]) and are driven by other phenotypes of the group mates that generally are not recorded or measurable but can have a genetic component. Consequently, genetics of group mates or the dam may influence the focal phenotype. Such effects are known as indirect genetic effects (IGE) [5,6] and have been reported to be important components of various traits in different species [2][3][4]7]. Their role in the evolutionary processes of wild species and in response to selection in livestock species have also been explored [8].
The intensity of social interactions can depend highly on external factors, such as the number of interacting individuals [9,10] and rearing conditions. For instance, Piles et al. [7] recently showed that IGE for average daily gain (ADG) were stronger for rabbits under a restricted feeding regime than for rabbits fed ad libitum. In addition, interactions between group members can also vary over time. Several studies have shown that aggressive behavior of animals in group housing conditions is generally stronger at mixing and tends to decrease over time [11][12][13], which suggests that IGE can vary over time. Nonetheless, very few genetic studies have investigated changes of IGE over time for traits under selection in livestock species. To overcome this, the objectives of our study were to analyze longitudinal records of ADG in rabbits using a structured antedependence (SAD) model [14] that includes IGE and to evaluate, by simulation, response to selection on longitudinal ADG using different selection strategies. The idea of the SAD approach is to model an observation at time t by regression on the preceding observations. We chose to model a longitudinal trait with IGE using a SAD model because of the following advantages over the two main classical approaches used in genetic studies for longitudinal data, i.e. character process (CP) and random regression (RR) models: (1) the SAD approach can account for non-stationarity at the level of variances and correlations [14], in contrast to the CP approach, while the RR approach cannot model a situation with stationarity at the variance level and correlation between time points different than 1; and (2) extension of the SAD approach to the multiple-trait case is straightforward, in contrast to the CP approach, and requires fewer parameters than the RR approach [15]. The same applies to the case of correlated random effects in the analysis of just one phenotype, which is the case of correlated direct and indirect effects in our study.

Animals and housing conditions
Animals were raised on the experimental farm of IRTA in Spain between July 2012 and June 2014. A detailed description of the experiment is in Piles et al. [7]. In short, after weaning (32 d of age), kits were housed in cages (0.38 m 2 ) of eight individuals and fed under a restricted feeding regime corresponding to 75% of ad libitum feed intake during 5 weeks. To obtain a feed restriction of 25%, the amount of food given during week j was computed as 0.75 times the average feed intake of contemporary kits that were on a full feeding regime during week j − 1, plus 10% to account for the estimated increase in feed intake of growing animals. For week 1, the amount of feed for the restricted feeding regime was computed from data recorded in previous experiments on the same line from animals on a full feeding regime raised during the same season (multiplied by 0.75, 10% increase not included). Feed (commercial pellets for rabbits with (weeks 1 to 4) or without (week 5) antibiotics, as detailed in Piles et al. [7]) was distributed once a day in a 3-place feeder. Water was available ad libitum (one nipple drinker per cage). Individual body weights (BW) were recorded at weaning and weekly after weaning. Average daily gain (g/d) for each week was calculated as the change in BW from the beginning to the end of the week divided by the number of days elapsed (7 ± 1 d). On the weighing day, information regarding the animal health status was also recorded. Groups with animals showing disease symptoms (not caused by antagonist behaviors) or groups who suffered death events during the week were discarded from the analyses. The final dataset comprised 11,255 ADG records from 3096 individuals born in 1106 litters. The pedigree included information on 7701 rabbits. Descriptive statistics for the weekly ADG records are in Table 1.

Data analysis
Let y ilm (w j ), be the ADG of animal i (1 ≤ i ≤ 3096), born in litter l (1 ≤ l ≤ 1106), raised in group m (1 ≤ m ≤ 387), during week (1 ≤ j ≤ 5). The linear mixed model used to study ADG was: where μ i (w j ) represents the fixed effects at week j, K i the set of the seven group mates of focal individual i, and DGE i (w j ), IGE k (w j ), l l (w j ), g m (w j ) and p i (w j ) are the direct genetic, indirect genetic, litter, group, and pseudopermanent animal effects for week w j , respectively. The litter, group and pseudo-permanent random effects were independent from each other and distributed as: l ∼ N (0, I l ⊗ l ), g ∼ N 0, I g ⊗ g , and p ∼ N 0, I p ⊗ p , where l , g and p are 5 × 5 covariance matrices corresponding to litter, group, and pseudo-permanent animal effects, respectively, for the 5-week period of observation, and I l , I g and I p are identity matrices of appropriate size. Conversely, the direct and indirect genetic effects were correlated (co)variance matrix for the genetic effects (5 weeks for the direct genetic effects and 5 weeks for the indirect genetic effects) and A is the additive relationship matrix based on pedigree. Possible non-null covariances between random effects at different times were taken into account using a SAD approach [16][17][18]. It should be noted that this model does not include a residual term to help convergence and to avoid identifiability problems between structured permanent and classical residual covariance matrices [19], as in previous studies using the SAD approach [15,17]. Thus, the residual variance was by definition included in the (co)variance matrix of the pseudo-permanent animal effects of the model. For a given random effect, independent from the other random effects of the model, p for instance, the general form of the SAD model of order α for animal i was: p i w j = α s=1 θ s,j p i w j−s + e p,i w j , where θ sj is the sth antedependence parameter for week j, and e p,i (w j ) is the error term, normally distributed with mean 0 and innovation variance σ p 2 (w j ). Parameters θ sj and σ p 2 (w j ) were assumed to be continuous functions of time: θ s,j = β s q=0 a sq w q j , for a function of degree β s , and The SAD models were then defined by the order of the antedependence (α), the degree of the polynomial for each antedependence parameter (β 1 to β α ), and the degree of the polynomial for the innovation variance (γ). We denote the different SAD models with those parameters as follows: SADα − β 1 . . . β α γ . For instance, a SAD1-12 model corresponds to a SAD model with antedependence of order 1, and degrees 1 and 2 for the polynomial functions of the antedependence parameter and the innovation variance, respectively. A detailed description of this model is in Additional file 1.
In the case of correlated random effects (direct and indirect genetic effects), the dependence between the two terms was also considered using a SAD model, which is a particular case of the multiple-trait SAD model that was proposed by David et al. [20]. The general form of the SAD model for correlated effects of order α and α′ can be written as (for two effects, DGE and IGE, and for j > max α, α ′ ) : θ sj DGE i w j−s + δ j IGE i w j + e DGE,i w j , where notations are the same as for the SAD model for independent effects and δ j is the cross-antedependence parameter at week j. The parameter δ j was assumed to be a polynomial function of time.
The fixed effects included in the model and described in Piles et al. [7] were first selected step by step by comparing nested models using the likelihood ratio test. The fixed effects were week * body size (10 levels) and week * batch (70 levels) combinations, litter size (7 levels), and parity level (4 levels). After selection of fixed effects, the order and degree of the antedependence parameters were selected for each random effect in a model that assumed independence between direct and indirect genetic effects. Selection was performed by comparing nested models using the likelihood ratio test. Then, the correlation between direct and indirect genetic effects was included in the model. The degree of the cross-antedependence function (δ j ) was also selected using the likelihood ratio test. The SAD model was fitted to the data using ASReml [21] and the OWN Fortran program that combines the single and multiple-SAD programs that we have recently developed (https ://zenod o.org/recor d/19203 6 [20]), such that the two approaches can be applied in the same model, making it possible to consider cross-antedependence between direct and indirect genetic effects. The OWN program is freely available online at https ://zenod o.org/recor d/12280 58.
Once the SAD model has converged, parameter estimates ( ω ) can be used to compute estimates of the (co) variance matrix of the different random effects. For a given random effect, ˆ = L ′D−1L −1 , where D is a diagonal matrix with innovation variance estimates as components, L is a lower triangular matrix with 1's on the diagonal and negatives of the cross-antedependence parameter estimates below the diagonal entries [22]. The total phenotypic variance was calculated for each week j and for a group size of 8 as [23]: where σ 2 DGEj , σ 2 IGEj and σ DGEj,IGEj are the direct and indirect genetic variances and their covariance for week j, respectively, σ 2 lj , σ 2 gj and σ 2 pj are the litter, group and pseudo-permanent variances for week j, respectively, and r is the average genetic relationship among cage mates, which was 0.16 in this study [7]. Estimates of direct, indirect and total heritabilities were computed for each week j as the ratio of the direct genetic, total indirect genetic ( 49σ 2 IGEj [1]) and total heritable variance IGEj [24]), respectively, relative to the total phenotypic variance. To obtain standard errors of the co-variance and heritability estimates, we performed a multivariate normal sampling approach of ω , as described in Houle and Meyer [25]: where H ω is the inverse of the information matrix at convergence. We sampled ω 10,000 times, computed ˜ and heritabilities for each sample, and removed samples that led to a non-positive semi-definite matrix of ˜ . Estimates of heritabilities and ˆ were the mean of the heritabilities and ˜ across samples and their standard errors were the standard deviation of the samples. In addition, we calculated correlations between elements of ˜ to assess the estimability of the covariance components in this complex model.

Simulation to assess response to selection
To evaluate the advantage of the SAD model with IGE for genetic selection of longitudinal traits that are affected by IGE in comparison with a model that ignores IGE, we performed a simulation study with genetic selection. The simulated population was a closed nucleus of discrete generations and constant size. The base population consisted of 30 unrelated sires and 120 unrelated dams. Each founder sire was mated to four founder females to give birth to 960 offspring (8 offspring per mating and a 50/50 sex ratio). Then, among the offspring, one male per sire and 120 females were randomly selected to be breeders of the next generation. During the genetic selection process over seven generations, the best male per sire (i.e. 30 males) and the best 120 females were selected at each generation to be breeders of the next generation. The choice of the "best" animals was based on their estimated breeding values (EBV) that were calculated according to the genetic evaluation strategies described below. Each selected male was randomly mated to four of the selected females, avoiding individuals from the same sire family. The population and phenotypes were generated to match the previously described experimental data design. Individuals in each generation were assigned to groups of eight animals, with each group originating from four fullsib families and each family contributing two progeny. Individual phenotypes (5 records per animal, corresponding to 5 weeks of fattening) were constructed according to a multiple trait model that considered ADG as a different trait for each week. Average daily gain of each animal was computed for each week as the sum of: a week effect, the animal's direct genetic effect for the corresponding week, the indirect genetic effects of the seven other pen mates for the corresponding week, a week-specific group effect, and a pseudo-permanent effect (correlated within animal): y im w j = µ i w j + DGE i w j + k∈K i IGE k w j + g m w j + p i w j , (= Eq. (1) without litter effects). For each generation, multivariate normal distributions with unstructured covariance matrices were used to simulate the random effect values for each time point (genetic effects, group effect and pseudo-permanent effect). To stick as close as possible to reality, (co)variance values used in the unstructured matrices were close to those obtained with the SAD model in the data analysis step (detailed values in Additional file 2). For the genetic effects, to better assess the impact of the correlation between direct and indirect effects on selection response, we considered three sets of parameters for the genetic (co)variance matrix to mimic weak, moderate, and strong genetic antagonism between direct and indirect effects. For these sets, the genetic covariance matrix was the same as that obtained in the data analysis part of the study but three different sets of values were assigned to the direct-indirect covariances: values similar to those estimated in the data analysis part of the study (strong antagonism), those values divided by 2 (moderate antagonism) or by 4 (weak antagonism). Direct and indirect breeding values of the founders were simulated as for the other random effects. Then, the direct and indirect breeding values of an offspring were simulated as the average breeding value of its parents plus a Mendelian sampling deviation, which was sampled from a multivariate normal distribution with covariance corresponding to half of the genetic covariance matrix.
Three selection scenarios were investigated. In scenario (1), a SAD model that ignored IGE was used to predict genetic effects [ y im w j = µ i w j + DGE i w j + g m w j + p i w j , SAD1-11 for all random effect functions] and the sum of the weekly direct EBV was used as the selection criterion. In scenarios (2) and (3), a SAD model (same orders and degrees for the polynomial functions as those selected in the data analysis part of the study) with IGE was used to predict direct and indirect genetic effects for each time point [Eq. (1) with litter effects excluded]. The estimate of the total breeding value for each time point (i.e. weekly TEBV) was then computed as the sum of the direct EBV plus 7 times the indirect EBV of the corresponding time point leading to five TEBV per animal. In scenario (2), the selection criterion was the TEBV of the first week, while the sum of the weekly TEBV was used as the selection criterion in scenario (3). Variance components were estimated with the SAD model that included IGE [selection strategies (2) and (3)] or not [selection strategy (1)] using the 4800 records of the 960 offspring of the base population. Then, variance components were fixed to their estimates and considered as known to predict EBV in the different generations of the selection population.
The simulation was replicated 150 times. After all replicates were run, we assessed the capability of the SAD model with IGE to correctly account for the longitudinal structure of the data by assessing the standard deviation and the mean bias of each estimated variance component. Responses to the three selection strategies were compared based on the mean ADG in each generation.

Parameter estimates
Mean ADG was 26.5 g/d for week 1, increased by 50% from week 1 to 3, and then remained constant (Table 1). This increase was associated with a small decrease in the coefficient of variation for ADG (38% for week 1 and 26% for week 5). Weekly ADG ranged from 0.7 to 89.0 g/d, which is quite large.
After selection of models, we retained the SAD1-11 model for all random effect functions and a polynomial function of degree 1 for the cross antedependence parameter δ that models the covariance between the direct and indirect genetic effects. The likelihood ratio test comparing this full SAD model to the same SAD model in which IGE were excluded was 40.52 (8 degrees of freedom).
Estimates of variances of the pseudo-permanent, litter, and group effects for each week and their correlations between weeks that were obtained from the parameter estimates of the SAD model are presented in Figs. 1, 2 and 3, respectively. The pseudo-permanent effect variance increased with time from 23.1 to 52.8 (g/d) 2 . The ratio of the standard error to the pseudo-permanent effect variance was quite stable over time (~ 0.03). Estimates of correlations between the pseudo-permanent effects for different weeks were rather low and not significantly different from 0, except between consecutive weeks from week 2 to week 5, for which they were low and negative[− 0.08 (± 0.02) for weeks 2 to 3, − 0.13 (± 0.02) for weeks 3 to 4, and − 0.18 (± 0.03) for weeks 4 to 5]. The variance of litter effects decreased with time from 8.7 to 4.2 (g/d) 2 . Estimates of correlations between the litter effects for the different weeks were always positive and tended to decrease with the time interval between weeks. These correlations were higher between the four last weeks (ranging from 0.37 to 0.84) than between week 1 and the other weeks (ranging from 0.14 to 0.38). The estimate of the variance of the group effect increased with time from 5.7 to 16.4 (g/d) 2 . Estimates of correlations between the group effects of consecutive weeks tended to be negative (ranging from − 0.20 to − 0.26) and null otherwise.
Estimates of genetic variances and correlations are in Fig. 4. Estimates of genetic variances were not significantly different between consecutive weeks for both direct and indirect genetic effects. However, the observed general trend was an increase of the direct genetic variance from weeks 1 to 3 [from 8.03 to 13.47 (g/d) 2 ] followed by a decrease [6.20 (g/d) 2 at week 5], and  Estimates of correlations between the direct genetic effects for different weeks were moderate to high (ranging from 0.46 to 0.86) and tended to decrease with the time interval between measurements. The same trend was observed for the indirect genetic effects for weeks 2 to 5 (correlations ranging from 0.62 to 0.91). Estimates of correlations of indirect genetic effects of the first week with those of other weeks did not follow the same pattern and were weaker (ranging from 0.33 to 0.47). Estimates of genetic correlations between direct and indirect effects were negative within (ranging from − 0.57 to − 0.89) and between weeks (ranging from − 0.33 to − 0.84).
Heritability estimates are in Table 2. Estimates of direct heritability were moderate for weeks 1 to 4, ranging from 0.15 to 0.24, and its associated standard error was stable over time (0.02 to 0.03). For week 5, the estimate of direct heritability was lower [0.08 (± 0.02)]. The estimate of the heritability of indirect genetic effects ( 49σ 2 IGEj /total phenotypic variance in week j) was moderate (ranging from 0.16 to 0.44). Although estimates of indirect heritability were not significantly different between weeks, we observed that this heritability was 2.5 times higher for week 1 than for the other weeks. Estimates of total heritability were low to moderate (ranging from 0.05 to 0.19) and were associated with large standard errors (0.04 to 0.10). Although estimates of total heritability were not significantly different between weeks, we observed that it tended to decrease from week 1 to 4.

Estimation efficiency of variance components
Estimates of variance components converged for both SAD models (with or without IGE) for only a percentage of the simulation replicates (for 79, 77 and 48% of the replicates for the scenarios with weak, moderate, and strong simulated antagonism between direct and indirect genetic effects, respectively). Consequently, replicates for which one of the two SAD models did not converge were removed from the analysis. Table 3 includes the bias and variability over the replicates that converged for each variance component from the SAD model that included IGE. Bias and variability of variance components were not affected by size of the genetic antagonism between direct and indirect effects, except for estimates of the indirect genetic variance and indirect genetic correlation. Coefficients of variation ( std(ω) E(ω) * 100% ) were in the same range for the direct genetic and group variances (~ 19 and 17%, respectively) but were higher for the indirect genetic variances and increased from 39 to 58% as the genetic antagonism decreased. The direct genetic variance was correctly estimated (the mean relative bias, E ω − ω /ω , ranged from 1 to 2%), while the group variance tended to be slightly underestimated (mean relative bias ~ − 4%) and the indirect genetic variance was overestimated (mean relative bias ranging from 18 to 27%). All correlations were correctly estimated, with a

Fig. 4 Estimates of genetic variances for each week (on the diagonal, SE in brackets) and of correlations between weeks and between direct and indirect genetic effects (below the diagonal)
bias close to 0, except the indirect genetic correlation, which was slightly overestimated (bias ranging from − 0.09 to − 0.03). The standard deviations of the group and the direct genetic correlations were in the same range (~ 0.07), while the standard deviation for the direct-indirect genetic correlation was higher (~ 0.15) and the standard deviation for the indirect genetic effects was the highest (~ 0.26).

Response to selection
Changes in the mean ADG per generation are shown in Fig. 5 for the three selection strategies and the three sets of genetic antagonism between direct and indirect effects. For the three scenarios, the mean ADG increased linearly with generation when a SAD model including IGE was used to predict EBV. For the three sets of genetic antagonism between direct and indirect effects, the increase was highest when selection was performed using the sum of the weekly TEBV obtained with a SAD model that included IGE. Selection on the TEBV of week 1 obtained with a SAD model that includes IGE or on the direct EBV obtained with a SAD model that ignores IGE resulted in similar increases in ADG per generation, except when a strong genetic antagonism between direct and indirect effects was simulated. In that case, the model with IGE outperformed the model without IGE, for which there was no response to selection. For all selection strategies, response in ADG declined as the genetic antagonism between direct and indirect effects increased.
The effects of the different selection strategies and scenarios across the different weeks of growth were also explored by assessing the mean ADG by week in the 7th generation, relative to the base generation (Fig. 6). When a weak to moderate genetic antagonism between direct and indirect effects was simulated, the same pattern was observed: ADG increased for all weeks and for all selection strategies. As expected, response to selection at week 1 was higher with selection on the TEBV of week 1 obtained with a SAD model that included IGE. The pattern was different with a strong genetic antagonism between direct and indirect effects. On the one hand, with selection on the sum of the

Table 3 Average and range of bias and dispersion of the variance component estimates obtained with a SAD model with IGE applied to the simulated longitudinal data
The average and the range is given for each criteria (average over five components for the variance (5 weeks  weekly TEBV obtained with a SAD model with IGE, the mean ADG increased for weeks 1, 3, 4 and 5, while no changes were observed for week 2. On the other hand, with selection on the TEBV of week 1 obtained with a SAD model with IGE, there was a high increase in the mean ADG for week 1, a clear decrease for weeks 2 and 3, no changes for week 4, and a moderate increase for the last week of observation. When a SAD model without IGE was used to select animals, the pattern was the opposite (increases in weeks 2 and 3 and reductions in weeks 1 and 5).

Discussion
One of the characteristics of the analysis of longitudinal data is the within-subject correlation of the measurements (in addition to between-subject correlations due to shared genetic or environmental factors). A large number of parameters (100) are necessary to model the covariance structure of our data using an unstructured model with the same random effects as described here. Several flexible approaches exist to model the covariance structure of the data with a reduced number of parameters, such as random regression (RR), SAD [26], and character process (CP) models [16,27]. Among these, we chose the SAD approach for several reasons: (1) it relaxes the stationary correlations assumption made in CP models [28]; (2) it has been shown to better fit the data compared to the RR model in many situations [16][17][18]20]; (3) it is less sensitive to the drawbacks reported for the RR model, such as border effects [29] and inability to properly estimate correlations that decrease rapidly over time [16]; and (4) it generally requires fewer parameters than the CP model to model covariance structures. The SAD approach has been used to perform genetic studies of several longitudinal traits, such as repeated measurements of weight, feed intake, reproduction traits [17,18,20]. It has also been used to model residual (co)variances in models for quantitative trait loci (QTL) detection in the framework of functional mapping [30] of longitudinal traits in animals [31] and plants [32]. To our knowledge, this is the first time that a SAD approach has been used in quantitative genetic mixed models to model a trait with indirect genetic effects. It should be noted that we tried to apply a RR model to these data using the REML approach but the estimation procedure never converged. The mean weekly ADG reported here are consistent with those reported over the whole fattening period under a feed restriction of 80% by Drouilhet et al. [33] (40.4 g/d). Because of the limited amount of food, the ADG of restricted animals was lower than that of animals fed ad libitum from the same population of rabbits [7]. In addition, there was an effect of feeding restriction on the growth pattern, with the deceleration phase and the inflexion point of the growth curve being delayed. Because animals were raised in collective cages, we included indirect genetic effects in the model to account for the competition between individuals to access feed. In a previous analysis of the same dataset that did not consider the longitudinal characteristics of the data, Piles et al. [7] confirmed that indirect genetic effects play an important role in ADG under restricted feeding. The significance of IGE was also confirmed in our study when the longitudinal characteristics of the data were accounted for by the likelihood ratio test that compared the SAD model with and without IGE. The importance of indirect genetic effects on the ADG of group-housed animals has also been reported in other species such as pigs [34,35] and in other breeding conditions (e.g. conventional barren vs enriched pens [36]).
Estimating IGE can be challenging [37] because the statistical power for detecting IGE is determined by the group-population structure [38,39] and there is a risk of confounding IGE with environmental (i.e. group) effects, leading to non-identifiability of the (co)variance components of the social model [40]. Even if the structure of the design in our study (small group size, large number of groups and 4 * 2 full-sibs per group) meets the requirements to detect IGE and a group effect was included in the model as a random effect to help identifiability of parameters, identifiability of the (co)variance components is not guaranteed [40]. Exploring the information matrix I ω helps to detect identifiability problems for an unknown parameter vector ω [41]. The condition number (square root of the ratio of the first to the last eigenvalue) of the information matrix I ω of the SAD model used here was large (512), which may lead to the occurrence of identifiability problems. Nonetheless, the condition number of the submatrix of I ω that considered only the first regression coefficients of the SAD functions for all random terms (11 * 11 matrix) was equal to 35, which indicates that dependencies of parameter estimates were not between parameters of different random effects. In fact, close inspection of the correlation matrix among parameter estimates (Fig. 7, H ω = I ω −1 ) showed that correlation between regression coefficients within polynomial functions were high, while correlations between SAD parameters related to different random effects were low. For instance, the correlations between parameter estimates for the genetic effects and those related to the group effect were moderately high (ranging from − 0.57 to 0.56, and averaging 0.14), which indicates that the information provided by the data was sufficient to disentangle common environment effects from indirect genetic effects in the SAD model. To better understand the relationship between variance components of the random effects obtained with estimates of the SAD model, we also computed the correlation between estimates at the level of the (co)variance components using Fig. 7 Correlation matrix between parameter estimates. For group, litter and permanent random effects ("perm"); parameters are presented in the following order: the two regression coefficients of the antedependence parameters, then the two regression coefficients of the innovation variance. For the genetic effects ("gen"), the parameters are presented in the following order: the first six parameters correspond to the regression coefficients of the antedependence parameters for the direct (first and second parameters), indirect (third and fourth parameters) effects, and the regression coefficients of the cross-antedependence parameters (fifth and sixth parameters). The four last parameters correspond to the regression coefficient for the innovation variance of the direct (seventh and eighth parameters) and indirect (ninth and tenth parameters) genetic effects. Cell with a (dot) indicates an absolute correlation value higher than 0.6 the 10,000 random samples of ω in ω ∼ N ω, H ω . The corresponding 100 * 100 matrix is shown in Additional file 3. Correlations between (co)variance component estimates of different random effects were generally low. Those of the five indirect genetic variances with the five group variance estimates were less than 0.27 (in absolute value), while those of the 25 direct-indirect genetic covariances with the group (co)variance estimates were less than 0.18 (in absolute value). Such results favor the possibility to separate the variance estimates of the different random effects in the model.
The relative importance of the litter effect was moderate, in accordance with previous studies [42,43]. It decreased with time (from 18 to 5% of the total variance), probably as a consequence of the well-known decrease over time of the maternal influence on progeny performance [44,45]. The relative importance of the group effect (approximately 15% of the total variance) was in line with previous results in other species (18% [35]). Because residuals were assumed to be independent among the individuals of a given group, the variance of the group effect included in the model originates from both social interactions between individuals (correlated residual effects) and physical differences between pens (location in the building for instance) [46]. Accounting for the environmental dependence between group mates in this manner is possible if, and only if, the covariance among residuals of group mates is positive [23]. This condition was verified week by week for this population in a previous study (correlation between residuals of group mates over weeks were 0.11, 0.17, 0.19, 0.18 for weeks 1 to 4, respectively) [7]. The slight increase in the relative importance of the group effect over the study period (from 11 to 20% of the total variance) could be related with this previously reported slight increase in the weekly correlation between residuals of group mates.
We obtained moderate estimates of the direct heritability, which is consistent with previous estimates reported in the literature for the same species 0.22 [33] and 0.18 [47]. The decrease in the direct genetic variance, and thus in heritability, at week 5 is probably related to the change in feed composition that occurred that particular week.
The moderate heritability of indirect genetic effects is in line with those reported in previous studies for pigs (~ 0.13 [48] and ~ 0.39 [23]). It is interesting to note that the indirect genetic effects were higher for week 1, which is in accordance with studies on the behavior of animals raised in groups that report a higher level of aggressiveness at mixing, which reflects the establishment of social hierarchy that generally occurs within 3 days after mixing. In the same population as used here, Dalmau et al. [13] observed signs of antagonistic behavior such as biting, displacement, and animals jumping on top of each other over the whole growing period but especially during week 1. In mice and pigs, indirect genetic effects are known to induce antagonistic behavior at mixing [49,50]. Furthermore, we observed that the genetic correlation of the IGE in week 1 with the IGE in other weeks was lower than the correlation between IGE in weeks 2, 3, 4, and 5, which also indicates that IGE expressed just after mixing is a different trait from IGE expressed later in life.
The genetic correlation between direct and indirect genetic effects was negative regardless of the week considered, which indicates that individuals with a positive (direct) breeding value for their own growth tend to have a negative (indirect) breeding value for the growth of their cage mates and vice versa. This is consistent with the idea of competition between animals for a fixed amount of feed. Nonetheless, although we found no evidence of lack of identifiability for the parameters in our analysis, we suspect that the estimates of the genetic correlations between direct and indirect genetic effects are higher than their real value. An antagonism probably exists between direct and indirect genetic effects but not with such a strong negative correlation. The antagonism of these two effects has also been reported in previous studies on restricted feeding. Muir [51] described a moderate to high negative (− 0.56) genetic correlation between direct and indirect genetic effects on weight at 6 weeks of age of Japanese quail in groups of 16 birds. However, Bergsma et al. [23] found a neutral relationship between direct and indirect genetic effects for pigs in groups of 6 to 12 pen mates of the same sex.
Based on our results, we hypothesize that the growth of animals under a restricted feeding regime can be improved and that the delay of their growth curve can be shortened by reducing the negative impact of indirect effects, especially during the first week after mixing. One way to do this is to provide an enabling environment that decreases indirect effects by limiting the feed restriction during the first week after mixing (85 instead of 75% of ad libitum feed intake instead, for instance). Currently, this strategy is applied on several rabbits farms (Tudela F., personal communication). Another strategy is to perform genetic selection to improve the ADG, as investigated here by simulation.
For the simulations, unstructured (co)variance matrices were used to simulate the data and the SAD model was used to estimate variance components for different weeks. To obtain convergence of complex SAD models, such as those used in this study; the most appropriate method involves starting with a simple SAD model (low antedependence orders and low degrees for the polynomial functions) and to increase the degree and order step-bystep by using estimates of the preceding reduced model as starting values. For the simulations, we did not use this