A promising resilience parameter for breeding: the use of weight and feed trajectories in growing pigs

Abstract

Background

Increasing resilience is a priority in modern pig breeding. Recent research shows that general resilience can be quantified via variability in longitudinal data. The collection of such longitudinal data on weight, feed intake and feeding behaviour in pigs has been facilitated by the development of technologies such as automated feeding stations.

The goal of this study was to investigate resilience traits, which were estimated as deviations from longitudinal weight, feed intake and feeding behaviour data during the finishing phase. A dataset with 324,207 records between the age of 95 and 155 days on 5,939 Piétrain pigs with known pedigree and genomic information was used. We provided guidelines for a rigid quality control of longitudinal body weight data, as we found that outliers can significantly affect results. Gompertz growth curve analysis, linear modelling and trajectory analyses were used for quantifying resilience traits.

Results

To our knowledge, this is the first study comparing resilience traits from longitudinal body weight, feed intake and feeding behaviour data in pigs. We demonstrated that the resilience traits are lowly to moderately heritable for deviations in body weight (h² = 2.9%–20.2%), in feed intake (9.4%–23.3%) and in feeding behaviour (16.2%–28.3%). Additionally, these traits have good predictive abilities in cross-validation analyses. Deviations in individual body weight and feed intake trajectories are highly correlated (r_g = 0.78) with low to moderate favourable genetic correlations with feed conversion ratio (r_g = 0.39–0.49). Lastly, we showed that some resilience traits, such as the natural logarithm of variances of observed versus predicted body weights (lnvar_weight), are more robust to lower observation frequencies and are repeatable over three different time periods of the finishing phase.

Conclusions

Our results will help future studies investigating resilience traits and resilience-related traits. Moreover, our study provides first results on standardization of quality control and efficient data sampling from automated feeding station data. Our findings will be valuable for breeding organizations as they offer evidence that pigs’ general resilience can be selected on with good accuracy. Moreover, this methodology might be extended to other species to quantify resilience based on longitudinal data.

Background

Resilience in livestock usually refers to the ability of animals to be minimally affected by (environmental) stressors and/or to cope with these stressors and quickly return to their optimal production level [1,2,3,4,5]. As such, resilience is becoming an important breeding goal in pig breeding [6]. Increasing resilience is particularly interesting as it can simultaneously tackle animal welfare concerns, reduce labor and treatment costs [1, 5,6,7,8]. Moreover, the need for robust, easy-to-handle animals rises with an increased number of animals per farmer. This is evidenced in the European Union, where the average size of pig farms keeps growing [9]. Although the need for more resilient pigs is evident, it has been difficult and/or costly to phenotype informative traits for pigs’ (general) resilience [7]. On one hand, most routinely phenotyped resilience indicators are scored as binary (e.g., ‘dead’ versus ‘alive’) or ordinal (e.g., ‘no’, ‘mild’ or ‘severe’ disease) traits. These traits often have low frequencies, with low variability and low heritabilities [10, 11]. On the other hand, immunological traits, such as viral load or antibody levels, show moderate to high heritabilities and a good association with animal health, but are costly to phenotype and in practice challenging to obtain [12, 13].

Recently however, several studies showed that increasing within-family and within-individual uniformity can improve animals’ general resilience. Blasco et al. [14] and Formoso-Rafferty et al. [15] independently executed two successful selection experiments on respectively litter size uniformity in rabbits and birth weight uniformity in mice. For lines with increased (within-family) uniformity, both studies found a correlated selection response with a higher survival in these uniform lines and a favorable association with disease susceptibility traits. Scheffer et al. [5] and Berghof et al. [1] proposed to derive resilience traits from longitudinal phenotypes by quantifying the variability in longitudinal data. Here, the hypothesis is that animals with a higher within-individual uniformity over time will have a higher resilience as they will show less deviations from their optimal production level in the presence of (environmental) disturbances [1, 5]. Recent studies have reported that these within-individual deviations of longitudinal data are lowly to moderately heritable (Table 1). Moreover, these studies generally found favourable genetic correlations between within-individual uniformity and resilience-related traits, such as mortality and disease incidence. Hence, less deviations in longitudinal data (higher within-individual uniformity), was linked with higher survival and lower disease incidence (Table 1). However, the number of studies investigating this relationship is currently limited.

Table 1 Overview of genetic studies on within-individual trait deviations based on longitudinal data

Thanks to technological developments, longitudinal data can be collected on a large scale in practice [5]. For instance, the use of automatic feeding stations (AFS) enables individual recording of pigs’ feed intake, feeding behaviour (duration and time of visits) and body weight. Despite the elevated cost of AFS, most pig breeding organizations have invested in this technology [1]. In addition, advances in wearables and computer vision systems may create longitudinal data in pigs for a variety of traits [4, 5] including body temperature, respiration rate [26] and activity levels [27]. The integration of genomics and other ‘omics’ techniques could further aid the development of efficient selection programs for increased resilience [7].

In this study we will investigate the genetic background of resilience proxies based on longitudinal body weight, feed intake and feeding behaviour data in a Piétrain pig population. This study is the first to examine the value of body weight deviations based on trajectory analysis as a novel proxy for resilience. Moreover, it is unique in its comparison and analysis of deviations in weight, feed intake and feeding behaviour over time, as previous studies have only focused on deviations in weight, feed intake or feeding behaviour. Lastly, we investigate the influence of observation frequency on the stability of resilience traits and the influence of observation period on the repeatability of resilience traits. Therefore, the repeatability of resilience traits over different stages of the finishing period (observation period) was studied as well as the impact of less data points per individual (observation frequency). Genetic parameters, such as heritability, genetic coefficient of variation and genetic correlations are estimated for the resilience traits. In an effort to better understand the value of genomics in selection for resilience, we assessed the predictive abilities using pedigree relationships or single-step genomic evaluation.

Methods

Animals and data collection

