<figure><img src="image_1.png" /></figure>

HAL
open science

# Impact of environmental disturbances on estimated genetic parameters and breeding values for growth traits in pigs

Vincent Le, Tom Rohmer, Ingrid David

► To cite this version:

Vincent Le, Tom Rohmer, Ingrid David. Impact of environmental disturbances on estimated genetic parameters and breeding values for growth traits in pigs. Animal, 2022, 16 (4), 9 p. 10.1016/j.animal.2022.100496. hal-03653106

HAL Id: hal-03653106

https://hal.inrae.fr/hal-03653106

Submitted on 10 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Impact of environmental disturbances on estimated genetic parameters and breeding values for growth traits in pigs

<figure><img src="image_3.png" /></figure>

Vincent Le<sup>a,b,*</sup>, Tom Rohmer<sup>a</sup>, Ingrid David<sup>a</sup>

## ARTICLE INFO

**Article history:**

Received 25 June 2021

Revised 16 February 2022

Accepted 21 February 2022

Available online 23 March 2022

**Keywords:**

Longitudinal performance

Resilience

Resistance

Robustness

Simulation

## ABSTRACT

Due to the diversification of farming systems and climate change, farm animals are exposed to environmental disturbances to which they respond differently depending on their robustness. Disturbances such as heat stress or sanitary challenges (not always recorded, especially when they are of short duration and low intensity) have a transitory impact on animals, resulting in changes in phenotypes of production (feed intake, BW, etc.). The aim of this study was to evaluate the impact of such unknown disturbances on the estimated genetic parameters and breeding values (BV) for production traits. A population of 6 120 individuals over five generations divided into eight batches of 10 pens was generated, each individual underwent an ~100-day test period. A longitudinal phenotype mimicking piglet weight during the fatening period was simulated for each individual in two situations: disturbed and non-disturbed. The disturbed phenotype was modified according to the robustness of the animal and the intensity and duration of the disturbance that the animal was subjected to. Various sets of simulations (1 000 replicates per set) were considered depending on the type of disturbance (at the level of the batch, pen, or individual), the genetic correlation (negative, neutral, or positive) between the two components of the robustness (resistance and resilience), the genetic correlation (negative, neutral, or positive) between growth and the components of robustness, and the heritability of the components of robustness (weak or moderate). An animal model was used to estimate the genetic parameters and BV for two production traits: the BW at 100 days of age (BW<sub>100</sub>) and average daily gain (ADG). The estimated heritability of the production traits was lower in the disturbed situation compared to the non-disturbed one (reduction of 0.08 and 0.05 points respectively for BW<sub>100</sub> and ADG). The correlations between estimated breeding values of the observed phenotypes (EBV) and BV for production traits in absence of disturbance were lower in the disturbed situation (reduction of 0.04 and 0.06 points for BW<sub>100</sub> and ADG respectively) while the partial correlation between EBV and BV for robustness was not significantly different from 0 in the two situations. These results suggest that selection in a well-controlled environment with random disturbances of low intensities does not allow to improve animal robustness while it is less effective for improving production traits than selection under no environmental disturbances.

© 2022 The Author(s). Published by Elsevier B.V. on behalf of The Animal Consortium. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

## Implications

The impact of low-intensity and random disturbances, encountered in nucleus farms, on estimated genetic parameters and breeding values for production traits is unknown. This study showed that, for production traits recorded at a given time, this impact is noticeable but is not sufficient to concomitantly select for animal robustness and production. To reach this last objective, it would be interesting to study the dynamics of the phenotype

over time to decipher production from robustness, which is made possible by the development of automatic devices in farms.

## Introduction

Conventional animal breeding conditions are becoming more diversified, in line with the reduction of pharmaceutical inputs, diversification of feed resources, and reduction of energy costs of livestock farms (Rauw and Gomez-Raya, 2015). All these measures are due to societal and farmer demands (aims to improve animal welfare and minimize costs). In addition, environmental conditions are also subject to greater variability due to climate change (heat

<footnote><sup>a</sup>GenPhySE, Université de Toulouse, INRAE, ENVT, F-31326 Castanet Tolosan, France</footnote>

<footnote><sup>b</sup>Alliance R&D, 35650 Le Rheu, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: giang-nam.le@inrae.fr (V. Le).</footnote>
waves) and new regulations regarding rearing conditions (abolition of individual cages) (Star et al., 2008). These changes generate a range of suboptimal production environments, to which animals respond differently depending on their robustness (Knap, 2005). This robustness can be described in terms of resistance and resilience (Nguyen-Ba et al., 2020), which are respectively defined as the ability of an animal to minimize the direct impact of the disturbance (de Goede et al., 2013) and its ability to quickly return to the state before the disturbance (Colditz et al., 2016). When the disturbances that the animals are subject to are known and recorded (temperature, disease, feed, farm management changes, etc.), it is possible to take them into account in the genetic evaluation model in order to obtain unbiased and accurate estimates of genetic parameters and breeding values (BV) for the trait of interest according to the environment. Conventional models that account for genotype by environment interaction used in such cases are the multiple-trait (Meyer, 2009) and the reaction norm models (Kolmodin et al., 2002). When the disturbances are unknown, a reaction norm model with unknown covariates has been proposed, although it can suffer from identifiability issues (i.e. no unique solution) (Shariati et al., 2009). The effect of unknown disturbances is generally considered in the evaluation model by the inclusion of a contemporary group effect (most often the herd-year-season effect in ruminants (Van Vleck, 1987) and the pen-within-batch effect in monogastric animals (Frey et al., 1997)). This model assumes that all animals in the same contemporary group are exposed to the same disturbances and that there is no genotype by environment interaction (i.e. response to change in the environment is the same for all genotypes, meaning that all animals have the same robustness) (Phocas et al., 2016). In nucleus farms where animals are raised in a highly controlled environment, disturbances are expected to be of low magnitude. It is however unclear what is the impact of these disturbances on the estimates of the genetic parameters and BV in this situation especially when the phenotypes of interest are measured at a given time point and thus not necessarily recorded concomitantly with the exposure of disturbances. The objective of this study was to evaluate the impact of random low-intensity disturbances that animals can face during a given period of time on genetic parameters' estimates and estimated breeding values of the observed phenotypes (EBV). Simulations of the BW of fattening pigs was used to meet this objective. Impact of disturbances on genetic evaluation for BW at 100 days of age (BW100) and for average daily gain (ADG) over the fattening period was evaluated.

## Material and methods

The aim of the simulations was to reproduce the effects of disturbances on pig's growth that may occur in a pig production system. They comprise the simulation of i) a simplified population of non-overlapping generations reared in different pens and batches ii) disturbances of different duration, intensity and starting date during the fattening period that occur at the individual, pen or batch level iii) the resistance and resilience (unobserved traits) of each individual iv) the observed trait which is the daily BW of each animal during the fattening period ignoring the effect of the disturbance (non-disturbed phenotype) and including the effect of the disturbance (disturbed phenotype). In the latter case, the dynamic of the phenotype over time is modified according to the disturbance's intensity and duration; as well as the animal's resistance and resilience.

### Population simulation

The simulated oversimplified pig population consisted of 5 non-overlapping generations without selection. The founder generation

