Heritability and genetic correlations of plasma metabolites of pigs with production, resilience and carcass traits under natural polymicrobial disease challenge

Metabolites in plasma of healthy nursery pigs were quantified using nuclear magnetic resonance. Heritabilities of metabolite concentration were estimated along with their phenotypic and genetic correlations with performance, resilience, and carcass traits in growing pigs exposed to a natural polymicrobial disease challenge. Variance components were estimated by GBLUP. Heritability estimates were low to moderate (0.11 ± 0.08 to 0.19 ± 0.08) for 14 metabolites, moderate to high (0.22 ± 0.09 to 0.39 ± 0.08) for 17 metabolites, and highest for l-glutamic acid (0.41 ± 0.09) and hypoxanthine (0.42 ± 0.08). Phenotypic correlation estimates of plasma metabolites with performance and carcass traits were generally very low. Significant genetic correlation estimates with performance and carcass traits were found for several measures of growth and feed intake. Interestingly the plasma concentration of oxoglutarate was genetically negatively correlated with treatments received across the challenge nursery and finisher (− 0.49 ± 0.28; P < 0.05) and creatinine was positively correlated with mortality in the challenge nursery (0.85 ± 0.76; P < 0.05). These results suggest that some plasma metabolite phenotypes collected from healthy nursery pigs are moderately heritable and genetic correlations with measures of performance and resilience after disease challenge suggest they may be potential genetic indicators of disease resilience.

www.nature.com/scientificreports/ addition, metabolomics offers the potential to identify new phenotypes or traits that can be used for selection of resilient animals. There is evidence in the literature that blood metabolites in humans and cattle are heritable [14][15][16][17] . Given the similarities that exist between porcine and human, blood metabolites in porcine might be expected to be heritable. Heritability is a prerequisite for the use of a trait in breeding strategies and, ideally, indicators of resilience should be expressed in healthy animals, such that they can be measured in high-health nucleus farms to provide prediction of resilience in commercial environments 18 . It is therefore important to understand the genetic architecture of these potential new phenotypes.
This study is part of a larger project entitled "Phenomics for genetic and genome-enabled improvement of resilience in pigs", which had identification of predictors of resilience on young healthy pigs as one of its main objectives. The project uses a natural polymicrobial disease challenge model 18 with collection of resilience related traits, including average daily gain (ADG), feed intake and feed intake duration (ADFI and ADFD), number of individual health treatments (nTRT) and mortality, residual feed intake (RFI), and feed conversion ratio (FCR). In addition, new resilience traits based on day-to-day variation in feeding patterns were proposed by Putz et al 18 .
The overall purpose of the current study was to use data and plasma samples collected on healthy pigs prior to pathogen/disease challenge to estimate heritability of 44 metabolites and 2 amino acid indices, the phenotypic and genetic parameters of plasma metabolite concentration in relation to their subsequent performance, disease resilience, and carcass traits under the natural disease challenge.
Moreover, we estimated the phenotypic correlations between the metabolites that are involved in the same pathway (glycine, serine, alanine and threonine metabolism) and the results showed that betaine was positively correlated with dimethylglycine and l-serine, and l-serine was positively correlated with l-methionine (Supplementary Table S2). In addition, betaine did have significant positive phenotypic correlations with l-glycine (Supplementary Table S3; P < 0.05).
Genetic correlations. Genotypic correlations of the metabolites with performance, resilience, and carcass traits were larger than phenotypic correlation (Supplementary Figure S1b) however with larger SE. Estimates of genetic correlation between metabolites and ADG in the three phases, qNurADG, cNurADG and FinADG are provided in Table 2. Overall, estimates of genetic correlations among plasma metabolites and ADG in the three phases were very low, with high SE. However, some metabolites showed significant correlation estimates, with the largest negative genetic correlation between plasma creatinine and qNurADG (− 0.60 ± 0.18). Dimethylglycine, betaine, L-methionine and L-serine showed positive genetic correlation estimates with qNurADG. Metabolites that were positively correlated with qNurADG are visualized in a compound network in Fig. 1. The metabolites that were genetically positively correlated with qNurADG are involved in metabolic pathways, including the glycine and serine, methionine and cysteine, glycerophospholipid, and glycosphingolipid pathways. We estimated the genetic correlations between the metabolites that are involved in the same pathway (glycine, serine, alanine and threonine metabolism) and the results showed no significant genetic correlations between these metabolites (Supplementary Table S2).
No significant genetic correlations were estimated between metabolites and cNurADG. Only two metabolites namely: l-glutamic acid and oxoglutarate were estimated to have significant (positive) genetic correlations with FinADG. Table 3 shows estimates of genetic correlations of metabolites with the feed intake traits ADFI, ADFD, FCR, and RFI. l-Glutamic acid was positively correlated with ADFI (P < 0.001), while several other metabolites, including dimethylglycine and l-aspartate, tended to be positively genetically correlated with ADFI (0.05 < P < 0.1). Four metabolites that were positively correlated with ADFD (Table 3; P < 0.05): dimethylglycine, l-glycine, betaine and citric acid were visualized in a compound network (Fig. 2). These metabolites are involved in metabolic www.nature.com/scientificreports/ pathways such as the TCA cycle and the glycine and serine, methionine, and cysteine pathways but they were not significantly genetically correlated with each other (Supplementary Table S3). In addition, the l-glutamine was positively correlated at the genotypic level with FCR (P < 0.05) and isobutyric acid was negatively correlated with RFI (P < 0.05). Table 1. Estimates and standard errors (SE) of heritability and litter effects, phenotypic variance SD and variance due to pen (σ 2 P ). Branched-chain amino acids (BCAA) was calculated as the sum of l-leucine, l-isoleucine and l-valine and ketogenic amino acids (ketoAA) was calculated as the sum of l-lysine and l-leucine. Significance of the heritability is indicated as ***, **, *, corresponding to P < 0.001, P < 0.01 and P ≤ 0.05. "-" indicates not estimable.