The study was carried out on Piétrain pigs from Hendrix Genetics (Hypor Maxter). The nucleus pig test barn (France) consisted of 14 compartments with 10 pens per compartment and on average 15 pigs per pen (1.0 m² per pig). Water was provided ad libitum in each pen from one nipple drinker and feed was provided with an automatic feeding system (AFS): Nedap pig performance testing feeding station (Nedap N.V.; Groenlo, the Netherlands). Individual recordings of weight (accuracy of 0.5 kg), feed intake (accuracy of 1 g), visit duration (accuracy of 1 s) and number of visits were obtained with the AFS per day. The daily records were calculated as summary statistics based on a pigs’ daily feeding station visits. Before data quality control (QC), the dataset comprised of 7,880 pigs born between May 2017 and September 2021. In total, these pigs had 522,122 AFS recordings for weight, feed intake, feeding duration and number of visits (on average 66 records per pig for each trait). Moreover, for all these pigs with AFS recordings, individual weights were also recorded by technicians at birth, 14 days of age, start of test (81 ± 5 d and 32.6 ± 7.6 kg) and end of test (161 ± 12 d and 114.3 ± 13.1 kg). At the end of test, muscle thickness and fat thickness were measured via ultrasound probing between 3^rd and 4^th last rib using Exago (IMV) device for these pigs.

Quality control

This study investigates variability in longitudinal data from AFS, and links this variability with underlying biological/genetic factors. Therefore, it is vital that variability due to technical errors and/or noise are removed as much as possible. In a first step of quality control, outlier correction limits were designed based on population statistics to identify and exclude gross weight recording inaccuracies. Specifically, AFS weight recordings below 10 kg (n = 874) or above 160 kg before an age of 160 days (n = 9,339) were set to missing (511,909 AFS weight records retained). Additionally, only pigs with at least twenty AFS weight recordings were retained to ensure a sufficient number of records per individual for the accurate estimation of resilience traits. After this first step of QC, 6,831 pigs and 505,990 AFS weight records were retained.

Next, inaccurate AFS weight recordings were identified on a pen level using the root mean square error (RMSE), similar to [20]. RMSE was obtained by linear regression of weight on age. Weight records of individuals in outlying pens were visually inspected (Fig. 1) and treated as follows: (i) erroneous weight recordings within a time period < 20 d were set to missing; (ii) Individuals with erroneous weight recordings over a longer time period were removed from the dataset (6,788 pigs and 501,320 AFS weight records retained). Next, a 10-day rolling median weight was calculated per individual. Weight recordings deviating more than 3 kg from this median rolling weight were considered as outliers and set to missing (Fig. 2; 6,728 pigs and 495,312 AFS weight records retained). Furthermore, pigs with gaps in weight recordings larger than ten days were removed. Hereafter, the RMSE of weight regressed on age was re-calculated, and outlying individuals were checked again. After these QC steps, the dataset contained 6,457 pigs, 439,963 weight recordings (82% of pigs and 84% of records before QC).

Fig. 1

Outlier detection on pen level by analyzing root mean squared error of weight (RMSE_weight). a Histogram of RMSE_weight in function of age on a pen level before quality control. Pens with high RMSE_weight estimates were visually inspected for (technical) errors. b Example of a pen with no severe outlying weights at the pen level, although some individual weight recordings are outlying. Weight evolution of individual pigs are represented with a specific color. c Example of a pen with outlying weights at start of trajectory. Such outliers are often due to an adaptation phase of the pigs, i.e., upon entering the automated feeding station, pigs tend to enter the station with their penmates, inflating the daily weight estimates. Weight evolution of individual pigs are represented with a specific color. d Example of technical issues causing outlying weights and high RMSE_weights. In these cases, outliers were set to missing, or outlying individuals were removed from the dataset. Weight evolution of individual pigs are represented with a specific color

Fig. 2

Example of 10-d rolling median approach combined with second order polynomial regression to detect outliers. Observed weights outlying predicted weight ± 3 kg were considered as outliers (red) and set to missing

For daily feed intake (FI), visit duration and number of visits, values exceeding the average plus four times the standard deviation were set to missing (5,550 g/d for FI; 3.3 h/d for visit duration; 33 visits/d for number of visits). Hereafter, individual and pen RMSE were obtained for these traits by regressing them on age. However, no outliers were detected using this method. After these QC steps, the dataset contained 6,457 pigs with 438,132 feed intake records, 437,753 visit duration records and 436,886 number of AFS visit records.

Next, only AFS records were kept between an age of 95–155 d to standardize age limits across animals. These thresholds were selected because most of our AFS recordings fall within this range (Fig. 3) and because most pigs show a learning curve after entering the pen with AFS, which disappears around d 95 in our dataset. Finally, data were further standardized by removing pigs with (i) starting age > 110 d (n = 75), (ii) maximum age < 120 d (n = 226 pigs), (iii) > 30% missing records for weight or feed intake (n = 240 pigs) and (iv) < 20 d with AFS records (n = 188 pigs). The final dataset after QC comprised 5,939 pigs (5,811 boars and 128 sows) with 324,478 AFS weight recordings, 323,775 feed intake recordings, 323,304 visit duration recordings and 322,910 number of visit recordings between 95 and 155 days of age (75% of pigs and ~ 62% of records before QC). The pigs originated from 1,273 dams and 130 sires (2,105 unique litters). Pedigree consisted of 9,369 pigs with a pedigree depth ranging from 13 to 19 generations. Genomic information (45,436 SNPs) was available for 6,726 pigs in total, of which 5,160 pigs (87% of dataset) had own phenotypic records. The evolution of weight in function of age for data after QC is shown in Fig. 4a.

Fig. 3

Distribution of number of weight records. a Number of individuals with automatic feeding station weight recordings in function of age (d). Red lines indicate thresholds of 95 and 155 d. b Histogram of number of records per pig after selecting the age range of 95–155 d. The maximum amount of records is 60

Fig. 4