comprised 12 sires and 150 dams that were randomly mated to give birth to 1 200 phenotyped offspring (8 offspring per dam, 4 males and 4 females). Among the progeny, 12 males and 150 females were sampled randomly to be the parents of the next generation. The same process was repeated from the first generation (G1) to the last generation (G5). The final population without the founders comprised N = 6000 individuals. To mimic a pig production system whereby animals are raised in small groups during the fattening period; animals of G1 to G5 were, within each generation, randomly distributed across batches and across pens within a batch. Eight batches and 10 pens per batch were considered, leading to 15 individuals in each pen for each generation.

### Disturbance simulation

Three different types of disturbances were simulated: batch disturbances (all animals in the same batch are subjected to the same disturbance), pen (all animals in the same pen are subjected to the same disturbance), and individual (the disturbance acts on a single animal). All batches were subject to a batch disturbance, while pen disturbances and individual disturbances occurred in a certain proportion of pens and animals. The pens and animals affected by these disturbances were randomly sampled. Two proportions of pen and individual disturbances were considered (20% and 40%). They were chosen based on input from experts in the field (Canario, personal communication). The intensity, starting point, duration of a given disturbance were randomly sampled according to the distributions provided in Table 1. The intensities were considered to be constant over time. A given batch could suffer from one batch disturbance only (same for pen and individual level: a given pen could suffer from one pen disturbance only...), but disturbances at different levels were not exclusive. Thus, an animal could be exposed throughout its life to one batch disturbance, one pen disturbance, and one individual disturbance that can occur at the same time.

### Resistance and resilience simulation

The resistance (resisᵢ) and resilience (resilᵢ) of each animal were considered on a [0,1] scale. The resistance corresponds to the ability of an animal to minimize the direct impact of the disturbance on growth performance. A resistance value of 0 corresponds to an absence of resistance, the direct impact of the disturbance will be maximal. Conversely, a resistance value of 1 indicates an animal insensible to disturbance, the perceived intensity will be null and the disturbance will have no impact on growth. The resilience corresponds to the ability to quickly return to the state before the disturbance. A resilience value of 0 indicates absence of resilience, the animal will remain in the state it was at the end of the disturbance (i.e. constant difference with the value of the phenotype it should have had if there had been no disturbance), while a resilience value of 1 indicates strong resilience (the growth of the animal will be maximal in order to quickly return to the BW it should have had if there had been no disturbance). Resistance and resilience were simulated, for each animal, i = 1, ..., N, according to the following model:

$$ 
\begin{aligned}
\operatorname{logit}(\mathrm{resis}_i) &= \mathbf{x}_{ri}\beta_r + u_{ri} + e_{ri} \\
\operatorname{logit}(\mathrm{resil}_i) &= \mathbf{x}_{Ri}\beta_R + u_{Ri} + e_{Ri},
\end{aligned}
\qquad (1) 
$$

where β_r and β_R are the vectors of fixed effects (chosen arbitrarily, see Table 1) affecting resistance and resilience (two levels for each), x_{ri} and x_{Ri} the corresponding incidence row vectors, u_{ri} and u_{Ri} are the resistance and resilience additive genetic effects (distribution given later in the text), and e_{ri} and e_{Ri} are the residuals for the resistance and resilience of animali, with joint distribu-
<table><caption>**Table 1**
Distribution of the random variables and the values of the parameters used in the simulations of the BW of pigs.</caption><thead><tr><th>Parameter/ variable</th><th>Description</th><th>Distribution/value</th></tr></thead><tbody><tr><td>enter_age<sub>i</sub></td><td>Age in day when animal <i>i</i> enters the test period</td><td>U({0, 1, ..., 7})</td></tr><tr><td>μ<sub>j</sub></td><td>Population mean for y<sub>wd</sub> at day <i>j</i></td><td>μ<sub>j</sub> = 15 + <i>j</i></td></tr><tr><td>β</td><td>Vector of fixed effects for y<sub>wd</sub> andy<sub>d</sub>, one fixed effect with three levels</td><td>β = [1, 2, 3]<sup>T</sup></td></tr><tr><td>x<sub>ij</sub></td><td>Vector linking the fixed effect to the observed phenotype for animal <i>i</i> at day <i>j</i></td><td>[1,0,0], [0,1,0] or [0,0,1] randomly sampled</td></tr><tr><td>β<sub>r</sub></td><td>Vector of fixed effects for resistance, one fixed effect with two levels</td><td>β<sub>r</sub> = [1, -1]<sup>T</sup></td></tr><tr><td>x<sub>ri</sub></td><td>Vector linking the fixed effect to the resistance for animal <i>i</i></td><td>[0,1] or [1,0] randomly sampled</td></tr><tr><td>β<sub>R</sub></td><td>Vector of fixed effects for resilience, one fixed effect with two levels</td><td>β<sub>R</sub> = [1, -1]<sup>T</sup></td></tr><tr><td>x<sub>Ri</sub></td><td>Vector linking the fixed effect to the resilience for animal <i>i</i></td><td>[0,1] or [1,0] randomly sampled</td></tr><tr><td>int</td><td>Intensity of the disturbance at the batch level<br>Intensity of the disturbance at the pen or individual level</td><td>Γ(0.3, 1)<br>U((0.1, 2))</td></tr><tr><td>start</td><td>Start time of a disturbance</td><td>U([2, 3, ..., 100])</td></tr><tr><td>dur</td><td>Duration of a disturbance</td><td>U([1, 2, ..., 25])</td></tr><tr><td>σ<sup>2</sup><sub>a0</sub></td><td>Variance of the genetic regression coefficient a<sub>0</sub></td><td>4</td></tr><tr><td>σ<sup>2</sup><sub>a1</sub></td><td>Variance of the genetic regression coefficient a<sub>1</sub></td><td>2.4</td></tr><tr><td>σ<sup>2</sup><sub>a2</sub></td><td>Variance of the genetic regression coefficient a<sub>2</sub></td><td>0.8</td></tr><tr><td>σ<sub>a0a1</sub></td><td>Covariance between a<sub>0</sub> and a<sub>1</sub></td><td>0.62</td></tr><tr><td>σ<sub>a0a2</sub></td><td>Covariance between a<sub>0</sub> and a<sub>2</sub></td><td>-0.36</td></tr><tr><td>σ<sub>a1a2</sub></td><td>Covariance between a<sub>1</sub> and a<sub>2</sub></td><td>-0.28</td></tr><tr><td>σ<sup>2</sup><sub>e,j</sub></td><td>Residual variance at day <i>j</i>,<br>σ<sup>2</sup><sub>e,j</sub> = 11φ<sub>0</sub>(j)<sup>2</sup> + 8φ<sub>1</sub>(j)<sup>2</sup> + 0.6φ<sub>2</sub>(j)<sup>2</sup></td><td></td></tr><tr><td>σ<sup>2</sup><sub>u<sub>i</sub></sub> + σ<sup>2</sup><sub>e<sub>i</sub></sub></td><td>Sum of genetic and residual variance of resistance</td><td>1</td></tr><tr><td>σ<sup>2</sup><sub>u<sub>R</sub></sub> + σ<sup>2</sup><sub>e<sub>R</sub></sub></td><td>Sum of genetic and residual variance of resilience</td><td>1</td></tr><tr><td>h<sup>2</sup>ADG</td><td>Heritability of the average daily gain</td><td>0.29</td></tr><tr><td>h<sup>2</sup>BW<sub>100</sub></td><td>Heritability of the BW at 100 days</td><td>0.26</td></tr></tbody></table>