Metabolite
Heritability (SE) Litter (SE) Phenotypic SD (SE) σ 2 P (SE) www.nature.com/scientificreports/ Table 2. Estimates (SE in parentheses) of genetic correlations between the most heritable metabolites and average daily gain in three stages, i.e. the quarantine nursery, the challenge nursery, and the finisher.
Significance of genetic correlations are indicated in bold as: **, *, x corresponding to P < 0.01, P ≤ 0.05 and 0.05 < P < 0.10 respectively; "-" indicates not estimable. www.nature.com/scientificreports/ Estimates of genetic correlations of metabolites with carcass traits are presented in Table 4. The only significant estimate was of citric acid with carcass back fat (CBF) (P < 0.05).
Finally, only oxoglutarate had a significant negative genetic correlation estimate with the number of treatments across the challenge nursery and finisher (Table 5; P < 0.05), while creatinine had a significant positive genetic correlation with mortality in the challenge nursery (Table 6; P < 0.05).

Discussion
The objectives of this study were to estimate (1) heritabilities of metabolites in plasma of young healthy pigs and (2) their genetic correlations with production, disease resilience, and carcass traits in growing pigs that were exposed to a natural polymicrobial disease challenge. The metabolomics samples were collected from healthy pigs at an average of 26 days of age in the quarantine nursery and the traits analyzed included those in the quarantine nursery, as well as in the challenge nursery and finisher, and carcass traits at slaughter. This study contributes to one of the main overarching objectives of the polymicrobial natural disease challenge model, which is to identify genetic predictors of resilience using samples from healthy pigs. Such predictors would be very useful for Table 3. Estimates of genetic correlations (SE in parentheses) between the most heritable metabolites and average daily feed intake (ADFI), duration (ADFD), feed conversion ratio (FCR) and residual feed intake (RFI). Significance of genetic correlations are indicated in bold as: ***, **, *, x corresponding to P < 0.001, P < 0.01, P ≤ 0.05 and 0.05 < P < 0.10 respectively. "-" indicates not estimable. www.nature.com/scientificreports/ breeding for disease resilience, as they can be measured in high health genetic nucleus herds as indicator traits for disease resilience in commercial farms. We measured the concentration of a panel of 44 metabolites in plasma, including amino acids, short chain fatty acids, sugars, alcohols, organic acids, amines, and TCA cycle intermediates, and urea cycle intermediates. Variation in metabolite concentration can be due to environmental effects, diet, gender, age, physiological conditions, and genetic effects. Literature indicates that in humans, approximately 50% of phenotypic differences Table 4. Estimates of genetic correlations (SE in parentheses) between the most heritable metabolites and carcass traits: carcass weight (CWT), backfat (CBF) and loin depth (CLD), lean yield (LYLD) and dressing percentage (DRS). Significance of genetic correlations are indicated in bold as: *, x corresponding to P ≤ 0.05 and 0.05 < P < 0.10 respectively. "-" indicates not estimable.  www.nature.com/scientificreports/ in metabolite levels is due to genetics, but heritability estimates differ across metabolite classes 16 . Metabolites can be grouped into primary and secondary metabolites. Primary metabolites are directly involved in primary metabolic processes, such as normal growth, development, reproduction, and immune response, e.g. amino acids and products derived from glycolysis and the TCA cycle. Primary metabolites are highly conserved across species and serve as precursors for the synthesis of secondary metabolites 19 . For example, amino acids, in addition to being the building blocks of proteins, are also regulators of innate and adaptive immune responses in living cells. Many studies have demonstrated that glutamic acid, glutamine, histidine, methionine, leucine, isoleucine, and valine are functional regulators of macrophages, dendritic cells, and T-cells. [20][21][22][23] . Our results showed that 31 metabolites had low to moderate heritability and two metabolites had relatively high heritability (> 0.4). In a study of beef cattle, only 11 of 33 metabolites measured (29 in common with this study) were reported to be heritable 17 . Similar heritability estimates to Li et al. 17 were reported here for betaine, creatinine, pyruvic acid and citric acid. In our study, other metabolites, such as ketone bodies (3-hydroxybutyric acid, acetoacetate and acetone), creatine, succinate, formate, and methylhistidine, showed negligible heritability estimates, suggesting that they are primarily influenced and manipulated by environmental effects such as diet, and/or age, health status etc. Branched-chain amino acids, was calculated as the sum of l-leucine, l-isoleucine and l-valine and BCAA did not have higher heritability than the amino acids that contributed to the index.
Phenotypic correlations of metabolites with production, disease resilience, and carcass traits were generally very low. These results are in line with those reported previously in dairy cows by Buitenhuis et al. 15 . However, genetic correlation analyses found that dimethylglycine (0.28), betaine (0.39), l-methionine (0.45), and l-serine (0.54) were positively correlated with ADG in the quarantine nursery, while creatinine (− 0.60) was negatively correlated with this trait. These metabolites are involved in glycine, serine, alanine and threonine pathway, which suggests that this pathway might be a target to improve ADG in young healthy piglets. Indeed, supplementation with dimethylglycine has been shown to improve growth performance, significantly increasing total body weight gain and feed intake, and improving feed efficiency in low-birth-weight piglets 24 . We found that l-methionine (0.15), l-glutamine (0.12) and betaine (0.18), had positive phenotypic correlation estimates with qNurADG. Methionine supplementation has been reported to improve growth rate in nursery pigs 25 . Positive effects of methionine supplementation on intestine structure have also been observed, with greater average daily gain from days 7-14 of age and improved feed efficiency 25 .
Creatinine was genetically highly negatively correlated (− 0.60) with ADG in the quarantine nursery. Variation in creatinine concentration can arise as result of environmental and genetic factors. For example, phenotypic variation in creatinine concentration can be expected due to batch, transportation and fighting or re-grouping of animals. Creatinine is considered a waste product produced by muscles from the breakdown of creatine and is removed from the blood and released into urine by the kidneys. Typically, 95% of creatinine is found in muscle and an increase of creatinine concentration in blood is often considered as an indicator of kidney malfunction 27 . Conversion of muscle creatine into creatinine can reflect protein degradation 28 and serum creatinine concentration has been proposed as an indicator of protein deposition 29 . Animals with low ADG are expected to have lower protein or muscle deposition and our results suggest that increased plasma creatinine concentration might be a Table 6. Estimates of genetic correlations (SE in parentheses) between the most heritable metabolites and mortality in the challenge nursery, finisher, and across the nursery and finisher. Significance of genetic correlations are indicated in bold as: *, corresponding to P ≤ 0.05. "-" indicates not estimable.

