**2. Results**

#### *2.1. Descriptive Statistics and Linear Model Analysis for genes and Metabolites*

The data on 749 metabolites and 25,880 genes from 40 samples were generated using untargeted metabolomics and transcriptomic approaches, respectively. We utilized data of 405 annotated metabolites (see methodology) for further analysis. For the transcriptomic data generated on 25,880 genes, we analyzed the data for each of the two groups (breed-specific and FE-specific), as described in the methodology. The genes with a gene count<1 were removed, resulting in 20,233 genes for both the groups. The gene count data for each group (breed-specific and FE-specific) was normalized, and the linear model was fitted into the data as given in the methodology. The genes were also screened for their chromosomal information from the Ensembl *Sus scrofa* database. After normalization, removal of values < 0 and obtaining the gene chromosomal location information, 17,726 (breed-specific), and 17,697 (FE-specific) genes were retained in each group. As a quality control for IntLIM, we filtered out genes with the lowest 5% of the variation, which gave 16,839 genes (breed-specific) and 16,812 genes (FE-specific) that were subjected to IntLIM analysis. A schematic representation of the study design and analysis steps are given in Figure 1. We performed the PCA analysis (Figure S1) on the filtered metabolome-transcriptome data, which included 405 metabolites, and 16,839 genes (breed-specific), and 16,812 genes (FE-specific). The results of PCA analysis for the metabolites-genes based on breed and FE groups are shown in Figure S1.

**Figure 1.** Schematic representation of the study design and analysis steps.

#### *2.2. Gene Metabolite Interaction of Breed-Specific and FE-Specific Groups*

From the IntLIM analysis, we identified gene-metabolite pairs that have a strong association with respect to the breed (Duroc and Landrace) and FE (low and high FE) groups. For the breed-specific groups, all possible combinations of gene-metabolite pairs (6,819,795 model runs) were evaluated, using Duroc and Landrace as the breed-group. Based on this approach, we identified 21 gene-metabolite associations (false discovery rate—FDR adjusted *p*-value ≤ 0.1 and correlation difference effect size > 0.1) (Table 1). Clustering these pairs by the direction of association (positive and negative correlations) within each breed group revealed two major clusters (Figure 2a) in each breed. First, the Duroc correlated/Landrace anti-correlated cluster consists of seven gene-metabolite pairs (three unique metabolites and five unique genes) with a high positive correlation in Duroc and low or negative correlation in Landrace (Figure 2a). Second, the Duroc anti-correlated/Landrace correlated cluster consists of 14 gene-metabolite pairs (10 unique metabolites and nine unique genes) with relatively high negative correlations in Duroc and positive correlations in Landrace. The two top-ranked gene-metabolite pairs (ranked in descending order of absolute value of Spearman correlation difference between Duroc and Landrace) in the Duroc correlated and anti-correlated clusters were ENSSSCG00000028124 (*SNRPN*)—Rhodamine B (Figure 2b) and ENSSSCG00000000401 (*GLS2*)—cystathionine ketimine (Figure 2c) respectively. *SNRPN* and Rhodamine B are positively correlated in Duroc (r = 0.7) but negatively correlated in Landrace (r = −0.5) (Figure 2b). *GLS2* and cystathionine ketimine are negatively correlated in Duroc (r = −0.9), but positively correlated in Landrace (r = 0.2) (Figure 2c).

Similarly, we used IntLIM for the FE-specific group and evaluated all possible combinations of gene-metabolite pairs (6,808,860 models), with low and high FE as a binary phenotype. With this approach, we identified 12 FE-specific gene-metabolite correlations (FDR adjusted interaction *p*-value ≤ 0.1, and a Spearman correlation difference > 0.1) involving eight unique gene and metabolites each (Table 2). The heat map with gene-metabolite Spearman correlation for low and high FE group showed a clear separation between the two groups (Figure 3a). The high FE-correlated cluster of 12 gene-metabolite pairs (eight unique genes and metabolites with high correlations in high-FE groups) were negatively correlated with the low-FE group. The two gene-metabolite pairs (ranked in descending order of absolute value Spearman correlation difference between high and low FE group) in high-FE correlated clusters were ENSSSCG00000025106 (*THNSL2*)—pyrocatechol (Figure 3b) and ENSSSCG00000036609 (*TBXT*)—ketoleucine (Figure 3c), respectively. Both pairs showed a positive correlation in high-FE group (r = 0.6, r = 0.5) while showed a negative correlation in the low-FE group (r = −0.7, r = −0.3) (Figure 3b,c).