U(): uniform distribution, Γ(): gamma distribution.

Resistance and resilience were simulated, for each animal, $i = 1, \dots, N$, according to equation [1]:

$$ \operatorname{logit}(\operatorname{resil}_i) = \mathbf{x}_{ri} \beta_r + u_{ri} + e_{ri} $$

$$ \operatorname{logit}(\operatorname{resil}_i) = \mathbf{x}_{Ri} \beta_R + u_{Ri} + e_{Ri} $$

The phenotype in the non-disturbed situation was simulated according to equation [2]:

$$ y_{wd,ij} = \mu_j + \text{enter\_age}_i + \mathbf{x}_{ij}\beta + \sum_{k=0}^{2} a_{ki} \varphi_k (\text{age}_{ij}) + e_{ij} $$

tion $[e_r \ e_R] N(0, [\sigma_{e_r}^2 \ 0 \ 0 \ \sigma_{e_R}^2] \otimes I_N)$, with $I_N$ the identity matrix of size $N$, $e_r = (e_{ri})_i$ and $e_R = (e_{Ri})_i$, and $\otimes$ the Kronecker product. Note that the resistance and the resilience in Equation [1] are considered to be constant over time for each animal.

## Simulation of longitudinal phenotypes

We considered that each animal *i* enters the fattening period at an animal-specific age that ends at 100 days of age, thus leading to a test period of length $n_i$ that varies from 93 to 100 days (the enter-age is unique to each animal regardless of animal group, see Table 1). During the fattening period, a longitudinal phenotype that mimics the daily BW was simulated for each individual of generation G1 to G5 in the following two situations: without disturbance vs disturbances possible, indexed with *wd* and *d* respectively. For individual *i* = 1, ..., *N* and day *j* ($j \in \{1, \dots, n_i\}$), the respective longitudinal phenotypes are denoted $\mathbf{y}_{wd,i} = (y_{wd,ij})_{j=1,...,n_i}$ and $\mathbf{a}_{d,i} = (y_{d,ij})_{j=1,...,n_i}$.

The phenotype in the non-disturbed situation was simulated by a random regression model (Schaeffer, 2004) using second-order orthogonal polynomials:

$$ y_{wd,ij} = \mu_j + \text{enter\_age}_i + \mathbf{x}_{ij}\beta + \sum_{k=0}^{2} a_{ki} \varphi_k (\text{age}_{ij}) + e_{ij}, \quad (2) $$

where $\mu_j = 15 + j$ represents the population mean at day *j* ($j \in \{1, \dots, 100\}$) that increased linearly with time. The entity enter_age_i is the enter-age of animal *i*. The vector $\beta = [\beta_1 \beta_2 \beta_3]^T$, with incidence row vector $\mathbf{x}_{ij}$, is the vector of fixed effects arbitrarily chosen to reflect the effect of one cross-classified factor with three levels (the daily social status of the animal for instance: aggressor, victim and none). The level of the factor was randomly sampled for each individual and each day *j* from a discrete uniform distribution $U([1, 2, 3])$ (i.e. if 2 is sampled then $\mathbf{x}_{ij} = [0, 1, 0]$) (see Table 1). The entity age_ij denotes the age of animal *i* at day *j*, that is to say: $\mathbf{age}_{ij} = \text{enter\_age}_i + j$. Functions $\varphi_k$ are the Legendre orthogonal polynomials of degree *k* (Robson, 1959) and the coefficients $a_{ki}$ are the random additive genetic coefficients for animali.

The joint distribution of the random coefficients $\mathbf{a}_k = (a_{ki})_i$, appearing in Eq. (2) and of the random effects $\mathbf{u}_r = (u_{ri})_i$ and $\mathbf{u}_R = (u_{Ri})_i$ of Eq. (1) was the following:

$$ \begin{bmatrix} \mathbf{a}_0 \\ \mathbf{a}_1 \\ \mathbf{a}_2 \\ \mathbf{u}_r \\ \mathbf{u}_R \end{bmatrix} \mathbf{N}(0, \mathbf{G} \otimes \mathbf{A}), \text{ with } \mathbf{G} = \begin{bmatrix} \sigma_{a_0}^2 & \sigma_{a_0 a_1} & \sigma_{a_0 a_2} & \sigma_{a_0 u_r} & \sigma_{a_0 u_R} \\ & \sigma_{a_1}^2 & \sigma_{a_1 a_2} & \sigma_{a_1 u_r} & \sigma_{a_1 u_R} \\ & & \sigma_{a_2}^2 & \sigma_{a_2 u_r} & \sigma_{a_2 u_R} \\ \text{sym} & & & \sigma_{u_r}^2 & \sigma_{u_r u_R} \\ & & & & \sigma_{u_R}^2 \end{bmatrix}, $$

where **A** is the genetic relationship matrix. Residuals $e_{ij}$ were independent and had centered Gaussian distributions with variances at day *j* equal to $11\phi_0(j)^2 + 8\phi_1(j)^2 + 0.6\phi_2(j)^2$. Variances and covariances of the symmetric **G** matrix and formula for the variance of the residuals were chosen in order to obtain moderate heritabilities of the undisturbed BW over time (from 0.20 to 0.35). A detailed description of the parameters/distributions used in the simulation program is provided in Table 1.

The BW in the disturbed situation was modeled by considering the dynamic nature of the response to a perturbation using a dynamic model. Let's consider the following notations: $y_{wd,ij}$ and $y_{d,ij}^*$ are the non-disturbed (wd = without disturbance) and the disturbed BW of animal *i* at day *j* corrected for fixed effects (i.e. $y_{wd,ij}^* = y_{wd,ij} - \mathbf{x}_{ij}\beta$ and $y_{d,ij} = y_{d,ij} - \mathbf{x}_{ij}\beta$), $\Delta_{d,ij} = y_{d,ij} - y_{d,ij-1}^*$ and $\Delta_{wd,i} = y_{wd,i}^* - y_{wd,i-1}^*$, $K_{ij}$ is the number of disturbances that animal *i* was subjected to on day *j* ($K_{ij} \in [0, 3]$), and $int_{k,i} = 0, \dots, K_{ij}$ are their intensities($int_0 = 0$), $int_{ij} = \sum_{k=0}^{K_{ij}} int_k$ (i.e. the sum of the intensities of the disturbances that animal *i* was subjected to on day *j*). At the first time point, $y_{d,11} = y_{wd,11}$ for all animals. Then, for *j* > 1, the disturbed phenotype is recursively defined as:

$$ \frac{\Delta_{d,ij}}{\Delta_{wd,ij}} = 1 + \text{resil}_i \frac{y_{wd,ij}^*}{y_{wd,11}^*} \left( 1 - \frac{y_{d,ij-1}^*}{y_{wd,ij-1}^*} \right) - (1 - \text{resil}_i) \text{int}_{ij}, \quad (3) $$