Metabolite
Challenge www.nature.com/scientificreports/ genetic indicator of low ADG in healthy nursery piglets. In addition, creatinine was phenotypically negatively correlated with ADG of young healthy pigs (− 0.21). Moreover, creatinine was genetically positively correlated with mortality in the challenge nursery, suggesting that young healthy pigs that have higher plasma creatinine, genetically have lower ADG as healthy nursery pigs but are more likely to die when challenged by disease. Further research is necessary to validate creatinine as a potential genetic marker for ADG in healthy nursery pigs and as an early genetic indicator trait for mortality under disease. Plasma oxoglutarate concentration in the quarantine nursery was genetically positively correlated with ADG in the finisher and negatively correlated with nTRT across the nursery and finisher. Moreover, l-glutamic acid showed positive genetic correlation estimates with ADG (0.72) and with ADFI (0.62) in the finisher. Oxoglutarate, also known as alpha-ketoglutarate, is a key organic acid of the TCA cycle and a source of glutamic acid and glutamine, and stimulates protein synthesis and inhibits protein degradation in muscles 30,31 . Positive effects of oxoglutarate on the protein synthesis and skeletal system have been reported in various farmed species, including turkeys 32 , pigs 33,34 , and sheep 35,36 . Our results suggest that young healthy pigs that had greater oxoglutarate and l-glutamic acid concentration, genetically have greater ADG in the finisher and received fewer health treatments across the challenge nursery and finisher. Interestingly, none of the metabolites were genetically correlated with ADG in the challenge nursery. For the data used here, Cheng et al. 37 reported that the estimate of heritability of ADG was lower in the challenge nursery (0.19) than in the quarantine nursery (0.31) and in the finisher (0.30), which might have impacted the ability to accurately estimate the genetic correlations of metabolites with ADG in the challenge nursery.
Four metabolites (dimethylglycine, l-glycine, betaine and citric acid) had moderate to high positive genetic correlation estimates with ADFD (0.30-0.52). These metabolites are involved in the TCA cycle and in glycine, serine, alanine, and threonine metabolism. Interestingly, dimethylglycine and betaine were also positively correlated with ADG in the quarantine nursery, which suggests that pigs with greater dimethylglycine and betaine concentration genetically have greater ADG as young healthy pigs and might spend more time eating in the finisher stage.
A moderate positive genetic correlation (0.34) between citric acid and carcass backfat was estimated, which is in line with the fact that citric acid is involved in fatty acid metabolism. Citrate, the conjugated base of citric acid, is formed from oxaloacetate and acetyl-coenzyme A (acetyl-CoA) by citrate synthase. Once transported to the cytosol, citrate is converted to acetyl CoA and oxaloacetate. Acetyl-CoA is then converted to malonyl CoA and can be used as a substrate for fatty acid synthesis 38 . The genetic correlation observed suggests that young healthy pigs that have higher plasma content of citric acid genetically have higher carcass backfat when disease is present.
Interestingly, isobutyric acid was estimated to be negatively correlated, genetically, with RFI (− 0.38), which suggests that young healthy pigs that have higher plasma isobutyric acid content genetically have better feed efficiency under disease conditions. Isobutyric acid is a carboxylic or short chain fatty acid (SCFA) that is generated via microbial (gut) metabolism. Isobutyric acid has been described and quantified in the faeces of human, rats, horses, and pigs 39 . It has been reported that SCFAs such as acetic, propionic, and butyric acids, derive mostly from carbohydrates, while other SCFAs such as isobutyric acid are mostly derived from proteins 40 , specifically from degradation of branched-chain amino acids, and greater concentration of isobutyric acid could be indicative of better utilization of dietary protein by the microbiota 41 . Literature indicates that low RFI animals (efficient) have greater ileal isobutyric acid concentrations 42 . Thus, our results suggest that low RFI pigs genetically might have a better utilization of dietary protein by the microbiota. Further research is necessary to investigate the role of isobutyric acid as a potential genetic biomarker of RFI in pigs with or without disease.
In conclusion, the results suggest that some metabolic phenotypes measured in plasma of young healthy pigs are moderately heritable. The present work contributes to our understanding of the genetic parameters of plasma metabolite concentrations in young healthy pigs and their relationships with production, resilience, and carcass traits under disease. To the best of our knowledge, this is the first study to report estimates of heritabilities of plasma metabolite concentrations and of their genetic correlations with production, resilience, and carcass traits in pigs following polymicrobial disease challenge. Metabolites have the potential to be used in high health genetic nucleus herds as indicator traits for disease resilience in commercial farms. Further studies are warranted to validate the identified possible genetic indicators of resilience.