Evolution of weight and standardized weight in function of age. a Evolution of weight in kg in function of age in d for the dataset after quality control. The mean weight per age (d) is shown in solid red line, a one standard deviation difference from the mean is shown in dashed red lines. b Standardized weights with a mean of zero and a standard deviation of 1 per age in days. For example, a score of ‘2’ indicates a pig had a weight which was two standard deviations above the mean of the population on that specific age. The mean standardized weight per age (d) is shown in solid red line, a one standard deviation difference from the mean is shown in dashed red lines

Derivation of traits

After QC, traits were operationalized. Average daily gain (ADG) was estimated as

$$ADG\;(kg/d)=\frac{weight\;\left(kg\right)\;at\;maximum\;age\;AFS-weight\;(kg)\;at\;minimum\;age\;AFS}{maximum\;age\;AFS\;\left(d\right)-minimum\;age\;AFS\;(d)}$$

Average feed intake (AFI) per individual was estimated as total feed intake divided by number of days with a feed intake record. Feed conversion ratio (FCR) was estimated by dividing ADG by AFI.

An overview of resilience trait definitions is given in Table 2. A number of resilience traits was operationalized based on deviations in weight trajectories. First, this was established by individually fitting a Gompertz growth curve [28] based on AFS weights between 95–155 days of age, supplemented with birth weight, weight at 14 d, weight at start and end of test. Expected weights were estimated with R [29] using the nls function and the Gompertz growth curve formula (Additional file 1: Fig. S1):

$${weight}_{ij}={A}_{i}\times{e}^{-{B}_{i}\times{e}^{{k}_{i}\times{t}_{ij}}}+{\varepsilon }_{ij}$$

where A_i, B_i and k_i are the growth curve parameters for individual i, t_ij is day j for individual i and ${\varepsilon }_{ij}$ is residual error. For every individual, we quantified lnvar_weight as the natural logarithm of the variance in the daily differences between observed weights versus expected weights via Gompertz modeling (calculated with ln function in R), as well as skewness (skew_weight; calculated with skewness function) and the lag-one autocorrelation (lag1_weight; calculated with acf function) (Fig. 5), following Berghof et al. [1].

Table 2 Trait definition for the resilience traits

Fig. 5

Example of trait construction for two pigs (a–d versus e–h). The upper pig (a–d) showed little deviations in observed versus expected body weight, whereas the lower pig (e–h) showed many deviations in observed versus expected body weight. These examples are the same animals as shown in Fig. 6. a and e Example of Gompertz growth curve modelling on automated feeding station data of individual pigs. The Gompertz growth curve is shown as a solid red line, observed daily weights are given as black dots. b and f Deviations of observed versus predicted weights after Gompertz modeling: lnvar_weight, lag1_weight and skew_weight are estimated based on these deviations. c and g Example of standardized weights with mean zero and standard deviation one for the population on a daily basis. The variance of these standardized weights for an individual was used to calculate lnvar_{weight_standardized}. d and h Trajectory analysis of weight. Here, weight gain/loss is seen as a trajectory from start until end, with age in d as x-coordinate and weight as y-coordinate. From this trajectory, mean speed and straightness were calculated as resilience traits

Next, similar to Putz et al. [16], a linear regression of weight on age was used to estimate the root mean squared error (RMSE) of observed versus expected body weight deviations (calculated with lm function in R). As growing pigs (95–155 d) are more or less in a linear phase of their growth curve [30], this linear approach seems justified (Fig. 5a). Hereafter, we calculated the natural logarithm of the MSE (lnMSE_weight). Here, we used MSE instead of RMSE to make lnMSE_weight equivalent to lnvar_weight.

A major challenge with modelling has to do with some circularity: expected weights are also estimated based on the observed weights, and these “expectations” might come from a biased curve [1, 22]. To circumvent this issue, two different approaches were used. First, following Berghof et al. [24], all weights were standardized by age with a mean of zero and a standard deviation of one for each single day of age (Fig. 4b, Fig. 5c and g). From these standardized weights per day, the natural logarithm of the variance was then calculated (lnvar_{weight_standardized}). Pigs with a high lnvar_{weight_standardized} hence showed great variations in weight over time, compared to the population mean. Second, additional deviation traits were derived from trajectory analysis using the trajr package in R [31] (Fig. 5d and h). Trajectory analysis can be used to estimate deviations from expected patterns. Here, we estimated mean speed (TrajDerivatives function) and the straightness (TrajStraightness function). The trajr package estimates mean speed as:

$$Mean\ speed=\frac{Total\ path\ length\ of\ weight\ trajectory}{Age\ difference\ between\ start\ and\ end\ point}$$

whereas straightness index was estimated as:

$$straightness =\frac{Euclidean\ distance\ between\ start\ and\ end\ point}{Total\ path\ length\ of\ weight\ trajectory}$$

Hence, an animal with more body weight deviations will have a higher mean speed and a lower straightness index, as the total path length of weight trajectory will increase due to more deviations. The maximum straightness index value is 1, with values below one indicating more deviations from a straight line. The straightness index and mean speed are related, but can differ due to different ADG between animals. For example, two animals with a straightness index of 1 might still differ in mean speed, as a faster growing animal will have a higher mean speed as it will have more’distance traveled’ over the same time.

For daily feed intake, visit duration and number of visits, the natural logarithm of MSE after linear modeling was calculated using the same methodology as for lnMSE_weight, respectively leading to the traits lnMSE_FI, lnMSE_dur and lnMSE_{n_visit} (Fig. 6). Moreover, following Putz et al. [16], the number of off-feed days was calculated as the number of days during which feed intake (QR_FI) and/or visit duration (QR_dur) was in the 5% lowest quantile using quantile regression (QR) on age over all pigs (Additional file 2: Fig. S2).

Fig. 6

Example of feed intake, visit duration and number of daily visits for two pigs (a–c versus d–f). These examples are the same animals as shown in Fig. 5, and were selected based on body weight deviations, where the upper pig showed little deviations in observed versus expected body weight, whereas the lower pig showed many deviations in observed versus expected body weight. Red lines indicates the regression line from linear modeling. a and d Evolution of feed intake (kg/d) versus age (d). Based on the linear regression, lnMSE_FI was quantified. b and e Evolution of visit duration (s/d) versus age in d. Based on the linear regression, lnMSE_dur was quantified. c and f Evolution of number of daily visits to feeder versus age (d). Based on the linear regression, lnMSE_{n_visit} was quantified