The right side of the Equation [3] can be split into two terms $A_{ij} = 1 + \text{resil}_i \frac{y_{wd,ij}^*}{y_{wd,11}^*} \left( 1 - \frac{y_{d,ij-1}^*}{y_{wd,ij-1}^*} \right)$ and $B_{ij} = (1 - \text{resil}_i) \text{int}_{ij}$ that depict the response of an animal that is subject to a disturbance. Specifically, $A_{ij}$ takes into consideration the effect of the resilience of animal *i* and $B_{ij}$ the effect of the resistance of animal *i* weighted by the
total intensity of disturbances *int*<sub>ij</sub> experienced by animal *i* at dayj. Indeed, throughout the life of the animal, three situations can arise:

1) Before the first disturbance occurs, $y_{d,ij(−1)} = y_{wd,ij(−1)}$ and $int_{ij} = 0$, thus $A_{ij} = 1$; and $B_{ij} = 0$ so $\Delta_{d,ij} = \Delta_{wd,ij}$ and consequently $y_{d,ij} = y_{wd,ij}$.

2) Then, when the animal is subjected to disturbance(s); during the disturbed period $int_{ij} > 0$, thus $B_{ij}$ is greater than 0. The disturbance is hence considered to have a negative impact on the phenotype proportional to the intensity of the disturbance(s) that is moderated by the resistance of the animal. At the same time, as soon as the disturbed phenotype is lower than the undisturbed one, a resilience mechanism is involved ($A_{ij} > 1$) to also limit the effect of the perturbation.

3) Once the disturbance(s) is/are over, only the resilience mechanism remains ($B_{ij} = 0$ and $A_{ij} > 1$) until the disturbed phenotype reaches the non-disturbed one.

This model to account the effect of disturbances (Eq. (3)) was inspired from the one proposed by Nguyen-Ba et al. (2020):

$$ \frac{\Delta_{d,ij}}{\Delta_{wd,ij}} = 1 + \text{resil}_{NGi} \left( 1 - \frac{y_{d,ij(−1)}}{y_{wd,ij(−1)}} \right) - \text{resis}_{NGi} $$
where $\text{resil}_{NGi}$ and $\text{resis}_{NGi}$ refer to the resistance and resilience, respectively, of animal *i* as defined by Nguyen-Ba et al. (2020).

From the longitudinal non-disturbed and disturbed phenotypes, we extracted two synthetic phenotypes that are measurable in practice in all farms and that were considered of interest for selection: BW<sub>100</sub> and ADG. The impacts of disturbances on genetic parameter estimates and EBV were evaluated for these two synthetic phenotypes of interest. In the non-disturbed situation, only genetic effects affecting growth are involved in the observed phenotype while in the disturbed situation, the observed phenotype may be the result of genetic effects for growth but also for the resistance and resilience.

## Sets of simulations

Various sets of simulations were considered. These different alternatives were obtained by considering:

* Different proportions of pen and individual disturbances (only pen disturbance, only individual disturbance, or both disturbance types),

* Different correlations (negative, neutral, or positive) between the production traits and the two components of the robustness (resistance and resilience),

* Different correlations (negative, neutral, or positive) between resistance and resilience,

* Different heritabilities for resistance and resilience (moderate or low).

The various situations listed above are summarized in Table 2, and they led to 54 different alternatives. For each alternative, 1 000 independent datasets were generated by Monte Carlo simulations. For simplicity, in the rest of the document, the term 'robustness' will refer to resistance and resilience. Then, a positive correlation between production and robustness means a positive correlation between production and resistance and a positive correlation between production and resilience. A moderate heritability of robustness corresponds to situation with moderate heritability for resistance and moderate heritability for resilience and so on.

## Statistical analysis

An animal model was applied for each synthetic phenotype of interest in order to obtain the estimates of heritability and BV in the absence or presence of disturbances: $y_{obs} = X\beta_{obs} + Zu_{obs} + e_{obs}$ where $y_{obs}$ is the vector of ADG or BW<sub>100</sub>, $\beta_{obs}$ the vector of fixed effects with incidence matrix $X$, $u_{obs}$ the genetic additive effects ($u_{obs} N(0, A\sigma_{u_{obs}}^2)$) with incidence matrix $Z$, and $e_{obs}$ the vector of independent residuals $e_{obs} N(0, I\sigma_{e_{obs}}^2)$. Fixed effects included in the model were the contemporary group effect (pen × batch interaction) for both traits and the cross-classified factor with three levels used for the simulation at the age of 100 days for BW<sub>100</sub>.

We then (i) compared the heritabilities obtained in the disturbed and the non-disturbed situations, (ii) calculated the Pearson correlations between BV for the production trait (BV<sub>p</sub>) and EBV; the BV<sub>p</sub> corresponds to the true genetic potential of the animal for growth in the absence of disturbances and it was computed, for individuals $i = 1, \dots, N$, as $\sum_{k=0}^{2} a_{k,i} \phi_k(100)$ for BW<sub>100</sub> and $\sum_{k=0}^{2} (a_{k,i} \phi_k(100) - a_{k,i} \phi_k(1)) / 99$ for ADG, (iii) calculated the partial correlation between EBV and BV for resistance given BV<sub>p</sub> and BV for resilience and the partial correlation between EBV and BV for resilience given BV<sub>p</sub> and BV for resistance, (iv) compared the percentage of animals in common among the best 10% of animals (PB) based on their EBV or BV<sub>p</sub> and the PB based on their EBV obtained in the disturbed or the non-disturbed situation. In addition, (v) considering that 10% of the population was used as the

<table><caption>**Table 2**

Description of the various alternatives of the simulation of the BW of pigs.</caption><thead><tr><th>Item</th><th>Alternatives</th></tr></thead><tbody><tr><td>Proportion of pen (P<sub>pen</sub>) and individual (P<sub>ind</sub>) disturbances</td><td>P<sub>pen</sub> = 40%, P<sub>ind</sub> = 0%, P<sub>pen</sub> = 0%, P<sub>ind</sub> = 40%, P<sub>pen</sub> = 20%, P<sub>ind</sub> = 20%</td></tr><tr><td>Correlation ρ between phenotype of interest and resistance or resilience<sup>1</sup></td><td>σ<sub>a<sub>0</sub>u<sub>r</sub></sub> = 0, σ<sub>a<sub>0</sub>u<sub>R</sub></sub> = 0, σ<sub>a<sub>1</sub>u<sub>r</sub></sub> = 0, σ<sub>a<sub>1</sub>u<sub>R</sub></sub> = 0, σ<sub>a<sub>2</sub>u<sub>r</sub></sub> = 0, σ<sub>a<sub>2</sub>u<sub>R</sub></sub> = 0</td></tr><tr><td>Neutral (ρ = 0)</td><td>ρ<sub>a<sub>0</sub>u<sub>r</sub></sub> = 0.4, ρ<sub>a<sub>1</sub>u<sub>r</sub></sub> = 0.4, ρ<sub>a<sub>2</sub>u<sub>r</sub></sub> = 0.6, ρ<sub>a<sub>2</sub>u<sub>R</sub></sub> = -0.4, ρ<sub>a<sub>2</sub>u<sub>R</sub></sub> = -0.4</td></tr><tr><td>Positive (ρ ≈ 0.4)</td><td>ρ<sub>a<sub>0</sub>u<sub>R</sub></sub> = -0.4, ρ<sub>a<sub>0</sub>u<sub>R</sub></sub> = -0.4, ρ<sub>a<sub>1</sub>u<sub>R</sub></sub> = -0.6, ρ<sub>a<sub>1</sub>u<sub>R</sub></sub> = -0.6, ρ<sub>a<sub>2</sub>u<sub>R</sub></sub> = 0.4, ρ<sub>a<sub>2</sub>u<sub>R</sub></sub> = 0.4</td></tr><tr><td>Negative (ρ ≈ -0.4)</td><td>Neutral ρ = 0</td></tr><tr><td>Correlation (ρ = σ<sub>u<sub>R</sub>u<sub>R</sub></sub>/σ<sub>u<sub>R</sub>u<sub>R</sub></sub>) between resistance and resilience<sup>1</sup></td><td>Positive ρ = 0.4</td></tr><tr><td>Heritability of resistance and resilience<sup>2</sup></td><td>Negative ρ = -0.4</td></tr><tr><td></td><td>Moderate: h<sup>2</sup><sub>resis</sub> = h<sup>2</sup><sub>resil</sub> = 0.29</td></tr><tr><td></td><td>Low: h<sup>2</sup><sub>resis</sub> = h<sup>2</sup><sub>resil</sub> = 0.16</td></tr></tbody></table>