Material and methods
Animals, production, resilience and slaughter traits. This study was carried out in accordance with the Canadian Council on Animal Care guidelines (CCAC) 43 . All procedures were carried in accordance with the Animal Research: Reporting of In Vivo Experiments (ARRIVE) guidelines 44 . The Animal Care protocol was approved by the Animal Protection Committee of the Centre de Recherche en Sciences Animales de Deschambault (15PO283) and Animal Care and Use Committee at the University of Alberta (AUP00002227).
This study is part of a larger research project which investigates the underlying genetic mechanisms of disease resilience in grow-finisher pigs exposed to a natural polymicrobial disease challenge. The details of the polymicrobial challenge and phenotypes/traits collected were previously described in Putz et al. 18 , Cheng et al. 37 , and Bai et al. 45 . The challenge was established in late 2015 at the Centre de développement du porc du Québec inc (CDPQ) test station in Québec, Canada, with the aim to mimic a commercial farm with high disease pressure to maximize expression of genetic differences in resilience. For purpose of clarity, here we briefly describe the three phases of the experiment, which included a pre-challenge quarantine nursery (19 days on average beginning at 3 weeks of age), the challenge nursery (27 days on average), and lastly, the finishing phase (100 days on average). The natural disease challenge protocol was established by introducing naturally infected animals with known diseases into the challenge nursery and finisher barn at CDPQ. Some of the introduced pathogens included: viruses (PRRSV and swine influenza virus A), bacterial pathogens such as Brachyspira hampsonii, Haemophilus www.nature.com/scientificreports/ parasuis, Mycoplasma hyopneumoniae, Salmonella enterica serovar Typhimurium, and Streptococcus suis), and two parasites (Ascaris suum and Cystoisospora suis) 45 . Moreover, environmental enrichment (inedible point source objects) was applied in 50% of quarantine nursery pens in some of the batches, with the purpose to evaluate the impact of environmental enrichment on disease resilience. Pigs that received enrichment in the quarantine nursery, continued to receive enrichment in the challenge nursery and finisher 37 . Throughout the study feed was available ad libitum (the quarantine nursery, challenging nursery and finisher). In the quarantine nursery all pigs were fed the same commercial diet appropriate for pigs' age and weight (Délice, Nourisson and Premier Age (Cie Alfred Couture ltée; Quebec, Canada)). Production, resilience, and carcass traits were collected from a total of 3205 F1 crossbred (Landrace × Yorkshire) barrows. The phenotypes used in this study included: average daily gain recorded in the quarantine nursery, challenge nursery, and finisher stages. Finisher traits considered in the present study included feed intake and duration, feed conversion ratio, and residual feed intake. The number of parenteral treatments provided to individual pigs were also tabulated separately for the challenge nursery, finisher and combined challenge nursery and finisher as the number of treatments per 180 days. Mortality (0 = survived, 1 = died) in the challenge nursery, finisher, and across the nursery and finisher were also included 37 . Finally, carcass traits included: carcass weight, backfat, and loin depth, dressing percentage, and lean yield. Details of the recording and derivation of the resilience, production, and carcass traits can be found in Putz et al. 18 and Cheng et al. 37 . Cheng et al. 37 used two data sets for analysis of traits recorded in the finisher but in the current study we only use the data set that included phenotypic finisher data on pigs that survived to slaughter (survivor data), except for mortality. For the challenge nursery, phenotypes on all pigs were used. Details of the number of animals and traits considered here are described in Cheng et al. 37 . For purpose of brevity: data from 958 young healthy pigs were included in the present study for metabolomics analysis, nearly 3200 animals for nursery traits, around 2500 animals for finisher traits, and 2000 animals for carcass traits.
Blood samples. For metabolomics analysis we used blood samples collected from 958 young healthy animals which were introduced in 15 batches of 60 or 75 pigs. Blood was collected from the jugular vein into K2 ethylenediaminetetraacetic acid (EDTA) tubes (BD Vacutainer, Blood Collection Tubes, United States), on all pigs in the quarantine nursery at an average age of 26 days, 5 days post-arrival from their farm of origin 45 . After collection, the blood samples were centrifuged at 3000 rpm at 4 °C for 10 min, plasma collected and immediately frozen and stored at − 80 °C and only thawed for the metabolomics analysis. Two weeks after the first sampling, pigs were transferred to the test station and naturally exposed to multiple pathogens, as described in Putz et al. 18 and Bai et al. 45 .
All animals (n = 3205) were genotyped using a 650k Affymetrix Axiom Porcine Genotyping Array by Delta Genomics (Edmonton AB, Canada). Raw Affymetrix SNP data were processed by Delta Genomics, separately for each cycle, with the Axiom Analysis Suite, using all defaults. Details of genotyping and quality control are described in Cheng et al. 37 and Bai et al. 45 . After quality control, a total of 417,443 SNPs in 3205 pigs remained and were used for analysis.
Nuclear magnetic resonance spectroscopy and quality control. Forty-eight plasma metabolites were quantified using NMR, following established protocols at The Metabolomics Innovation Center at University of Alberta (TMIC), AB, Canada (https:// www. metab olomi cscen tre. ca/). Plasma samples were thawed on ice and a deproteinization step, involving ultra-filtration was performed as previously described 46 , in order to remove plasma macromolecules. Prior to filtration, 3 kDa cut-off centrifugal filter units (Amicon Microcon YM-3), were rinsed five times each with 0.5 mL of H 2 O and centrifuged (10,000 rpm for 10 min) to remove residual glycerol bound to the filter membranes. Aliquots of each plasma sample were then transferred into the centrifuge filter devices and spun (10,000 rpm for 20 min) to remove macromolecules (primarily protein and lipoproteins) from the sample. The filtrates were collected and the volumes for each sample were recorded. If the total volume of the sample was under 250 µL an appropriate amount of 150 mM KH 2 PO 4 buffer (pH 7) was added and the dilution factor was annotated and taken into account in the analysis. Subsequently, 46.5 µL of a standard buffer solution (54% D 2 O:46% 1.75 mM KH 2 PO 4 pH 7.0 v/v containing 5.84 mM DSS (2,2-dimethyl-2-silcepentane-5-sulphonate), 5.84 mM 2-chloropyrimidine-5 carboxylate, and 0.1% NaN 3 in H 2 O) was added to the sample.
The plasma samples (250 µL) were transferred in 3 mm SampleJet NMR tubes for spectral analysis. All 1 H-NMR spectra were collected on a 700 MHz Avance III (Bruker) spectrometer equipped with a 5 mm HCN Z-gradient pulsed-field gradient (PFG) cryoprobe. 1 H-NMR spectra were acquired at 25 °C using the first transient of the NOESY pre-saturation pulse sequence (noesy1dpr), chosen for its high degree of quantitative accuracy 47 . All FID's (free induction decays) were zero-filled to 250 K data points. The singlet produced by the DSS methyl groups was used as an internal standard for chemical shift referencing (set to 0 ppm) and for quantification all 1 H-NMR spectra were processed and analyzed using an in-house version of the MAGMET automated analysis software package using a custom metabolite library. MAGMET allows for qualitative and quantitative analysis of an NMR spectrum by automatically fitting spectral signatures from an internal database to the spectrum 48 . Each spectrum was inspected by an NMR spectroscopist in order to minimize compound misidentification and misquantification.
Prior to statistical analysis, a quality control step was applied to the metabolite data. Four metabolites that were frequently (> 20%) below the limit of detection or with at least 20% missing values were removed from consideration. A total of 44 metabolites and two amino acid indexes remained in the dataset. Other missing values (15 data points) were replaced by the median value of each metabolite in the original data. First, we assessed the significance of each fixed (batch , enrichment), covariable (age), and random effects (pen and litter) using linear www.nature.com/scientificreports/ regression models implemented in R statistical software 49 . The residuals of the model were plotted and visually inspected for the presence of outliers, which were excluded from the dataset. Data normalization (log10) of metabolite concentrations that were not normally distributed (2-hydroxybutyrate, ethanol, 3-hydroxybutyric acid, l-alpha-aminobutyric acid, methanol and creatine) was done prior to statistical analysis. Two indexes were also computed for statistical analysis: (1) branched-chain amino acids (BCAA), which was calculated as the sum of l-leucine, l-isoleucine and l-valine and (2) ketogenic amino acids (ketoAA), calculated as the sum of l-lysine and l-leucine. A summary of descriptive statistics of all metabolites after quality control is in Supplementary Table S4.
Variance component analyses. Variance components were estimated by GBLUP using the BLUPF90 programs 50 . The general following mixed linear model was used to estimate the heritability of each metabolite: where Y ijk is the trait (metabolite); Batch i is the fixed batch effect (i = 1, …, 15); Age ijk is the covariate of age when the pig entered the quarantine nursery; Pen j is the random effect of pen by batch corresponding the different phases (quarantine nursery, challenging nursery, or finisher), with Pen j ~ N (0, σ 2 P ) where σ 2 P is pen variance; Litter ijk is the litter environmental effect, with Litter ijk ~ N (0, σ 2 L ) where σ 2 L is the litter environmental variance; u ijk is the random additive genetic effect, with the vector u ~ N (0, Gσ 2 A ), where G is the genomic relationship matrix and σ 2 A is the additive genetic variance; and e ijk is the residual effect, with e ijk ~ N (0, σ 2 e ) where σ 2 e is the residual variance. The genomic relationship matrix, G, was created separately for each of the seven companies supplying pigs using the software preGSf90 50 and the method described by VanRaden 51 , and then combined into one G matrix, with genetic relationships between companies set to zero in order to focus on within-company variance components, as described by Cheng et al. 37 . For six metabolites namely: l-ornithine, l-leucine, l-valine, l-asparagine, 3-methyl 2-oxovaleric acid and formate, environment enrichment was included as a fixed effect because the effect of enrichment was significant (P ≤ 0.05).
For metabolites with heritability estimates greater than 0.20, genetic correlations with production, disease resilience, and carcass traits were estimated using bivariate models. For these traits, phenotypes from batches that were not measured for metabolites were also included. Models used for these traits were as described by Cheng et al. 37 .
Heritability was estimated as σ 2 a / (σ 2 a + σ 2 L + σ 2 e ) and the standard error (SE) for the heritability estimate was calculated according to the Monte-Carlo method suggested by Meyer and Houle 52 . The proportion of variance explained by sow or maternal effects, referring to effects that are common to individuals with the same mother, was estimated as σ 2 L / (σ 2 L + σ 2 e ). Genetic correlations between two traits were estimated as the estimate of the genetic covariance from the bivariate analysis divided by the product of the genetic standard deviations for the two traits. Significance of estimates of heritabilities and genetic correlations were determined using likelihood ratio tests with 1 degree of freedom. For heritability estimates, the resulting P-values were divided by 2 because the estimates are restricted to be positive 53 .

Network visualization.
For exploration and visualization of the biochemical pathway that metabolites are involved in, the Metscape plugin 54 in Cytoscape 3.8.2 55 was used. The metabolites that had significant (P ≤ 0.05) genetic correlations with qNurADG and ADFD were used for pathway visualization. The file with the list of KEGG elements was loaded into Metscape to generate a compound network. Only metabolites with KEGG IDs were considered for compound network and pathway analysis. In a compound network, metabolites are represented as nodes and reactions are represented as edges. A compound node with an outgoing edge is a substrate, while a compound with an incoming edge is the product of a specific biochemical reaction. Finally, we estimated phenotypic and genetic correlations among metabolites that belong to the same pathway.

Data availability
The data analyzed in this study are not publicly available, because the data were generated on samples from commercially owned animals, but they can be made available by the corresponding author on reasonable request.