Finally, after estimating these traits per pig from the daily AFS recordings, estimates deviating by more than four standard deviations from the mean were set to missing (184 for A; 55 for B; 30 for k; 27 for FI; 35 for ADG; 73 for FCR; 2 for lag1_weight; 11 for skew_weight; 0 for lnvar_weight; 1 for lnMSE_weight; 2 for lnvar_{weight_standardized}; 0 for straightness; 7 for mean speed; 2 for lnMSE_FI; 1 for lnMSE_dur; 0 for lnMSE_{n_visit}).

Genetic modelling

The blupf90 suite of programs [32] was used to estimate genetic parameters. Heritability (h²) was estimated as the proportion of additive genetic variance divided by total variance. Likewise, the common environmental effect (c²) was estimated as the proportion of variance explained by random contemporary group effects (c), divided by total variance. For the resilience traits, the genetic coefficient of variation (GCV) was estimated as a measure of evolvability or possible selection response of a trait [33]. For the lnvar and lnMSE traits, GCV was estimated as: $GCV=\sqrt{{\sigma }_{a}^{2}}$, as described by [18, 34, 35] for an exponential model.

All single trait animal models were of the form:

$${\varvec{y}}={\varvec{X}}{\varvec{b}}+{\varvec{Z}}{\varvec{a}}+{\varvec{W}}{\varvec{c}}+{\varvec{e}}$$

where y is the vector with phenotypes for the studied trait; b is the vector containing the fixed effects (sex, 2 levels; farm, 2 levels) and covariates (maximum age); a is the vector of additive genetic effects (9,371 animals in pedigree, 6,723 with genotype information), which is assumed to follow a normal distribution for the pedigree matrix (A) using only pedigree relationships:

$${\varvec{a}}\sim N(0,{{\varvec{A}}\sigma }_{a}^{2})$$

Or a normal distribution for the H matrix, combining both pedigree (A) and genomic (G) relationship matrices following [36,37,38] using single-step genomic evaluation:

$${\varvec{a}}\sim N(0,{{\varvec{H}}\sigma }_{a}^{2})$$

c is the vector of contemporary group effects (113 levels), assumed to follow a normal distribution $c\sim N(0,{{\varvec{I}}\sigma }_{c}^{2})$, with I the identity matrix; e is the vector of residual effects, assumed to follow a normal distribution $e\sim N(0,{{\varvec{I}}\sigma }_{e}^{2})$; X, Z and W are incidence matrices for respectively fixed effects, random animal effects and random contemporary group effects. The random contemporary group effect c is a combination of farm, compartment and date of entrance at farm. Contemporary groups with less than ten pigs were combined in a remainder group (165 pigs).

Likewise, genetic correlations (r_g) between traits were estimated via bivariate animal models of the form:

$$\left[\begin{array}{c}{\varvec{y}}1\\ {\varvec{y}}2\end{array}\right]=\left[\begin{array}{cc}{\varvec{X}}1& 0\\ 0& {\varvec{X}}2\end{array}\right]\left[\begin{array}{c}{\varvec{b}}1\\ {\varvec{b}}2\end{array}\right]+\left[\begin{array}{cc}{\varvec{Z}}1& 0\\ 0& {\varvec{Z}}2\end{array}\right]\left[\begin{array}{c}{\varvec{a}}1\\ {\varvec{a}}2\end{array}\right]+\left[\begin{array}{cc}{\varvec{W}}1& 0\\ 0& {\varvec{W}}2\end{array}\right]\left[\begin{array}{c}{\varvec{c}}1\\ {\varvec{c}}2\end{array}\right]+\left[\begin{array}{c}{\varvec{e}}1\\ {\varvec{e}}2\end{array}\right]$$

Similar to the single-trait animal model, y1 and y2 are the vectors with phenotypes for the studied traits; b1 and b2 are the vectors containing the fixed effects and covariates; a1 and a2 are the vectors of additive genetic effects, which is assumed to follow a normal distribution for the H matrix using single-step genomic evaluation:

$$\left[\begin{array}{c}{\varvec{a}}1\\ {\varvec{a}}2\end{array}\right]\sim N(\left[\begin{array}{c}0\\ 0\end{array}\right],\left[\begin{array}{cc}{\sigma }_{a1}^{2}& {\sigma }_{a1,a2}\\ {\sigma }_{a1,a2}& {\sigma }_{a2}^{2}\end{array}\right]\otimes {\varvec{H}})$$

c1 and c2 are the vectors of contemporary group effects (113 levels), assumed to follow a normal distribution $\left[\begin{array}{c}{\varvec{c}}1\\ {\varvec{c}}2\end{array}\right]\sim N(\left[\begin{array}{c}0\\ 0\end{array}\right],\left[\begin{array}{cc}{\sigma }_{c1}^{2}& {\sigma }_{c1,c2}\\ {\sigma }_{c1,c2}& {\sigma }_{c2}^{2}\end{array}\right]\otimes {\varvec{I}})$; e1 and e2 are the vector of residual effects, assumed to be independently normal distributed $\left[\begin{array}{c}{\varvec{e}}1\\ {\varvec{e}}2\end{array}\right]\sim N(\left[\begin{array}{c}0\\ 0\end{array}\right],\left[\begin{array}{cc}{\sigma }_{e1}^{2}& 0\\ 0& {\sigma }_{e2}^{2}\end{array}\right])$; X1, X2, Z1, Z2, W1 and W2 are incidence matrices for respectively fixed effects, random animal effects and random contemporary group effects.

Cross validation