Abbreviations: $a_0$, $a_1$, $a_2$ = genetic regression coefficients of the Legendre polynomials; $u_r$, $u_R$ = breeding values for resistance and resilience; $h_{resis}^2$, $h_{resil}^2$ = heritability of resistance (resis) and resilience (resil).

<footnote n="1">¹ When the correlation between resistance and resilience is negative, the sign of the correlation between the trait of interest and the resistance is the inverse of the sign of the correlation between the trait of interest and the resilience.</footnote>

<footnote n="2">² Heritability on the underlying scale</footnote>
parent of the next generation, for the G5 of each simulated data set (G5) we compared the mean of the true BV<sub>p</sub> in the disturbed and the non-disturbed case for the selected animals (i.e. the best 10% of animals based on their EBV). One-sided paired Z-tests were used at the α risk of 5% for the comparisons of criteria obtained in the disturbed and the non-disturbed situation, and t-tests were used to compare the average difference of the various criteria (without disturbance-with disturbance) between the various simulated alternatives.

## Results

To illustrate the simulated response of animals to a disturbance, the median BW over time of groups of animals undergoing three different situations are depicted in Fig. 1: the group of animals does not face disturbance, the group of animals undergoes a batch disturbance of 'moderate' intensity (1.81) from day 20 to 35 and the group of animals undergoes a batch disturbance of 'weak' intensity (0.52) from day 61 to 82. When there is a disturbance, the maximal difference between the disturbed and the non-disturbed phenotype was observed on the last day of the perturbation (it was 27.54% on day 35 for the 'moderate' case and 6.37% on day 82 for the 'weak' case). Once the perturbation was over, due to

the resilience of the animal, ADG was greater for disturbed animals than ADG of undisturbed animals so that their BW could return to their non-disturbed curve. At the population level, with the parameters used in this study to simulate phenotypes, on average, the impact of the disturbance led to a reduction of the slope of the curve by 40% during the disturbed period. Consequently, at the end of the disturbed period, the value of the disturbed phenotype was, on average, 90% of the non-disturbed one (increasing with the age of the animal), which represents a moderate impact on the phenotype. In all of the simulations sets, the correlation between the disturbed and the non-disturbed phenotype was high for ADG and BW<sub>100</sub> (0.89 ± 0.01 and 0.81 ± 0.02, respectively).

All statistic criteria are presented in Table 3 for the disturbed and the non-disturbed situations in order to understand the impact of disturbances on the estimated genetic parameters and BV. The estimated heritability obtained in the disturbed situation was significantly lower than that obtained in the non-disturbed situation for both traits. The average difference (without disturbance-with disturbance) between heritabilities was0.08 ± 0.02, and 0.05 ± 0.02 for the BW<sub>100</sub> and ADG, respectively. The Pearson correlations between BV<sub>p</sub> and EBV were significantly weaker for the disturbed phenotypes than for the non-disturbed ones for the two phenotypes of interest: the average differences

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Illustration of the effect of disturbances on the growth curve of pigs. In black: median BW of the population that is not disturbed. In blue: median BW of a group of animals undergoing a batch disturbance of 'moderate' intensity (1.81) from day 20 to day 35 (blue dotted lines). In red: median BW of a group of animals undergoing a batch disturbance of 'low' intensity (0.52) from day 61 to day 82 (red dotted lines).</figcaption></figure>
<caption>**Table 3**</caption>

Impact of disturbances on estimated genetic parameters and breeding values for production traits in pigs.

<table><thead><tr><th>Phenotype</th><th>Criterion<sup>1</sup></th><th>Without disturbance<sup>2</sup></th><th>With disturbance<sup>2</sup></th><th>P-value paired Z-test</th></tr></thead><tbody><tr><td rowspan="6">BW<sub>100</sub></td><td>h<sup>2</sup></td><td>0.29 ± 0.03</td><td>0.21 ± 0.03</td><td>&lt; 0.01</td></tr><tr><td>ρ<sub>EBV BV<sub>p</sub></sub></td><td>0.69 ± 0.02</td><td>0.65 ± 0.03</td><td>&lt; 0.01</td></tr><tr><td>P<sub>EBV BV<sub>p</sub></sub></td><td>0.47 ± 0.03</td><td>0.43 ± 0.03</td><td>0.02</td></tr><tr><td>μBV<sub>p</sub></td><td>3.07 ± 1.09</td><td>2.87 ± 1.11</td><td>0.22</td></tr><tr><td>ρ<sub>EBV u<sub>i</sub> BV<sub>p</sub> u<sub>R</sub></sub></td><td>0.00 ± 0.05</td><td>0.04 ± 0.06</td><td>0.28</td></tr><tr><td>ρ<sub>EBV u<sub>i</sub> BV<sub>p</sub> u<sub>R</sub></sub></td><td>0.00 ± 0.05</td><td>0.04 ± 0.06</td><td>0.24</td></tr><tr><td rowspan="6">ADG</td><td>h<sup>2</sup></td><td>0.26 ± 0.03</td><td>0.21 ± 0.03</td><td>&lt; 0.01</td></tr><tr><td>ρ<sub>EBV BV<sub>p</sub></sub></td><td>0.71 ± 0.03</td><td>0.65 ± 0.04</td><td>&lt; 0.01</td></tr><tr><td>P<sub>EBV BV<sub>p</sub></sub></td><td>0.46 ± 0.03</td><td>0.44 ± 0.04</td><td>0.15</td></tr><tr><td>μBV<sub>p</sub></td><td>0.04 ± 0.01</td><td>0.04 ± 0.01</td><td>0.42</td></tr><tr><td>ρ<sub>EBV u<sub>i</sub> BV<sub>p</sub> u<sub>R</sub></sub></td><td>0.00 ± 0.06</td><td>0.03 ± 0.07</td><td>0.20</td></tr><tr><td>ρ<sub>EBV u<sub>i</sub> BV<sub>p</sub> u<sub>R</sub></sub></td><td>0.00 ± 0.06</td><td>0.03 ± 0.07</td><td>0.19</td></tr></tbody></table>

Abbreviations: BW<sub>100</sub> = weight at 100 days of age; ADG = average daily gain.