**Figure 2.** Results of IntLIM applied to breed-specific groups. (**a**) Clustering of 21 identified gene-metabolite pairs (FDR adjusted *p*-value of interaction coefficient < 0.1, Spearman correlation difference > 0.1 in Duroc and Landrace breeds, (**b**) Gene-metabolite difference in ENSSSCG00000028124 (*SNRPN*)—rhodamine B (FDR adjusted *p*-value = 0.1, Duroc Spearman correlation = 0.7, Landrace Spearman correlation = −0.5), (**c**) Gene-metabolite difference in ENSSSCG00000000401 (*GLS2*)—cystathionine ketimine (FDR adjusted *p*-value = 0.01, Duroc Spearman correlation = −0.9, Landrace Spearman correlation = 0.2).

**Figure 3.** Results of IntLIM applied to FE-specific groups. (**a**) Clustering of 12 identified gene-metabolite pairs (FDR adjusted *p*-value of interaction coefficient < 0.1, Spearman correlation difference > 0.1 in high and low FE groups, (**b**) Gene-metabolite difference in ENSSSCG00000025106 (*THNSL2*)—Pyrocatechol (FDR adjusted *p*-value = 0.06, High-FE Spearman correlation = 0.6, Low-FE Spearman correlation = −0.7), (**c**) Gene-metabolite difference in ENSSSCG00000036609 (*TBXT*)—ketoleucine (FDR adjusted *p*-value = 0.08, High-FE Spearman correlation = 0.5, Low-FE Spearman correlation = −0.3).



Duroc\_cor: Spearman correlation in Duroc; Landrace\_cor: Spearman correlation in Landrace; Abs diff.corr: the absolute difference in correlation between Duroc and Landrace; Pval: *p*-value (*p* < 10− was selected as the cut-off value); FDRadjPval: FDR adjusted *p*-value.


**Table 2.** Gene-metabolite interaction pairs from IntLIM for the FE-specific groups.

High\_cor: Spearman correlation in high FE group; Low\_cor: Spearman correlation in low FE group; Abs diff.corr: the absolute difference in correlation between the high and low FE group; Pval: *p*-value (*p* < 10− was selected as the cut-off value); FDRadjPval: FDR adjusted *p*-value.

#### *2.3. Pathway and Gene Ontology Over-Representation Analysis*

We identified the pathways associated with the unique metabolites in each cluster identified from breed-specific (21) and FE-specific (12) interactions. The three unique metabolites from Duroc correlated/Landrace anti-correlated clusters were associated with arginine and proline metabolism (*p*-value = 0.02). Furthermore, the ten unique metabolites from Duroc anti-correlated/Landrace correlated cluster were associated with valine-leucine-isoleucine biosynthesis (*p*-value = 0.01) and valine-leucine-isoleucine degradation (*p*-value = 0.07) along with arginine and proline metabolism (*p*-value = 0.07). The eight unique metabolites from high-FE correlated/low-FE anti-correlated cluster were associated with valine-leucine-isoleucine biosynthesis (*p*-value = 0.01) and degradation (*p*-value = 0.07). The pathways associated with the metabolites in breed-specific and FE-specific clusters for unique metabolites are given in Supplementary Table S1a.

Unique and mappable genes from each group (breed-specific—each cluster, and FE-specific) were screened by using GeneMania to generate a composite functional association network that includes all the evidence of co-functionality. From the breed-specific group, unique genes (4 genes) from Duroc correlated/Landrace anti-correlated cluster (Table 1) mapped to 20 genes based on co-functionality from GeneMania (Table S1b). The gene-ontology enrichment analysis of the identified 24 genes (unique genes from Table 1 and co-functional genes from GeneMania) revealed enrichment of the regulation of hemopoiesis, response to thyroid hormone, and the sphingolipid catabolic process (Table S1c). These genes were enriched for the sphingolipid metabolism KEGG pathway (adjusted *p*-value corrected with Bonferroni step down ≤ 0.05) (Supplementary Table S1d). Unique and mappable genes (6 genes) identified from Duroc anti-correlated/Landrace correlated cluster (Table 1) were co-functional with 20 genes based on GeneMania (Table S1b). The gene ontology enrichment analysis of these 26 genes (unique genes from Table 1 and co-functional genes from GeneMania) revealed the ER to Golgi vesicle-mediated transport and membrane fusion (Table S1c) as an enriched biological process. Butanoate metabolism, alanine-aspartate-glutamate metabolism, and valine-leucine-isoleucine degradation were significantly enriched KEGG pathways from the Duroc anti-correlated/Landrace correlated cluster (Supplementary Table S1d). Similarly, from the FE-specific group, unique mapped genes (7 genes) from high-FE correlated/low-FE anti-correlated clusters (Table 2) were co-functional with 20 genes identified from GeneMania (Table S1b). These genes were involved with the cGMP metabolic process, purine nucleotide biosynthesis, and phosphorus-oxygen lyase activity pathways (Table S1c). The top significant KEGG pathway enriched was the cGMP-PKG signaling pathway (Supplementary Table S1d).