Cross validation was performed using three data masking strategies: within family masking, across family masking and temporal masking. For within and across family masking, we decided to use 5-fold cross-validation with random masking of 20% of the data (based on [39]), resulting in validation datasets of ~ 1,200 pigs. Moreover, ten replications were used to avoid random sampling effects [39], resulting in 50 models (10 × 5-fold validation) per trait. For within family masking strategy, one out of five offspring was randomly masked per sire. In our dataset, every sire had a mean of 45.7 offspring (range = 1–182, SD = 42.9), leading to a mean of 9 masked offspring per sire. For the across family masking strategy, all progeny from one out of five sires was randomly masked. As a result, the within-family strategy is valuable to estimate predictive ability from close relationships, whereas across-family masking allows to estimate predictive ability of distant relationships [39]. For the temporal cross-validation strategy, animals born after 01-10-2020 were masked (~ 30% of dataset). Temporal cross-validation allows to estimate forward (future) predictive ability.

The process of cross-validation was as follows. First, a univariate animal model (as specified before) was used on the full dataset using the remlf90 software. Observed phenotypes were adjusted for fixed and non-genetic random effects based on these results using the predictf90 software:

$${{\varvec{y}}}^{*}={\varvec{y}}-(\widehat{{\varvec{b}}}+\widehat{{\varvec{c}}} )=\widehat{{\varvec{a}}}+\widehat{{\varvec{e}}}$$

Predictive abilities were estimated as the Pearson correlation between breeding values of a validation dataset (with masked phenotypes) and the adjusted phenotypes ($y^\ast$ ):

$$Predictive\ ability=r({EBV}_{masked},{y}^{*})$$

Next, predictive abilities were expressed as a cross validation accuracy. This was done by dividing the predictive ability by the square root of the estimated h²:

$$Predictive\ ability\ accuracy=\frac{r\left({EBV}_{masked},{y}^{*}\right)}{h}$$

Evaluating the impact of observation frequency and observation period

Finally, the impact of observation frequency and observation period were evaluated. First, the influence of observation frequency on parameter stability was evaluated. Based on the full dataset, subsets were made with 1 out of 4 records per animal (~ 2 records per week), 1 out of 7 records per animal (1 record per week) and 1 out of 14 records per animal (1 record every two weeks) (Fig. 7). Phenotypic and genetic correlations were estimated for the full model versus reduced datasets using bivariate animal models. These (genetic) correlations indicate to what extent traits are sensitive to changes in observation frequency. Second, to assess influence of observation period, the full 60-day dataset with all records was divided in three age groups of twenty days: (i) 95–115 days of age (early), (ii) 115–135 days of age (middle) and (iii) 135–155 days of age (late). Based on these subsets, all traits were recalculated leading to, for example, lnvar_weight-early, lnvar_{weight-middle} and lnvar_weight-late. Hereafter, bivariate animal models were run within each trait to estimate phenotypic and genetic correlations between periods. These (genetic) correlations indicate the repeatability of a trait and whether a given trait genetically shifts over time.

Fig. 7

Example of different observation frequency and observation period settings for an individual pigs’ weight data. a All daily weight records within the 95–155 days of age interval, colored per observation period 95–115 d (early, red), 115–135 d (middle, orange), 135–155 d (late, green). b A subset sampled from the full dataset with only 1 out of 4 data points, which corresponds to about two records per week. c A subset sampled from the full dataset with only 1 out of 7 data points, which corresponds to about one record per week. d A subset sampled from the full dataset with only 1 out of 14 data points, which corresponds to about one record every two weeks

Results

An overview of the main trait distributions and their estimated h² and variance components is given in Table 3. All estimated phenotypic and genetic correlations are given in Table 4. Heritabilities for Gompertz growth curve parameters A, B and k were low (6.8%–10.3%). The body weight deviation traits skew_weight (2.9%) and lag1_weight (6.2%) were also lowly heritable. Standardizing weights before estimating lnvar increased h² (12.1% for lnvar_{weight_standardized} versus 11.0% for lnvar_weight). However, h² estimates for body weight deviations were highest for trajectory parameters straightness (15.5%) and mean speed (20.2%), which were moderately heritable. Deviations related to feed intake and feeding behaviour had higher h² (20.7%–28.3%) than body weight deviations (8.9%–20.2%). Despite low to moderate h², the resilience trait indicators had high genetic coefficients of variation: 20.5%–30.2% for body weight deviations and 29.1%–33.4% for feed intake and feeding behaviour deviations. QR_FI (h² = 9.4%) and QR_dur (h² = 16.1%) had low to moderate h² estimates.

Table 3 Descriptive statistics and genetic parameters

Table 4 Correlation table with phenotypical (below diagonal) and genetic correlations (above diagonal)

Phenotypic and genetic correlations (Table 4) between lnvar_weight and most other body weight deviation traits were high (r_p = 0.58–0.88; r_g = 0.53–0.93), except for skew_weight (r_p = 0.01; r_g = 0.32). Furthermore, lnvar_weight was phenotypically and genetically also moderately to highly correlated with deviations in feeding duration (lnMSE_dur; r_p = 0.49; r_g = 0.36) and feed intake (lnMSE_FI; r_p = 0.69; r_g = 0.78). Deviations in feed intake (lnMSE_FI) were moderately correlated with deviations in feeding duration (lnMSE_dur; r_p = 0.69, r_g = 0.49) but lowly correlated with deviations in number of daily visits (lnMSE_{n_vis}; r_p = 0.16, r_g = −0.34). ADG was moderately correlated with straightness (r_p = 0.32; r_g = 0.38) and mean speed (r_p = 0.32; r_g = 0.50). Additionally, ADG was negatively correlated with the number of days with a very low feed intake (QR_FI; r_p = −0.51, r _g = −0.70), indicating that pigs with high ADG have less off-feed days. A similar pattern was observed for AFI. For FCR, a low to moderate favourable correlation was found with lnvar_weight (r_p = 0.17; r_g = 0.37) and lnMSE_FI (r_p = 0.20; r_g = 0.49), indicating that pigs with more deviations in weight and feed intake have a higher FCR.