¹ h², ρ<sub>EBV BV<sub>p</sub></sub>, P<sub>EBV BV<sub>p</sub></sub>, μBV<sub>p</sub>, ρ<sub>EBV u<sub>i</sub> BV<sub>p</sub> u<sub>R</sub></sub> and ρ<sub>EBV u<sub>i</sub> BV<sub>p</sub> u<sub>R</sub></sub> are the heritability, the Pearson correlations between the simulated breeding values (BV) for the production trait (BV<sub>p</sub>) and the estimated breeding values of the observed phenotypes (EBV), the percentage of animals in common among the best 10% of animals based on their EBV or BV<sub>p</sub>, the mean of the simulated BV<sub>p</sub> for the best 10% of animals of the last generation of the phenotypes of interest, the partial correlation between EBV and the simulated BV for resistance given BV<sub>p</sub> and the simulated BV for resilience, and the partial correlation between EBV and the simulated BV for resilience given BV<sub>p</sub> and the simulated BV for resistance.

² μ ± sd calculated on all simulation sets

being 0.04 ± 0.01, and 0.07 ± 0.02 for the BW<sub>100</sub> and ADG, respectively. The partial correlation between EBV and BV for resistance or BV for resilience did not differ significantly between disturbed and non-disturbed situations for either traits. The PB based on their EBV or BV<sub>p</sub> was significantly higher in the non-disturbed situation compared to the disturbed one for the BW<sub>100</sub>, the difference being 0.04 ± 0.02. For both traits, we did not observe any significant differences between the disturbed and the non-disturbed situation for the mean BV<sub>p</sub> of the animals selected in the G5 based on their EBV. The PB based on their EBV between the non-disturbed and the disturbed situation were 0.73, and 0.81 for the BW<sub>100</sub> and ADG, respectively, and were significantly lower than one.

The impact of the disturbance on the estimated heritability and the EBV varied according to the correlation between robustness and production (Table 4). The decreases in heritability, correlations between EBV and BV<sub>p</sub>, and the PB based on EBV or BV<sub>p</sub> between the non-disturbed and the disturbed situation were significantly higher when the genetic correlation between robustness and production was negative compared to when it was positive for both

traits. More specifically, the heritability in the disturbed case was lower in the case of a negative correlation between robustness and production (0.19 and 0.20, for BW<sub>100</sub> and ADG respectively) compared to the positive case (0.22 and 0.23, for BW<sub>100</sub> and ADG respectively). The same trend was observed for the correlation between the BV<sub>p</sub> and EBV (0.65 and 0.64 in the negative case, and 0.67 and 0.66 in the positive case, for BW<sub>100</sub> and ADG respectively), and the PB based on their BV<sub>p</sub> or EBV (0.43 and 0.43 in the negative case, and 0.44 and 0.45 in the positive case, for BW<sub>100</sub> and ADG respectively). The partial correlation between EBV and BV for resistance or BV for resilience tended to increase with the correlation between robustness and production in the disturbed situation for BW<sub>100</sub> (by 4 points between the positive and negative correlation). No significant differences in the change in the mean BV<sub>p</sub> of animals selected in the G5 between the disturbed and the non-disturbed situation were observed depending on the correlation between robustness and production for both traits.

The changes in the different criteria depending on the correlation between resistance and resilience, the type of disturbances,

<caption>**Table 4**</caption>

Average difference between the non-disturbed and the disturbed situation for the various criteria according to the correlation between robustness and production of pigs.

<table><thead><tr><th rowspan="2">Criterion<sup>1</sup></th><th rowspan="2">Phenotype</th><th colspan="3">Correlation between robustness and production</th></tr><tr><th>Negative</th><th>Neutral</th><th>Positive</th></tr></thead><tbody><tr><td rowspan="2">h<sup>2</sup></td><td>BW<sub>100</sub></td><td>0.09 ± 0.01 <sup>a</sup></td><td>0.08 ± 0.01 <sup>ab</sup></td><td>0.07 ± 0.01<sup>b</sup></td></tr><tr><td>ADG</td><td>0.06 ± 0.00 <sup>a</sup></td><td>0.05 ± 0.00<sup>b</sup></td><td>0.03 ± 0.00<sup>c</sup></td></tr><tr><td rowspan="2">ρ<sub>EBV BV<sub>p</sub></sub></td><td>BW<sub>100</sub></td><td>0.04 ± 0.00 <sup>a</sup></td><td>0.03 ± 0.00<sup>b</sup></td><td>0.02 ± 0.00<sup>c</sup></td></tr><tr><td>ADG</td><td>0.07 ± 0.00 <sup>a</sup></td><td>0.06 ± 0.00 <sup>ab</sup></td><td>0.05 ± 0.00<sup>b</sup></td></tr><tr><td rowspan="2">P<sub>EBV BV<sub>p</sub></sub></td><td>BW<sub>100</sub></td><td>0.05 ± 0.01 <sup>a</sup></td><td>0.04 ± 0.00 <sup>ab</sup></td><td>0.03 ± 0.00<sup>b</sup></td></tr><tr><td>ADG</td><td>0.03 ± 0.00 <sup>a</sup></td><td>0.02 ± 0.00 <sup>ab</sup></td><td>0.01 ± 0.00<sup>b</sup></td></tr><tr><td rowspan="2">μBV<sub>p</sub></td><td>BW<sub>100</sub></td><td>0.20 ± 0.04 <sup>a</sup></td><td>0.17 ± 0.04 <sup>a</sup></td><td>0.14 ± 0.04 <sup>a</sup></td></tr><tr><td>ADG</td><td>0.00 ± 0.00 <sup>a</sup></td><td>0.00 ± 0.00 <sup>a</sup></td><td>0.00 ± 0.00 <sup>a</sup></td></tr><tr><td rowspan="2">ρ<sub>EBV u<sub>i</sub> BV<sub>p</sub> u<sub>R</sub></sub></td><td>BW<sub>100</sub></td><td>-0.01 ± 0.01 <sup>a</sup></td><td>-0.04 ± 0.01<sup>b</sup></td><td>-0.05 ± 0.01<sup>b</sup></td></tr><tr><td>ADG</td><td>-0.02 ± 0.01 <sup>a</sup></td><td>-0.03 ± 0.01 <sup>a</sup></td><td>-0.02 ± 0.01 <sup>a</sup></td></tr><tr><td rowspan="2">ρ<sub>EBV u<sub>i</sub> BV<sub>p</sub> u<sub>R</sub></sub></td><td>BW<sub>100</sub></td><td>-0.02 ± 0.01 <sup>a</sup></td><td>-0.04 ± 0.01<sup>b</sup></td><td>-0.06 ± 0.01<sup>b</sup></td></tr><tr><td>ADG</td><td>-0.02 ± 0.01 <sup>a</sup></td><td>-0.03 ± 0.01 <sup>a</sup></td><td>-0.02 ± 0.01 <sup>a</sup></td></tr></tbody></table>

Abbreviations: BW<sub>100</sub> = weight at 100 days of age; ADG = average daily gain.

