Defining indicators of resilience:
This study aimed to leverage the concept of "variability" in feed consumption as an indicator of animals' resilience. We hypothesized that animals exhibiting higher day-to-day variability in feed consumption might be more susceptible to external disturbances, including environmental and social stressors or disease, which could subsequently influence their feeding patterns. To obtain a comprehensive understanding of resilience, we investigated multiple indicators of this "variability".
Figure 1 illustrates the rationale adopted to derive the various resilience indicators in this study. In order to facilitate understanding of this process, we chose to present two distinct animals: one with higher resilience, represented by lower LnVar values (Panel A), and another with lower resilience and higher LnVar (Panel A). The first column of Fig. 1 displays the observed values, depicted by the green line, which represents a moving median calculated over a five-day period of feed consumption data (FCD). As a predictive model, linear regression was employed, and the resulting prediction is reported by the red line. Subsequently, we derived the Lag1 and LnVar indicators from the residual values (i.e. predicted values subtracted from the observed values) obtained through the regression analysis.
In the middle column of each panel (A and B) of Fig. 1, we illustrate the occurrence of ‘negative periods’ in red, representing instances where there were more than two consecutive days of negative residual values, indicating when the predictor exceeded the observed values. The largest area of negative period was identified and designated as the third resilience indicator, referred to as MaxArea.
Finally, the last column of each panel highlights the local minima (in red) and local maxima (in yellow). We utilized the sum of the local minima as the last resilience indicator, labeled as SumMin.
Figure 1
Figure 2 illustrates the distribution of the four resilience indicators. The normality assumption of these phenotypes is paramount for the subsequent analysis. In all cases normality assumption was not violated. To investigate potential distributional differences among the breeds considered, we conducted a post-hoc test using the Tukey method within a linear mixed model framework. The model incorporated effect of breed of the animal, as well as the nested effect of pen within room and the random effect of family (sire) of the animal. The results indicated that there were no significant differences observed among the breeds (P > 0.05). For more detailed information, please refer to Supplementary Material 1.
Figure 2
We performed a Spearman correlation analysis to examine the association of the resilience indicators with various performance traits of the pigs. These phenotypes included the percentage of muscle, backfat, and intramuscular fat, which were obtained at the end of the performance test, as well as the body weight measures at the end of the performance test, at an age of 156 ± 2.68 days. Additionally, we incorporated two traits measured during the period in which the resilience indicators were calculated: average daily feed consumption and average daily body weight.
The Spearman correlation values estimated using all individuals are presented in Fig. 3 (ALL), while the correlations estimates subsetting individuals from each breed Duroc (DR),Landrace (LR), Large White (LW) are depicted at the bottom of Fig. 3 (LW, LR, DR).
From the correlation between resilience indicators, a few clear patterns emerged. The Lag1, LnVar, and MaxArea exhibited positive correlations with each other, while SumMin revealed negative correlations with all other indicators, ranging from − 0.32 to -0.39. The strongest correlation was observed between LnVar and MaxArea, with a coefficient of 0.79. The Lag1 showed a correlation of 0.39 with LnVar and a stronger correlation of 0.60 with MaxArea.
When examining the correlations between resilience indicators and other productive traits, we discovered a negative correlation between the three positively correlated indicators— Lag1, LnVar, and MaxArea — and the percentage of muscle, with an average correlation coefficient of -0.15. The strongest correlation was observed between MaxArea and the percentage of muscle, with a coefficient of -0.21. Our analysis did not reveal any significant associations between the resilience indicators and the percentage of backfat, intramuscular fat, or final body weight.
When examining the correlation between resilience indicators and traits measured during the same time window, we observed that only LnVar showed a positive correlation with average FCD and average weight. On the other hand, despite the high correlation between MaxArea and LnVar, MaxArea exhibited an almost null correlation with FCD.
The correlation results suggest that resilience indicators are traits that show low or null collinearity with traditionally recorded traits in swine production. Correlation patterns within breeds remained consistent, as observed in Fig. 3 (LW, LR, DR). While the numeric values of the correlations may have varied slightly, the overall patterns remained unchanged.
Figure 3
Association between microbial composition and resilience
To investigate the potential association between the microbiota and resilience, we conducted a PERMANOVA analysis to determine whether variations in the composition of the microbiota were linked to the four resilience indicators. To account for possible confounding effects, we included the effects of Room and Breed in the analysis, along the resilience indicators.
The results of the PERMANOVA model are presented in Table 1. We found that all resilience indicators exhibited a significant association with microbial composition. Specifically, LnVar and MaxArea had the most pronounced influence on the microbiota composition (P < 0.001), followed by SumMin (P = 0.003), and Lag1 (P = 0.023). It is noteworthy that, although not the primary focus of this study, we also identified a significant impact of Room and Breeds (P < 0.001), which is further detailed in Supplementary Material 2.
Table 1
PERMANOVA analysis for the effect each of the four resilience phenotypes on microbial composition. The four phenotype was lag of one day of residual (Lag1), natural logarithm of residual variance (LnVar), area under the curve for periods with the largest consecutive negative errors (MaxArea), and sum of residual's local minima (SumMin).
Indicator | SumSq | F | Pr(> F) |
Lag1 | 1595 | 1.698 | 0.023 |
LnVar | 3208 | 3.416 | < 0.001 |
MaxArea | 3408 | 3.629 | < 0.001 |
SumMin | 2550 | 2.715 | 0.003 |
Indicator: Resilience phenotype; SumSq the total variation between the group means and the overall mean, F test Pr(> F) p value of the F statistic |
Table 1
We investigated the association between resilience and microbiota richness by performing a regression analysis, where we used the α-diversity of bacterial communities as the predictor variable and the resilience measures as the outcome variables. To quantify microbiota richness, we employed two diversity indices: the inverse of Simpson's index and the Shannon index. These indices capture both the richness (number of unique bacterial species) and evenness (distribution of species abundances) of the microbiota.
A significant (P < 0.01) and consistent decrease in both α-diversity measurements was observed as LnVar and MaxArea values increased. This suggests that animals with higher variability, indicating lower resilience, exhibited reduced microbiota richness. Similarly, negative trends were observed for Lag1, although statistical significance was not reached for the inverse Simpson index (P = 0.32). However, when richness was expressed using the Shannon diversity measure, the negative trend for Lag1 was slightly above the significance threshold (P = 0.07). These findings further support the notion that animals with higher variability and lower resilience tend to exhibit lower microbiota richness.
Conversely, SumMin showed almost negligible trends, suggesting that as local minimum values increased, α-diversity remained relatively stable. Although the trends between Lag1 and microbial composition were not statistically significant, a slight decline was observed when the Shannon index was used as the richness measure.
Figure 4
Consequently, we aimed to identify through differential abundance (DA) analysis which ASVs and their corresponding KEGG pathways were associated with these changes in the microbiome. By identifying these ASVs, we can gain insights into the specific microbial components that contribute to the observed changes in the microbiome associated with resilience.
Figure 5 presents the results of the ASV differential analysis. In the model, each ASV was treated as an independent factor. However, for clarity, the results in Fig. 5 are aggregated at the genus level, grouping together significant ASVs belonging to the same genus. In instances where multiple ASVs were associated with the same genus, the mean Log Fold Change (LFC) of the ASVs within each genus was reported.
Among the resilience indicators, LnVar and MaxArea exhibited notable differences in ASV abundance. Specifically, LnVar showed significant variations in the abundance of 7 ASVs across 4 different genera. On the other hand, MaxArea displayed significant differences in the abundance of 10 ASVs belonging to 7 different genera. These ASVs were associated with six distinct families: Lactobacillaceae, Ruminococcaceae, Prevotellaceae, Veillonellaceae, Lachnospiraceae, and Streptococcaceae. It is important to note that the last two families were only observed in relation to the MaxArea indicator.
Regarding the genera associated with resilience indicators, Lactobacillus exhibited a positive LFC value of 0.46 for LnVar and 0.51 for MaxArea, indicating its higher abundance in animals with greater variability. Lactobacillus was represented by three different ASVs in LnVar and by four ASVs (including the three from LnVar) in MaxArea. Rumminococcacea_UCG-014 displayed the most negative LFC values for both indicators, with − 0.35 for LnVar and − 0.6 for MaxArea. This genus was represented by two distinct ASVs in LnVar. Then, Misuokella and Prevotella_9 exhibited LFC values of -0.34 and − 0.28 for LnVar, and − 0.38 and − 0.33 for MaxArea, respectively. In addition, for MaxArea, ASVs belonging to the genera Streptococcus, Rumminococcacea_UCG-003, and Oscillospira were identified, with LFC values ranging from approximately − 0.29 to -0.27.
Figure 5
In addition to the analysis of individual ASVs, we also examined the abundances of KEGG pathways associated with enzymatic activities using a similar approach with the ANACOM-BC procedure. To maintain clarity and conciseness, in Fig. 7 we present the top 10 most influential microbiological enzymatic activities based on the absolute values of LFC. Supplementary Material 3 includes all KEGG pathways that surpassed the threshold for corrected P-values (Bonferroni threshold of 0.05).
Specifically, we identified 65 different classes of enzymatic activities for LnVar and 44 for MaxArea, indicating the potential involvement of a diverse range of microbiological enzymatic functions in relation to the resilience indicators.
In both indicators, LnVar and MaxArea, the most significant enzymatic activities were associated with the malonate cycle. Specifically, Malonate decarboxylase holo-[acyl-carrier protein] synthase (EC:2.7.7.66 as well as Malonyl-S-ACP decarboxylase (EC:4.1.1.87) and Acetyl-S-ACP:malonate ACP transferase (EC:2.3.1.187), were identified. These enzymatic activities involved in the malonate cycle contribute primarily to lipid metabolism.
NADH peroxidase (EC:1.11.1.1) exhibited higher abundance in less resilient animals, with an LFC value of approximately 0.8 for MaxArea and 0.3 for LnVar. This enzymatic activity along with Aryl-alcohol dehydrogenase (EC:1.1.1.90) belonged to the class of oxide reductase.
Sugar phosphatase (EC:3.1.3.23) involved in various metabolic pathways, showed higher abundance in less resilient animals, followed by UDP-N-acetylglucosamine kinase (EC:2.7.1.176).
Enzymatic activities more abundant in more resilient animals included those involved in the pathways of EC 5.2.1.1 - maleate isomerase and cycle, as well as 1-aminocyclopropane-1-carboxylate deaminase (EC:3.5.99.7), both associated with Cysteine and methionine metabolism. These enzymatic activities exhibited LFC values of approximately − 0.4 for MaxArea and − 0.25 for LnVar.
Finally, Sulfoacetaldehyde dehydrogenase (acylating) (EC:1.2.1.81) showed a positive LFC and was exclusively present in LnVar, whereas Sulfopropanediol 3-dehydrogenase (EC:1.1.1.308) was observed only when considering MaxArea as a resilience indicator.
Figure 6
In the subsequent phase of our investigation, we extended those analyses (i.e., α-diversity, PERMANOVA and DA) by categorizing the resilience measures into three distinct classes: highly-resilient (HR), lowly-resilient (LR), and a middle group (MR) acting as the control. This categorization was performed to amplify the differences observed at the extremes of the distribution of resilience phenotypes. By binning the resilience measures into these classes, we aimed to enhance the resolution and capture more pronounced distinctions between animals with high and low resilience.
We defined the HR animals as those exhibiting phenotypic values below the 95th percentile of the remaining animals within their respective breeds. Conversely, the LR animals were defined as those with values exceeding the 95th percentile threshold. Animals falling within the range of 47.5–52.5% were classified as the MR, which served as the control group in our analyses.
Consistent with our previous PERMANOVA analysis, we observed significant variations in microbiota composition among the different classes derived from all resilience indicators (refer to Supplementary Material 3). Interestingly, when comparing the MR to either the LR and HR using PERMANOVA, we found that there were no significant differences when transitioning from the MR to the HR for all resilience indicators. However, when transitioning from the MR to the LR group, significant differences in microbiota composition were observed for Lag1, LnVar, and MaxArea (P < 0.001). This suggests that the changes in microbiota composition may be more pronounced in animals with lower resilience compared to the more resilient ones, and the relationship between resilience and microbial composition might be non-linear.
The analysis of α-diversity further confirms the observed trends. Figure 7 show that significant differences in α-diversity (P < 0.001) were observed among resilience classes. Notably, these differences were observed solely when transitioning from MR or HR classes to the LR class. To simplify the presentation, only the Shannon index for α-diversity was reported in Fig. 7, while the results for the inverse Simpson index can be found in Supplementary Material 5.
These findings align with the results obtained when resilience was considered as a continuous covariate, indicating that a decline in α-diversity is associated with lower resilience in animals.
Figure 7
The DA analysis was also performed using class-based expression of different resilience indicators instead of linear variables. The results overlapped with the previous findings, both for individual ASVs and KEGG pathways, see Supplementary Material 6. Interestingly, the observed differential abundance of ASVs was primarily driven by comparisons between the less resilient animals and the control groups, rather than comparisons between the more resilient animals and the control group. For example, in Fig. 9 is possible to observed that that no significant ASVs were identified when comparing the HR and MR groups.
We also DA to investigate the abundance differences of ASVs among different classes of the four resilience indicators. The results were consistent with the previous analysis when resilience was considered as a linear variable. These analyses provided additional insights: the observed differential abundance of ASVs was primarily driven by comparisons between the less resilient animals and the control groups, rather than comparisons between the more resilient animals and the control group. In fact, no significant ASVs were identified when comparing the HR and MR groups
Figure 8
After identifying the connection and understanding the biological basis between gut microbiome composition and resilience, our next objective was to quantify the extent to which the microbiome contributes to the regulation of resilience. We aimed to determine the proportion of phenotypic variance in resilience that can be attributed to the microbial composition, which is commonly refer to as "microability" (m2) [28].
To address the influence of both environmental conditions and genetic background on resilience, we incorporated them into our analysis. Environmental conditions were accounted for by including the combination of the room and pen as a factor in the model. Additionally, the genetic background was considered by incorporating a sire effect in the model. While our primary focus was to quantify the impact of microbiota composition on resilience, we also provide a brief overview of the influence of these other factors on resilience.
The percentage contribution of each factor to the overall variance in the analysis is depicted in Fig. 9. Microbiota compositions explained a similar amount of variance for all resilience indicators, ranging from 10% for Lag1 and MaxArea to 11% for LnVar and SumMin. The highest posterior density intervals (HPD95) of m2 were reported in Part B of Fig. 4. For the first three indicators, the HPD95 intervals ranged from 0.04 to 0.16, while a wider HPD95 interval of 0.04 to 0.19 was observed for SumMin. Across all resilience indicators, the genetic background of the animals accounted for approximately 5–6% of the phenotypic variance. However, the influence of the environment varied depending on the specific indicator. For Lag1 and SumMin, the environment explained approximately 14–20% of the total phenotypic variance. In the case of MaxArea and LnVar, the environment contributed significantly more, accounting for approximately half of the total phenotypic variance. Specifically, the environment explained 37% of the phenotypic variance for MaxArea and 47% for LnVar.
Figure 9
Finally, we evaluated the discriminative capability of the microbiological data across the different classes of resilience indicators, we conducted a Partial Least Squares-Discriminant Analysis (PLS-DA). This analysis is of practical significance as it allows us to utilize microbiome composition as a potential biomarker for distinguishing and predicting animals based on their resilience levels.
Overall, our findings align with previous observations. The PLS-DA analysis demonstrated significant P-values (P < 0.05) for all resilience indicators when assessing the ability to discriminate less resilient animals. However, there were no significant differences in the discriminative ability when comparing the medium or high resilience classes.
The results of the analysis, illustrated in Fig. 10(A), revealed distinct clustering patterns for the HR and LR resilience classes in relation to the LnVar and MaxArea indicators. This suggests a clear differentiation between animals with high and low resilience levels. However, for SumMin, the separation between the classes was less pronounced, and no clear differentiation was observed for SumMin. Furthermore, the control class MR did not exhibit complete separation and appeared to be positioned between the High and Low classes in terms of LnVar and MaxArea. However, a clearer distinction was observed between the control class and the High and Low classes in the case of Lag1.
Additionally, individuals of different breeds were equally distributed among the different resilience classes, indicating the absence of breed-specific patterns of resilience.
The Receiver Operating Characteristic (ROC) plot (Fig. 10(B)) reinforced our previous observations, indicating that microbiome composition is more effective in discriminating between resilience classes. This effect was particularly evident for MaxArea and SumMin. Moreover, the confusion matrix further supported (Fig. 10C) these findings, with a higher number of correctly identified instances in the low resilience group across all indicators. The values presented in Table 2 support the previous statement. With the exception of SumMin, the other resilience indicators demonstrated significant ability to discriminate low resilience animals. This was particularly evident for LnVar and MaxArea, which had P-values of < 0.001. Lag1 showed a P-value of 0.01, along with AUC values of 0.726, 0.826, and 0.891 for Lag1, LnVar, and MaxArea, respectively. However, the AUC for discriminating the low-class using SumMin was 0.555.
Table 2
Performance of the Partial Least Squares Discriminant Analysis (PLS-DA). The table provides the classification results for the four resilience phenotypes: lag of one day of residual (Lag1), natural logarithm of residual variance (LnVar), area under the curve for periods with the largest consecutive negative errors (MaxArea), and sum of residual's local minima (SumMin). The "Class" column represents the classification of animal groups within each indicator, indicating Low resilience (L), High resilience (H), and Medium resilience (M). The "AUC Class" column displays the Area Under the Curve (AUC) values for distinguishing each class of resilience. For each class, the table reports the number of positive events (animals belonging to that class) and negative events (animals not belonging to that class)."P-values" column indicates the significance level of the classification. The "Overall AUC" column represents the AUC performance in discriminating all classes combined.
Indicator | class | AUC class | yes events | no events | P-values | AUC all classes |
Lag1 | L | 0.726 | 13 | 27 | 0.010 | 0.693 |
M | 0.682 | 15 | 25 | 0.028 |
H | 0.672 | 12 | 28 | 0.044 |
LnVar | L | 0.826 | 13 | 27 | > 0.001 | 0.726 |
M | 0.629 | 15 | 25 | 0.091 |
H | 0.729 | 12 | 28 | 0.011 |
MaxArea | L | 0.891 | 13 | 27 | > 0.001 | 0.691 |
M | 0.637 | 15 | 25 | 0.078 |
H | 0.544 | 12 | 28 | 0.336 |
SumMin | L | 0.555 | 13 | 27 | 0.293 | 0.617 |
M | 0.658 | 15 | 25 | 0.050 |
H | 0.639 | 12 | 28 | 0.086 |
Interestingly, despite the P-values for classifying the MR and HR classes being close to significance at 0.05 and 0.08, respectively, SumMin was the only indicator in which the microbiome exhibited superior ability to discriminate between the MR and HR classes compared to the LR class. This finding suggests a negative correlation between SumMin and the other indicators. Consequently, higher values of HR resilience could potentially indicate the contrary i.e., lower overall resilience in this context.
Lag1 was the only indicator in which microbiome could significantly discriminated the MR and HR classes, while in LnVar microbiome significantly discriminated the HR class from the others with a P-value of 0.011. When considering the AUC for discriminating all categories, LnVar had the highest value of 0.726, followed by Lag1 and MaxArea with values of 0.693 and 0.691, respectively. Lag1 had the lowest AUC value of 0.617.
In summary, we can conclude that all indicators of microbiome composition were effective in discriminating low resilience animals from the rest of the population, except for SumMin, for which higher values represent higher animal resilience. Among these indicators, MaxArea supported the best performance, while LnVar exhibited the best overall classification performance in terms of microbiome composition.
Figure 10