An overview of estimated predictive abilities per trait using both pedigree relationships and single-step genomic evaluation for three cross-validation strategies is given in Table 5 as cross validation accuracy and in Additional file 4: Table S2 as correlation.

Table 5 Predictive ability accuracy for cross validation scenarios: masking across or within family and temporal masking

Predictive ability accuracies of skew_weight were low for all strategies (0.00–0.23). For the body weight deviation traits, the trajectory parameters mean speed and straightness showed the highest predictive ability accuracies with single-step genomic evaluation (0.38–0.60). Feed intake deviations showed higher predictive ability accuracies than body weight deviations, and single-step genomic evaluation seemed to relatively increase predictive abilities for feed intake deviations more. Predictive abilities for lnMSE_dur, for example, increased by 54%, 21% and 33% respectively when adding genomics to across, within and temporal masking strategy.

Phenotypic correlations per trait for different observation frequencies are given as pairwise correlation plots in Additional file 5: Fig. S3. An overview of genetic correlations within traits estimated by using different observation frequencies ranging from 1 in 4 to 1 in 14 is provided in Table 6. As expected, ADG does not change substantially with lower data density (r_p = 0.85 and r_g = 0.95 with 1 in 14 density), as it is estimated as the average gain over a long period. FCR fluctuates more with lower observation frequency: when considering only 1 in 14 data points, the phenotypic and genetic correlations with the full dataset drop (r_p = 0.45 and r_g = 0.76) and h² drops from 22.1% to 10.1%. For the resilience traits, lnvar and lnMSE estimates were least dependent on observation frequency with r_p = 0.44–0.76 and r_g = 0.79–0.96 in the most extreme scenario, although the h² estimates decreased substantially from h² = 10.6%–23.3% to h² = 5.1%–10.1%. Skew_weight and lag1_weight were strongly impacted by differences in observation frequency, with r_p = 0.05–0.08 and r_g = 0.02–0.14 in the most extreme scenario. The trajectory parameters mean speed and straightness were moderately affected by data density (r_p = 0.29–0.33 and r_g = 0.50 with one in 14 data points), but showed a smaller decrease in h² estimate from h² = 15.0%–21.4% to h² = 13.1%–17.7%.

Table 6 Genetic parameters of full dataset versus reduced datasets (1 in x data points)

The influence of observation period was studied by dividing the full 60-day dataset (95–155 days of age) in three 20-day time periods during the finishing phase (early, middle and late). Phenotypic correlations per trait over time periods are given as pairwise correlation plots in Additional file 7: Fig. S4. Genetic correlations for each time period versus the full dataset within each trait are given in Additional file 6: Table S3. For ADG and FCR, early, middle and late estimates are moderately to highly correlated with the full dataset (respectively r_p = 0.59–0.63; r_g = 0.71–0.82 and r_p = 0.48–0.55; r_g = 0.65–0.85), although genetic correlations are low to moderate within time periods (respectively r_p = 0.04–0.17; r_g = 0.24–0.41 and r_p = −0.03–0.05; r_g = 0.34–0.66). This is in contrast to AFI, where correlations were also moderate to high between time periods (r_p = 0.45–0.68; r_g = 0.65–0.87). The body weight deviation traits lnvar_weight, lnMSE_weight, lnvar_{weight_standardized}, straightness and mean speed show high correlations between time periods and the full dataset (r_p = 0.47–0.74; r_g = 0.63–0.92) and moderate to high correlations within time periods (r_p = 0.23–0.50; r_g = 0.45–0.78). In contrast, lag1_weight and skew_weight show in general low correlations over time periods (r_p = −0.02–0.11; r _g = −0.15–0.35). Feed intake deviations lnMSE_FI, lnMSE_duration and lnMSE_{n_visit} showed moderate to high (genetic) correlations (r_p = 0.32–0.65; r_g = 0.73–0.97).

Discussion

Increasing resilience is becoming priority in modern pig breeding [1, 6]. Therefore, this study investigated resilience traits based on weight, feed intake and feeding behaviour in pigs which were estimated as perturbations in longitudinal data. We demonstrate that these resilience traits are lowly to moderately heritable and have good predictive abilities in cross-validation analyses. Moreover, deviations in individual body weight and feed intake trajectories are genetically highly correlated and show low to moderate favourable genetic correlations with feed conversion ratio. Lastly, we show that the observation frequency and observation period impact some resilience traits more severely than others. lnvar_{weight_standardized} and lnMSE_FI, for example, were more robust to low observation frequencies (as low as one data point in fourteen days) and showed moderate repeatability over three 20-day time periods of the finishing phase.

In the first part of our study, we quantified and evaluated several resilience traits. The body weight deviation traits lnvar_weight, skew_weight, lag1_weight were based on [1, 5] after Gompertz growth curve modelling, whereas lnvar_{weight_standardized} was based on Berghof et al. [24] after standardizing weights per age. The main difference between the two lnvar traits is that lnvar_weight uses the pigs’ individual data as a reference (based on growth curve modelling), whereas lnvar_{weight_standardized} takes the population statistics as a reference. The deviations in feed intake and behaviour (lnMSE_FI, lnMSE_dur, lnMSE_{n_visit}, QR_FI, QR_dur) were based on Putz et al. [16], although we chose to use MSE instead of RMSE, as this allowed us to directly estimate GCV [18]. In addition to these previously described resilience traits, we deducted resilience traits from linear modelling and trajectory analysis to our weight data in the finishing phase of pigs (lnMSE_weight, straightness, mean speed). We believe this approach is justified, as an expected weight evolution in the finishing phase of pigs is more or less linear [30]. Our hypothesis is that any deviation from this linear trajectory is probably due to an external challenge which can impact a pigs’ optimal production potential and challenges its resilience. Although trajectory analysis was developed for the analysis of (wild) animals’ actual trajectories in time and space [31], we believe this methodology could be translated to weight patterns of finishing pigs. The start weight of a pig can be regarded as the starting point, following a specific path over time to reach an end weight. Moreover, trajectory analysis is appealing as it does not require any complex modelling of expected (weight) trajectories. The main issue with modelling is that the predicted values tend to follow the observed values, complicating the prediction of the optimal production curve for challenged animals [1, 22]. Figure 5e, for example, shows that the modelled Gompertz growth curve is more or less the mean of the observed values, which results in an overestimation of positive deviations, and an underestimation of negative deviations [1].