¹h², ρ<sub>EBV BV<sub>p</sub></sub>, P<sub>EBV BV<sub>p</sub></sub>, μBV<sub>p</sub>, ρ<sub>EBV u<sub>i</sub> BV<sub>p</sub> u<sub>R</sub></sub> and ρ<sub>EBV u<sub>i</sub> BV<sub>p</sub> u<sub>R</sub></sub> are the heritability, the Pearson correlations between the simulated breeding values (BV) for the production trait (BV<sub>p</sub>) and the estimated breeding values of the observed phenotypes (EBV), the percentage of animals in common among the best 10% of animals based on their EBV or BV<sub>p</sub>, the mean of the simulated BV<sub>p</sub> for the best 10% of animals of the last generation of the phenotypes of interest, the partial correlation between EBV and the simulated BV for resistance given BV<sub>p</sub> and the simulated BV for resilience, and the partial correlation between EBV and the simulated BV for resilience given BV<sub>p</sub> and the simulated BV for resistance.

μ ± sd of the average of the sets concerned, sets with a negative correlation between resistance and resilience were not taken into account for these comparisons.

a,b Values within a row with different superscripts differ significantly at *P* < 0.05.
and the heritability of robustness are presented in the *Supplementary Materials* (Table S1, S2, and S3). For both traits, we did not observe any significant differences or trend in the impact of disturbance for the various criteria depending on these sets of parameters.

## Discussion

To model the effect of a disturbance on a longitudinal phenotype, the dynamic model used (Eq. (3)) was inspired from the one proposed by Nguyen-Ba et al. (2020). Both models consider that when a disturbance occurs, it will trigger both the mechanism of resistance and resilience of the animal at the same time. Once there is no longer a disturbance, as long as the ratio between the observed and theoretical curve is smaller than 1, the resilience of the animal will always have an impact on the evolution of the phenotype. It should be noted that both models assume that the negative impact of the disturbance on the phenotype induces a decrease in the value of the phenotype. The model proposed by Nguyen-Ba et al. (2020) is limited in its number of parameters to estimate to prevent identifiability issues. Since the purpose of our model is to simulate data, additional parameters were included, thereby leading to more flexibility. The resistance and resilience of an animal have been used in the present study on a 0–1 scale and included in the model such that the higher the resistance or resilience, the lower the negative response to the perturbation. In Eq. (3), the resilience of the animal was multiplied by $\frac{y_{wd,yj}}{y_{wd,j1}}$ to account for the change in the mean value of the phenotype over time. Indeed, at constant resilience, the capacity of an animal to reach the state before the manifestation of the disturbance is not the same at the beginning as at the end of the test period. It is widely known for example that, it will be harder for a pig to gain 1 kg of weight at the beginning than at the end of the test period. Furthermore, our model considers that the resistance mechanism involved in the response to a perturbation can vary when a disturbance occurs, as it is a time-specific function of both the intrinsic resistance potential of the animal and the intensity of the perturbation, while the resistance mechanism is supposed to be constant during the time window in which the perturbation occurs in the model proposed by Nguyen-Ba et al. (2020). In our model, the value of the resistance and resilience were considered as stable throughout the life of the animal since the duration of the test period of each animal is relatively short (~100 days). Nevertheless, the model can be expanded to account for varying resistance and resilience that may be dependent on the animal's age or experience. Indeed, the resistance (or resilience) of the animal may decrease over time due to exhaustion if it is subject to repeated disturbances, or conversely, in the case of disease responses, it may increase through the development of acquired immunity. We modeled the robustness by the combination of resistance and resilience in order to depict the response of the animal for a given phenotype over a short period of time. However, robustness is a complex concept and several other definitions have been proposed in the literature. For example, Friggens et al. (2017) defined the robustness of an animal by 'the ability, in the face of environmental constraints, to carry on doing the various things that the animal needs to do to favour its future ability to reproduce'. Using this framework, modeling the robustness should consider multiple phenotypes throughout the life of the animal.

In our study, the ADG and BW₁₀₀ were extracted from longitudinal BW for each individual. This choice was guided by the main (non-longitudinal) phenotypes under selection encountered in breeding (Martin et al., 2014; Tortereau et al., 2020). Nonetheless, the results obtained are certainly applicable to any other phenotypes measured at a given time point, as long as, following a disturbance

(whatever the impact of the disturbance on the trait, decrease or increase), the animal can eventually recover its original trajectory. Even if the results were obtained in an oversimplified pig system, it is unlikely that the results would have been different when considering a more complex system (overlapping breeding scheme) since the perturbations essentially add noises to the phenotype measurement.

We have simulated three different types of disturbances that mimic those that can occur in practice. Disturbances at the batch level comprise, for example, a change in temperature, feeding regime, or group prophylaxis. The problem of cannibalism and moderately contagious lameness can be considered to affect close conspecifics and thus occurring at the pen level. Finally, metabolic and non-contagious lameness should affect single individuals. In all three sets of disturbance types, 100% of the batches were subjected to batch disturbance but with an intensity that can be very low. Indeed, it is reasonable to consider that the overall environmental conditions of an animal cannot be optimal during the entire test period (i.e. the temperature can change, the caretakers intervene, etc.).

Different correlations between resistance (resilience) and production were investigated in the present study to cover the various possibilities that can occur. Indeed, there has been no consensus to date regarding the genetic relationship between these traits. For instance, Friggens et al. (2017) have deduced from the resource allocation theory that the correlation between the components of robustness and production should be negative. However, several experiments have shown that it is not always true. The relationship between production and response of the animal to disturbance has been reported in studies evaluating the response of animals to sanitary challenges. By working on the response of pigs to porcine reproductive and respiratory syndrome, Hess et al. (2016) have demonstrated a negative correlation between growth and viremia level (which implies a favorable relationship between production and robustness). Whereas Heckendorn et al. (2017) have demonstrated the opposite by working on the resistance of goats to gastrointestinal nematodes. There is also no consensus regarding the heritability of measurements related to the robustness: Hess et al. (2016) showed a moderately high heritability for the viral load (between 0.31 and 0.51) by working on pigs' response to porcine reproductive and respiratory syndrome, while Mazé-Guilmo et al. (2014) demonstrated that the heritability for resistance (or tolerance) should be low (between 0.18 and 0.19) by working on Leuciscus burdigalensis (a freshwater fish) parasitized by Tracheiastes polycolpus.

The simulated disturbances were of low intensities in order to evaluate the practical impact of disturbances on actual selection. The results obtained show that, between the disturbed and non-disturbed case, there is a significant difference in the heritability, EBV, and in the selected candidates. The impact of the disturbance on the selection of a production trait is exacerbated in the case of a negative correlation between robustness and production compared to the positive case due to a higher reranking of animals (production of the best animals is more affected by a disturbance than the worst animals). Inclusion of a contemporary group effect in the model that corrects for common environmental factors in a group is not sufficient to account for disturbances that occur at the group (batch or pen) level (i.e. there was no difference between the type of perturbation sets on the various criteria) because, since this effect only corrects for the deviation of the mean performance of the group, it does not take into account the variability in the response of the animals to a disturbance due to their variable robustness. Thus, even in a well-controlled environment as observed in conventional breeding program, the selection for production traits can be impacted by disturbances. Considering that the EBV in such a situation corresponds to a combination of the
robustness and production, the breeders could then select for these two terms in the same time. However, this study showed that the EBV in the disturbed situation is not a satisfactory combination of production and robustness: the partial correlation between the EBV and the simulated BV of resistance (and resilience) did not differ from the non-disturbed and disturbed case. In addition, these values were not significantly different from zero. Therefore, we do not select animals for production potential and robustness at the same time in that case. The reason for this is probably because the animals are subjected to random disturbances of small magnitude and that may appear long before the measurement of the phenotype. Additional simulations, to be presented elsewhere, have been carried out to confirm these hypotheses. Because animals raised in production herds are in less controlled and more various environments (lower-quality feed sources, higher sanitary pressure, older facilities, conventional and organic herds, etc.) compared to nucleus herd, it is important to select animals that will be able to adapt to these conditions, which implies an improvement of robustness while maintaining a high level of production. Thus, production and robustness have to be included in the breeding goal. Given our results, to simultaneously improve production and robustness, one should create stronger disturbed environment common to all animals in nucleus herds. However, this may not be ethically acceptable and too expensive for the breeders (Loïc Flatres-Grall, AXIOM company, personal communication). Thus, a more desirable alternative is to obtain separate EBV for growth and robustness. In order to do so, dealing with longitudinal data now available thanks to the development of electronic phenotyping tools should be useful. Such approach has already been investigated by several authors. For instance, Revilla et al. (2021) proposed to quantify the robustness of the animals by the deviation from the observed and theoretical curve of growth of piglets. However, reconstruction of the theoretical production curve is difficult when the time of the perturbation is unknown and further research is needed.

