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

Wim Gorssen , Carmen Winters, Roel Meyermans, Léa Chapard, Katrijn Hooyberghs, Steven Janssens,Abe Huisman, Katrijn Peeters, Han Mulder and Nadine Buys*

Abstract

Keywords Deviations, Genetics, Gompertz growth curves, Heritability, Pigs, Predictive ability, Resilience, Trajectory analysis

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-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-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 ‘a(chǎn)live’) 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.

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 m2per 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 3rdand 4thlast 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).

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.1 Outlier detection on pen level by analyzing root mean squared error of weight (RMSEweight).a Histogram of RMSEweight in function of age on a pen level before quality control.Pens with high RMSEweight 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 RMSEweights.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

Derivation of traits

After QC, traits were operationalized.Average daily gain(ADG) was estimated as 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 thenlsfunction and the Gompertz growth curve formula(Additional file 1: Fig.S1):

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

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

Average feed intake (AFI) per individual was estimated as total feed intake divided by number of days with a feed whereAi,Biandkiare the growth curve parameters for individuali, tijis dayjfor individualiand εijis residual error.For every individual, we quantified lnvarweightas the natural logarithm of the variance in the daily differences between observed weights versus expected weights via Gompertz modeling (calculated withlnfunction in R),as well as skewness (skewweight; calculated withskewnessfunction) and the lag-one autocorrelation (lag1weight; calculated withacffunction) (Fig.5), following Berghof et al.[1].

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 withlmfunction 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 (lnMSEweight).Here, we used MSE instead of RMSE to make lnMSEweightequivalent to lnvarweight.

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 (lnvarweight_standardized).Pigs with a high lnvarweight_standardizedhence showed great variations in weight over time, compared to the population mean.Second, additional deviation traits were derived from trajectory analysis using thetrajrpackage in R [31] (Fig.5d and h).Trajectory analysis can be used to estimate deviations from expected patterns.Here, we estimated mean speed (TrajDerivativesfunction) and the straightness (TrajStraightnessfunction).Thetrajrpackage estimates mean speed as:

whereas straightness index was estimated as:

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 straightnessindex 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.

Table 2 Trait definition for the resilience traits

Table 2 (continued)

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 lnMSEweight, respectively leading to the traits lnMSEFI,lnMSEdurand lnMSEn_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(QRFI) and/or visit duration (QRdur) was in the 5% lowest quantile using quantile regression (QR) on age over all pigs (Additional file 2: Fig.S2).

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 forA; 55 forB; 30 fork; 27 for FI; 35 for ADG; 73 for FCR; 2 for lag1weight; 11 for skewweight; 0 for lnvarweight; 1 for lnMSEweight; 2 for lnvarweight_standardized; 0 for straightness; 7 for mean speed; 2 for lnMSEFI; 1 for lnMSEdur; 0 for lnMSEn_visit).

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: lnvarweight, lag1weight and skewweight 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 lnvarweight_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

Genetic modelling

whereyis the vector with phenotypes for the studied trait;bis the vector containing the fixed effects (sex, 2 levels; farm, 2 levels) and covariates (maximum age);ais 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:

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, lnMSEFI was quantified.b and e Evolution of visit duration (s/d) versus age in d.Based on the linear regression, lnMSEdur was quantified.c and f Evolution of number of daily visits to feeder versus age (d).Based on the linear regression, lnMSEn_visit was quantified

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

cis the vector of contemporary group effects (113 levels), assumed to follow a normal distribution c ～N(0,Iσ2c) , withIthe identity matrix;eis the vector of residual effects, assumed to follow a normal distribution e ～N(0,Iσ2e) ;X,ZandWare incidence matrices for respectively fixed effects, random animal effects and random contemporary group effects.The random contemporary group effectcis 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 (rg) between traits were estimated via bivariate animal models of the form:

Similar to the single-trait animal model,y1andy2are the vectors with phenotypes for the studied traits;b1andb2are the vectors containing the fixed effects and covariates;a1anda2are the vectors of additive genetic effects, which is assumed to follow a normal distribution for theHmatrix using single-step genomic evaluation:

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 off-spring (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 theremlf90software.Observed phenotypes were adjusted for fixed and non-genetic random effects based on these results using thepredictf90software:

Predictive abilities were estimated as the Pearson correlation between breeding values of a validation dataset (with masked phenotypes) and the adjusted phenotypes ( 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 estimatedh2:

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, lnvarweight-early, lnvarweight-middleand lnvarweight-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.

Results

An overview of the main trait distributions and their estimatedh2and variance components is given in Table 3.All estimated phenotypic and genetic correlations are given in Table 4.Heritabilities for Gompertz growth curve parametersA,Bandkwere low (6.8%-10.3%).The body weight deviation traits skewweight(2.9%) and lag1weight(6.2%) were also lowly heritable.Standardizing weights before estimating lnvar increasedh2(12.1% for lnvarweight_standardizedversus 11.0% for lnvarweight).However,h2estimates 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 higherh2(20.7%-28.3%) than body weight deviations (8.9%-20.2%).Despite low to moderateh2, 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.QRFI(h2= 9.4%) and QRdur(h2= 16.1%) had low to moderateh2estimates.

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

Phenotypic and genetic correlations (Table 4) between lnvarweightand most other body weight deviation traits were high (rp= 0.58-0.88;rg= 0.53-0.93), except for skewweight(rp= 0.01;rg= 0.32).Furthermore, lnvarweightwas phenotypically and genetically also moderately to highly correlated with deviations in feeding duration(lnMSEdur;rp= 0.49;rg= 0.36) and feed intake (lnMSEFI;rp= 0.69;rg= 0.78).Deviations in feed intake (lnMSEFI)were moderately correlated with deviations in feeding duration (lnMSEdur;rp= 0.69,rg= 0.49) but lowly correlated with deviations in number of daily visits (lnMSEn_vis;rp= 0.16,rg= -0.34).ADG was moderately correlated with straightness (rp= 0.32;rg= 0.38) and mean speed (rp= 0.32;rg= 0.50).Additionally, ADG was negatively correlated with the number of days with a very low feed intake (QRFI;rp= -0.51,rg= -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 lnvarweight(rp= 0.17;rg= 0.37) and lnMSEFI(rp= 0.20;rg= 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.

Predictive ability accuracies of skewweightwere 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 lnMSEdur, 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 (rp= 0.85 andrg= 0.95with 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 (rp= 0.45 andrg= 0.76) andh2drops from 22.1% to 10.1%.For the resilience traits, lnvar and lnMSE estimates were least dependent on observation frequency withrp= 0.44-0.76 andrg= 0.79-0.96 in the most extreme scenario, although theh2estimates decreased substantially fromh2= 10.6%-23.3% toh2= 5.1%-10.1%.Skewweightand lag1weightwere strongly impacted by differences in observation frequency, withrp= 0.05-0.08 andrg= 0.02-0.14 in the most extreme scenario.The trajectory parameters mean speed and straightness were moderately affected by data density(rp= 0.29-0.33 andrg= 0.50 with one in 14 data points),but showed a smaller decrease inh2estimate fromh2= 15.0%-21.4% toh2= 13.1%-17.7%.

Table 3 Descriptive statistics and genetic parameters

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 (respectivelyrp= 0.59-0.63;rg= 0.71-0.82 andrp= 0.48-0.55;rg= 0.65-0.85),although genetic correlations are low to moderate within time periods (respectivelyrp= 0.04-0.17;rg= 0.24-0.41 andrp= -0.03-0.05;rg= 0.34-0.66).This is in contrast to AFI, where correlations were also moderate to high between time periods (rp= 0.45-0.68;rg= 0.65-0.87).The body weight deviation traits lnvarweight, lnMSEweight,lnvarweight_standardized, straightness and mean speed showhigh correlations between time periods and the full dataset (rp= 0.47-0.74;rg= 0.63-0.92) and moderate to high correlations within time periods (rp= 0.23-0.50;rg= 0.45-0.78).In contrast, lag1weightand skewweightshow in general low correlations over time periods (rp= -0.02-0.11;rg= -0.15-0.35).Feed intake deviations lnMSEFI,lnMSEdurationand lnMSEn_visitshowed moderate to high(genetic) correlations (rp= 0.32-0.65;rg= 0.73-0.97).

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

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.lnvarweight_standardizedand lnMSEFI, 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 lnvarweight, skewweight, lag1weightwere based on[1, 5] after Gompertz growth curve modelling, whereas lnvarweight_standardizedwas based on Berghof et al.[24]after standardizing weights per age.The main difference between the two lnvar traits is that lnvarweightuses the pigs’ individual data as a reference (based on growth curve modelling), whereas lnvarweight_standardizedtakes the population statistics as a reference.The deviations in feed intake and behaviour (lnMSEFI, lnMSEdur, lnMSEn_visit,QRFI, QRdur) 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 (lnMSEweight,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].

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

To our knowledge, this is the first study to reporth2for body weight deviation traits lnvarweight, lnvarweight_standard-izedand lnMSEweightin pigs (Table 3).Ourh2estimates range 8.9%-12.1% for these traits, which is similar to theh2estimate of 9%-11% in similar body weight deviations in layer chickens [24].h2and GCV for lnvarweight_standard-izedwere higher (h2= 12.1%; GCV = 30.2%) compared to lnvarweight(h2= 11.0%; GCV = 21.6%) and lnMSEweight(h2= 8.9%; GCV = 20.5%).This might be because lnvarweight_standardizedcorrects 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 higherh2estimates (15.5% and 20.2%),whereash2of lag1weight(2.9%) and skewweight(6.2%) was very low, similar to Poppe et al.[23].The estimatedh2of feed intake deviations (QRFI,lnMSEFI;h2= 9.4%-23.3%)and feeding behaviour deviations (QRdur,lnMSEdur,lnMSEn_visit;h2= 16.2%-28.3%) were also comparable to previous studies in pigs by Putz et al.[16] (h2= 8%-26%for feed intake), Homma et al.[18] (h2= 31% for feed intake;h2= 36%-40% for feeding behaviour), and Kavlak and Uimari [19] (h2= 7%-11% for feed intake;h2= 16%-20% for feeding behaviour).Estimated GCV for lnvarweightand lnMSEweightwere 21%-22% and were lower than GCV estimates of 29%-33% for lnvarweight_standardized,lnMSEFI, lnMSEdurand lnMSEn_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, estimatedh2for lnvarweight,straightness and mean speed were very low (h2= 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, increasedh2estimates toh2= 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 lnvarweight＞ 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 lnvarweightand lnMSEFI.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 lnMSEweightor lnMSEFI,and FCR (rg= 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 lnvarweightand most other body weight deviation traits were high (rp= 0.58-0.88;rg= 0.53-0.93) except for skewness (rp= 0.01;rg= 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 higherh2and 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 (rg= -0.41) and ADG (rg= 0.38).Notably, lnMSEn_visitswas lowly to negatively correlated with lnMSEFI(rp= -0.05,rg= -0.34) and lnMSEdur(rp= 0.16,rg= 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 (rpc)of these resilience traits in pigs is not yet known.Research on this topic is essential for pig breeding programs, as anrpc＜ 0.80 indicates crossbred information should be taken into account [42].For example, in a study on egg production data in layer chicken, anrpcwas 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, lnvarweight_standardizedseemed to be very stable withrp＞ 0.76 andrg＞ 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 highesth2and predictive ability for body weight deviations, these traits are also more sensitive to low data densities, withrp= 0.29-0.33 andrg= 0.50 for 1 in 14 data density compared to the full dataset.Further, lnMSEFIshowed to be quite stable with lower data densities withrp= 0.44 andrg= 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, withrp= 0.48 andrg= 0.80 between lnMSEFI_earlyand lnMSEFI_late, andrp= 0.73-0.87 andrg= 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 skewweightand lag1weight.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 highesth2, 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 highesth2and 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(h2= 3%-28%) with good predictive abilities.Moreover,we suggested new, promising resilience indicators based on trajectory analysis with higherh2and 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, lnvarweight_standardizedseemed 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.

Abbreviations

A,BandkGompertz growth curve parameters

ADG Average daily gain

AFI Average feed intake

AFS Automated feeding station

BLUP Best linear unbiased prediction

c2Common environmental effect

EBV Estimated breeding value

FCR Feed conversion ratio

GCV Genetic coefficient of variation

h2Heritability

lag1weightLag1 autocorrelation of observed versus predicted weight distribution

lnMSEdurNatural logarithm of mean squared error of visit duration in function of age

lnMSEFINatural logarithm of mean squared error of feed intake in function of age

lnMSEn_visitNatural logarithm of mean squared error of number of daily visits in function of age

lnMSEweightNatural logarithm of mean squared error of weight in function of age

lnvarweightNatural logarithm of variance of observed versus predicted weights

lnvarweight_standardizedNatural logarithm of variance of standardized weights

Mean speed Mean speed of weight in function of age after trajectory analysis

QC Quality control

QRdurNumber of days with visit duration below 5% of quantile after quantile regression

QRFINumber of days with feed intake below 5% of quantile after quantile regression

rgGenetic correlation

RMSE Root mean square error

rpPhenotypic correlation

skewweightSkewness of observed versus predicted weight distribution

SNP Single nucleotide polymorphism

ssGBLUP Single-step genomic best linear unbiased prediction

straightness Straightness index of weight in function of age after trajectory analysis

y*Adjusted phenotype

σaAdditive genetic standard deviation

σcCommon environmental standard deviation

σeResidual standard deviation

Supplementary Information

The online version contains supplementary material available at https:// doi.org/ 10.1186/ s40104- 023- 00901-9.

Additional file 1: Fig.S1.Gompertz growth curve distribution and parameters.

Additional file 2: Fig.S2.Quantile regressionof feed intake and visit duration.

Additional file 3: Table S1.Genetic correlations between all trait combinations using bivariate models.

Additional file 4: Table S2.Predictive abilities as Pearson correlation between masked breeding values and corrected phenotype.

Additional file 5: Fig.S3.Pairwise correlation plots for all evaluated traits with full datasets and reduced datasets.

Additional file 6: Table S3.Estimated genetic correlations between the full dataset and reduced datasets.

Additional file 7: Fig.S4.Pairwise correlation plots for all evaluated traits with full datasets and reduced datasets.

Acknowledgements

We would like to acknowledge Manuel Revilla and Jordi Vila Teixidor for their help on this study.We also thank the editor and reviewers for their valuable input.

Authors’ contributions

WG analysed the data and wrote the manuscript.WG, SJ, HM, AH, KP and NB designed and conceived this study.WG, CW, LC, KH, RM, SJ, HM, AH, KP and NB critically reviewed the analyses and the manuscript.All authors read and approved the final manuscript.

Funding

This study was partially funded by an FR PhD fellowship (1104320N; WG) and two SB PhD fellowships (1S05818N (CW) and 1S37119N (RM)) of the Research Foundation Flanders (FWO).Moreover, RM and LC were also partly funded by a KU Leuven C2 project (C24/18/036).KH was funded by the UNIPIG project of VLAIO(HBC.2019.2866).The funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Availability of data and materials

The datasets generated and/or analysed during the current study are not publicly available due to data restriction from Hendrix-Genetics but are available from the corresponding author on reasonable request and with permission of Hendrix Genetics.

Declarations

Ethics approval and consent to participate

Data on the pigs were collected according to Hendrix Genetics protocols,under the supervision of Hendrix Genetics employees.Data were collected as part of routine animal data collection in a commercial breeding program for pigs, and therefore ethical approval was not necessary.

Consent for publication

Not applicable.

Competing interests

KP and AH are employees of Hendrix-Genetics and provided the data for this study, although Hendrix Genetics did not fund this study.Moreover,the funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.All authors declare that the results are presented in full and as such present no conflict of interest.

Author details

1Center for Animal Breeding and Genetics, Department of Biosystems, KU Leuven, Kasteelpark Arenberg 30 - Box 2472, 3001 Leuven, Belgium.2Laboratory for Biological Psychology, KU Leuven, Tiensestraat 102 - Box 3714, 3000 Leuven, Belgium.3Hendrix Genetics, P.O.Box 114, 5830 AC Boxmeer, The Netherlands.4Research Animal Breeding and Genomics, Wageningen University, P.O.Box 338, 6700 AH Wageningen, the Netherlands.

Received: 7 March 2023 Accepted: 31 May 2023