To our knowledge, this is the first study to report h² for body weight deviation traits lnvar_weight, lnvar_{weight_standardized} and lnMSE_weight in pigs (Table 3). Our h² estimates range 8.9%–12.1% for these traits, which is similar to the h² estimate of 9%–11% in similar body weight deviations in layer chickens [24]. h² and GCV for lnvar_{weight_standardized} were higher (h² = 12.1%; GCV = 30.2%) compared to lnvar_weight (h² = 11.0%; GCV = 21.6%) and lnMSE_weight (h² = 8.9%; GCV = 20.5%). This might be because lnvar_{weight_standardized} corrects for a scaling effect, since changes in mean levels tend to change variance levels as well [1, 40, 41]. Remarkably, straightness and mean speed had slightly higher h² estimates (15.5% and 20.2%), whereas h² of lag1_weight (2.9%) and skew_weight (6.2%) was very low, similar to Poppe et al. [23]. The estimated h ² of feed intake deviations (QR_FI, lnMSE_FI; h² = 9.4%–23.3%) and feeding behaviour deviations (QR_dur, lnMSE_dur, lnMSE_{n_visit}; h² = 16.2%–28.3%) were also comparable to previous studies in pigs by Putz et al. [16] (h² = 8%–26% for feed intake), Homma et al. [18] (h² = 31% for feed intake; h² = 36%–40% for feeding behaviour), and Kavlak and Uimari [19] (h² = 7%–11% for feed intake; h² = 16%–20% for feeding behaviour). Estimated GCV for lnvar_weight and lnMSE_weight were 21%–22% and were lower than GCV estimates of 29%–33% for lnvar_{weight_standardized}, lnMSE_FI, lnMSE_dur and lnMSE_{n_visit}, but in the same range as (22%–39%) [18]. These high genetic coefficients of variation indicate a large potential for genetic improvement of these traits [1, 35].

There are no standard guidelines yet on how to perform quality control of weight data from AFS. The quality control procedure in the current paper was based on the structure and identified issues from our dataset, combined with the methodology from previous work [20]. We would like to stress the importance of rigid quality control on an individual level when quantifying resilience traits, especially for body weight deviations. In contrast to feed intake and feeding behaviour, weight is accumulated over time, i.e. you can only gain or lose weight gradually. However, erroneous weights, such as sudden drops and rises, do often appear in raw data from AFS. These errors can be technical (machine error) or due to a learning curve of the pigs after introduction to AFS [20] (Fig. 1). Without any quality control, estimated h² for lnvar_weight, straightness and mean speed were very low (h² = 1.8%–4.0%; results not shown). Applying a limited quality control on a population level, for example applying minimum and maximum thresholds for weight as a function of age, increased h² estimates to h² = 5.7%–7.1% (results not shown). However, these estimates are still considerably lower than what is achieved in a dataset with a rigid, individual quality control. Here, standard guidelines on quality control of AFS data might be valuable, although there might be no “one size fits all” approach. Our advice is to always visually check the weight trajectories of individual animals with outlying resilience traits, for example lnvar_weight > 3 standard deviations from mean, even after quality control.

The data in Table 4 suggests a strong connection between resilience traits for feed intake and weight, as shown by the estimated genetic correlation of 0.78 between lnvar_weight and lnMSE_FI. This correlation implies that individual deviations in feed intake are rapidly reflected in weight perturbations. However, the correlation does not equal one, indicating that these various indicators of resilience may signify different aspects of pigs’ resilience. Here, changes in feed intake might be considered as a short term response to a challenge, as a challenged animals’ appetite is usually directly affected [1, 2]. Variations in weight can be considered as a moderate term response since weight gain/loss is mainly determined by food and water intake and several other factors over time. Moreover, we estimated a favourable genetic correlation between lnMSE_weight or lnMSE_FI, and FCR (r_g = 0.39–0.49). As feed efficiency is one of the most important traits in pig breeding, this favourable correlation would facilitate an implementation of resilience traits into breeding programs. Correlations between lnvar_weight and most other body weight deviation traits were high (r_p = 0.58–0.88; r_g = 0.53–0.93) except for skewness (r_p = 0.01; r_g = 0.32). These correlations indicate different traits mostly capture the same genetic variation, but some differences exist between traits. Since the weight trajectory parameters straightness and mean speed showed higher h² and do not rely on complex modeling, these traits might be more interesting to implement in breeding programs. Additionally, straightness has a favourable genetic correlation with FCR (r_g = −0.41) and ADG (r_g = 0.38). Notably, lnMSE_{n_visits} was lowly to negatively correlated with lnMSE_FI (r_p = −0.05, r_g = −0.34) and lnMSE_dur (r_p = 0.16, r_g = 0.00). Similar genetic correlations were found by [18]. These findings might imply that more deviations in daily visits to feeding station do not necessarily lead to more variation in the time spent at the AFS and might even reduce deviations in feed intake which is counterintuitive.

It should be noted that our data were collected in purebred pigs in a high health breeding farm. This is in contrast to commercial crossbred finishing pigs, which are typically raised in a more challenging environment with, for example, a higher disease pressure and more social stressors such as a higher pig density. The commercial conditions might elicit more easily differences in resilience [1]. Nonetheless, our data show considerable heritable variation for resilience traits with reasonable predictive ability. However, the purebred-crossbred correlation (r_pc) of these resilience traits in pigs is not yet known. Research on this topic is essential for pig breeding programs, as an r_pc < 0.80 indicates crossbred information should be taken into account [42]. For example, in a study on egg production data in layer chicken, an r_pc was estimated ranging from 0.16–0.47 (lnvar of egg production) to 0.56–0.63 (lag1 autocorrelation) [25]. Furthermore, the main limitation of the present study is that we could not corroborate our resilience traits with resilience related factors such as mortality, disease prevalence, treatments, etc., as done by Putz et al. [16].