## Conclusion

The results of this study show that random weak disturbances, typically encountered in a breeding farm, have an impact on the estimates of the genetic parameters and BV for production traits and on selection. However, conversely to expected, the EBV obtained in such situation do not correspond to a combination of production and robustness abilities that would be interesting to select. Thus, to select accordingly to a predefined breeding goal for production and robustness, it would be interesting to separate these two components by studying the dynamic of the phenotype over time instead of measuring the trait at a given time point. This subject is currently undergoing further research.

## Supplementary material

Supplementary data to this article can be found online at https://doi.org/10.1016/j.animal.2022.100496.

## Ethics approval

Not applicable.

## Data and model availability statement

The simulation program was not deposited in an official repository.

It is available upon request to the authors.

## Author ORCIDs

Vincent Le: https://orcid.org/0000-0001-7645-5878

Ingrid David: https://orcid.org/0000-0002-2514-6693

Tom Rohmer: https://orcid.org/0000-0002-4751-2324

## Author contributions

Vincent Le: Formal analysis, Writing - Original Draft.

Ingrid David: Conceptualization, Supervision, Methodology, Software, Writing - Reviewing and Editing.

Tom Rohmer: Formal analysis, Supervision, Writing- Reviewing and Editing.

## Declaration of interest

The authors declare that there are no conflicts of interest.

## Acknowledgements

None.

## Financial support statement

This study was funded 50% by the genetic animal division of INRAE and 50% by Alliance R&D

## References

* Colditz, I.G., Hine, B.C., Colditz, I.G., Hine, B.C., 2016. Resilience in farm animals: biology, management, breeding and implications for animal welfare. Animal Production Science 56, 1961-1983.

* Frey, M., Hofer, A., Künzi, N., 1997. Comparison of models with a fixed or a random contemporary group effect for the genetic evaluation for litter size in pigs. Livestock Production Science 48, 135-141.

* Friggens, N.C., Blanc, F., Berry, D.P., Puillet, L., 2017. Review: deciphering animal robustness. A synthesis to facilitate its use in livestock breeding and management. Animal 11, 2237-2251.

* de Goede, D.M., Gremmen, B., Blom-Zandstra, M., 2013. Robust agriculture: balancing between vulnerability and stability. NJAS - Wageningen Journal of Life Sciences 64-65, 1-7.

* Heckendorn, F., Bieber, A., Werne, S., Saratsis, A., Maurer, V., Stricker, C., 2017. The genetic basis for the selection of dairy goats with enhanced resistance to gastrointestinal nematodes. Parasite (Paris, France) 24, 32.

* Hess, A.S., Islam, Z., Hess, M.K., Rowland, R.R.R., Lunney, J.K., Doeschl-Wilson, A., Plastow, G.S., Dekkers, J.C.M., 2016. Comparison of host genetic factors influencing pig response to infection with two North American isolates of porcine reproductive and respiratory syndrome virus. Genetics, Selection, Evolution 48, 43.

* Knap, P.W., 2005. Breeding robust pigs. Australian Journal of Experimental Agriculture 45, 763-773.

* Kolmodin, R., Strandberg, E., Madsen, P., Jensen, J., Jorjani, H., 2002. Genotype by environment interaction in Nordic dairy cattle studied using reaction norms. Acta Agriculturae Scandinavica, Section A-Animal Science 52, 11-24.

* Martin, P., Raoul, J., Bodin, L., 2014. Effects of the FecL major gene in the Lacaune meat sheep population. Genetics Selection Evolution 46, 48.

* Mazé-Guilmo, E., Loot, G., Páez, D.J., Lefèvre, T., Blanchet, S., 2014. Heritable variation in host tolerance and resistance inferred from a wild host-parasite system. Proceedings. Biological Sciences 281, 20132567.

* Meyer, K., 2009. Factor-analytic models for genotype × environment type problems and structured covariance matrices. Genetics Selection Evolution 41, 21.

* Nguyen-Ba, H., van Milgen, J., Taghipoor, M., 2020. A procedure to quantify the feed intake response of growing pigs to perturbations. Animal 14, 253-260.

* Phocas, F., Belloc, C., Bidanel, J., Delaby, L., Dourmad, J.Y., Dumont, B., Ezanno, P., Fortun-Lamothe, L., Foucras, G., Frappat, B., González-García, E., Hazard, D., Larzul, C., Lubac, S., Mignon-Grasteau, S., Moreno, C.R., Tixier-Boichard, M., Brochard, M., 2016. Review: towards the agroecological management of ruminants, pigs and poultry through the development of sustainable breeding programmes: I-selection goals and criteria. Animal 10, 1749-1759.

* Rauw, W.M., Gomez-Raya, L., 2015. Genotype by environment interaction and breeding for robustness in livestock. Frontiers in Genetics 6, 310.

* Revilla, M., Lenoir, G., Flatres-Grall, L., Muñoz-Tamayo, R., Friggens, N.C., 2021. Quantifying growth perturbations over the fattening period in swine via mathematical modelling. Peer Community Journal 2, e9.

* Robson, D.S., 1959. A simple method for constructing orthogonal polynomials when the independent variable is unequally spaced. Biometrics 15, 187-191.
Schaeffer, L.R., 2004. Application of random regression models in animal breeding. Livestock Production Science 86, 35-45.

Shariati, M.M., Korsgaard, I.R., Sorensen, D., 2009. Identifiability of parameters and behaviour of MCMC chains: a case study using the reaction norm model. Journal of Animal Breeding and Genetics = Zeitschrift Fur Tierzuchtung Und Zuchtungsbioologie 126, 92-102.

Star, L., Ellen, E., Uitdehaag, K., Brom, F., 2008. A plea to implement robustness into a breeding goal: poultry as an example. Journal of Agricultural and Environmental Ethics 21, 109-125.

Tortereau, F., Marie-Etancelin, C., Weisbecker, J.-L., Marcon, D., Bouvier, F., Moreno-Romieux, C., François, D., 2020. Genetic parameters for feed efficiency in Romane rams and responses to single-generation selection. Animal 14, 681-687.

Van Vleck, L.D., 1987. Contemporary groups for genetic evaluations. Journal of Dairy Science 70, 2456-2464.