Predictive ability analysis using three masking strategies indicated good prospects for selection on most resilience traits (Table 5). The across family masking strategy generally yielded lower predictive abilities compared to the within family masking strategy, as family relationships are more distant in the across family masking strategy. Moreover, adding genotypes to the analysis in general improved predictive abilities. Interestingly, trajectory parameters straightness and mean speed yielded the highest predictive abilities for body weight deviations, demonstrating their potential use for breeding programs. Moreover, resilience indicators for feed intake and feeding behaviour yielded higher predictive abilities. Using single-step genomic evaluation generally improved predictive ability, mainly for the across-family (average increase of +62.2%) and temporal (+67.8%) masking strategy compared to within-family masking (+13.2%). This was expected, as these masking strategies use more distant family information, without own phenotypes and, hence, adding extra genomic information relatively improves predictive ability more [39]. Sae-Lim et al. [40] previously showed that predictive ability of (untransformed) body weight uniformity in salmon could be improved by adding genomic information.

As indicated by Berghof et al. [1], the frequency of observations and observation length are crucial to determine good resilience traits. In our study, we used daily recordings from AFS over a 60-day period within a pigs’ finishing phase (95–155 d). However, AFS may be used more efficiently and/or a limited number of manual weight recordings might be a suitable alternative. Moreover, AFS have not yet been developed and generally used for many livestock species. Therefore, we examined the influence of frequency of observations (Table 6 and Additional file 5: Fig. S3) and length of observation period (Additional file 6: Fig. S4 and Additional file 7: Table S3) by using different data densities and by splitting the dataset in three 20-day periods. If only one record every two weeks or daily records for a short time period would be informative for some resilience traits, these observations could also be collected manually. For example, Berghof et al. [24] used seven weight recordings with a 4-week interval in layer chickens, whereas [43] only had five manual weight recordings of Nile tilapia over a 162-day period. Another option would be to more efficiently use the expensive technology (e.g., AFS), by rotating it over animals so it can be used more efficiently, or by only recording a shorter observation period, although this might pose practical/sanitary issues in pigs. Interestingly, lnvar_{weight_standardized} seemed to be very stable with r_p > 0.76 and r_g > 0.96 between full dataset and only one weight recording every two weeks (±5 records in total). These results reiterate the need for data standardization, particularly for traits with a changing average and variance over time such as weight. Whereas trajectory parameters straightness and mean speed seem to have highest h² and predictive ability for body weight deviations, these traits are also more sensitive to low data densities, with r_p = 0.29–0.33 and r_g = 0.50 for 1 in 14 data density compared to the full dataset. Further, lnMSE_FI showed to be quite stable with lower data densities with r_p = 0.44 and r_g = 0.84 between full data and 1 in 14 scenario. Moreover, phenotypic and genetic correlations for deviations in feed intake were high over different time periods, with r_p = 0.48 and r_g = 0.80 between lnMSE_{FI_early} and lnMSE_{FI_late}, and r_p = 0.73–0.87 and r_g = 0.90–0.97 between 20-day time periods and the total 60-d period. These results show that, similar to FI, feed intake deviations are moderately repeatable over time: pigs with a high variability in feed intake at the start of the finishing phase, will generally also have a high variability in feed intake at the end of the finishing phase. Observational period and frequency had a large impact on skew_weight and lag1_weight. Therefore, these indicators might not be useful for data with a low observation frequency and/or a short observation period.

In light of our findings, we provided suggestions on the choice of resilience traits to include in a breeding program. The inclusion of resilience traits based on feed intake and feeding behaviour deviations show to be most promising, with highest h², GCV and predictive ability. Additionally, these traits seem to be robust to changes in observation frequency and period. However, our study also suggests to include body weight deviations as resilience indicator in breeding programs, as the (genetic) correlations with feed intake and feeding behaviour resilience traits substantially differed from one. We hypothesize that body weight deviations reflect more moderate term responses to external challenges, whereas feed intake and feeding behaviour better reflect short term responses to external stressors. For body weight deviation traits, we recommend to perform a rigid quality control of body weights, as we found that outliers can significantly affect results. Although we provide some guidelines for QC of AFS body weight data, most studies currently still perform an ad-hoc QC. Future work on (more) uniform guidelines for QC could further improve standardization and replicability of results across studies. Regarding quality control of body weights based on AFS data, future studies should focus on more uniform guidelines. We also recommend standardizing weights over time. Finally, the trajectory analysis traits straightness and mean speed showed promise as body weight resilience traits as they had the highest h² and predictive ability and a favourable (genetic) correlation with FCR. However, these traits seem more sensitive to observation frequency.

Conclusions

To our knowledge, this is the first study comparing resilience traits from longitudinal body weight, feed intake and feeding behaviour data in pigs. We showed these resilience traits are lowly to moderately heritable (h² = 3%–28%) with good predictive abilities. Moreover, we suggested new, promising resilience indicators based on trajectory analysis with higher h² and predictive ability, although these traits were more sensitive to observation frequency. Next, we were the first to report the influence of observation frequency and observation period on these resilience traits and showed that feed intake and feeding duration deviations are very robust to low data density and moderately repeatable over time. Within body weight deviation traits, lnvar_{weight_standardized} seemed most robust to low data density, stressing the need for weight standardization over age when quantifying body weight deviations. Our results can help the design of future studies to look at the relationship between these resilience traits and resilience-related traits such as mortality and disease incidence, and to estimate the purebred-crossbred correlation. We believe our findings will be very useful for pig breeding programs, and will aid in the improvement of pigs’ general resilience by selective breeding. We recommend the inclusion of resilience indicators from both feed intake and body weight deviations in breeding programs, as they could offer valuable insights into different aspects of pigs' resilience. Moreover, we are confident our methodology can be extended to other species as well.