**2. Results**

#### *2.1. First Component Score and Significant Metabolic Pathways of 45 Metabolites*

The partial least squares-discriminant analysis (PLS-DA) results revealed that the first component score (component 1) explained more than 75% variation of all 45 metabolites (Figure 1A). It showed that metabolite values of Duroc were higher than those of Landrace, the same as metabolites from second sampling time higher than those from first sampling time. In addition, the Duroc and Landrace pigs were clearly stratified, especially using the metabolite values between Duroc from first sampling time and Landrace from second sampling time (Figure 1A). The most significant metabolic pathways were the aminoacyl-tRNA biosynthesis; following by the arginine biosynthesis; the arginine and proline metabolism; and the alanine, aspartate, and glutamate metabolism (Figure 1B). As the pathway impact of the aminoacyl-tRNA biosynthesis was zero, we discarded this significant pathway and only used the metabolites enriched in the other three significant pathways for GWAS (Table 1). Thus, the metabolite means for 5 compounds in the arginine biosynthesis (arginine, aspartic acid, citrulline, glutamic acid, and ornithine); 5 compounds in the arginine and proline metabolism (arginine, glutamic acid, ornithine, proline, and pyruvic acid); and 4 compounds in the alanine, aspartate, and glutamate metabolism

(alanine, aspartic acid, glutamic acid, and pyruvic acid) metabolites were calculated and shown in Table 1.

**Figure 1.** (**A**) Partial least squares-discriminant analysis (PLS-DA) of 45 metabolites. Note: D/L with first/second indicates the sampling time of Duroc/Landrace pigs. (**B**) Metabolic pathways for 45 metabolites using *Homo sapiens* as the library. Note: The size and color of the circles for each pathway indicate the matched metabolite ratio and the −log (*p*-value), respectively.

#### *2.2. Genome-Wide Significant SNPs and Gene Annotation*

Metabolite based GWAS for first, second, and combined two sampling times revealed 152 genome-wide significant SNPs (Supplementary Table S1) in association with 17 metabolites (Supplementary Table S2). Unfortunately, no significant SNP was detected in association with first component scores (*p*-values ≥ 2.78 × <sup>10</sup>−6) and metabolites enriched in the significant metabolic pathways (*p*-values ≥ 1.74 × <sup>10</sup>−4); thus, GWAS results of these two scenarios were not listed. Manhattan plots of genome-wide association for isovalerylcarnitine and propionylcarnitine are shown in Figure 2, and Manhattan plots for the other 43 metabolites are shown in the Supplementary Figure S1. Along the whole genome, significant SNPs associated with isovalerylcarnitine and propionylcarnitine from the second sampling time were mainly located on the chromosome one (Figure 2). The overlapped significant SNPs associated with more than two different metabolites were shown in Table 2, where 57 significant SNPs on genome level were associated with isovalerylcarnitine and propionylcarnitine from the second sampling time. In addition, another 3 metabolites (1-hexadecyl-sn-glycero-3-phosphocholine, 1-myristoyl-sn-glycero-3-phosphocholine, lysoPC (16:0)) were also significantly associated with 10 SNPs (Table 2).


**Table 1.** Significant metabolic pathways (False Discovery Rate (FDR) < 0.1) using *Homo sapiens* as the library.


**Table 2.** Common significant single-nucleotide polymorphisms (SNPs) of genome-wide association for more than two different metabolites from first, second,

**Figure 2.** Manhattan plots of genome-wide association for (**A**) isovalerylcarnitine and (**B**) propionylcarnitine. Note: *Y*-axis indicates the log10 (*p*-value). Blue dotted and red solid lines indicate the genome-wide threshold of 0.05 and 0.01 after Bonferroni multiple testing, respectively. The three tracks indicate the metabolites from first sampling time, second sampling time, and combined two sampling times from outside to inside.

After annotation of significant SNPs to the neighboring genes and gene components (Sscrofa10.2/susScr3), we found that 90 significant SNPs were within a 500-kb window of 52 neighboring genes (Supplementary Table S1) and that 6 significant SNPs were directly located in the gene components of 5 genes (Table 3). For example, if we only consider the SNPs on chromosome one, we found 29 significant SNPs were near 9 genes (Supplementary Table S1), whereas ALGA0004000, ALGA0004041, and ALGA0004042 were located in the introns of *FBXL4* and *CCNC* (Table 3). These results show that these genes may be involved in regulating abundance of the metabolites that are significantly different between high and low RFI pigs. Between using porcine RefSeq database of Sscrofa10.2/susScr3 and Sscrofa11.1/susScr11, the results of significant SNPs annotated to the genes overlapped greatly, but SNPs had different distances to the annotated genes between two versions (Supplementary Table S1). In Table 3, we found that the annotations of ALGA0004042 and ALGA0061605 to *CCNC* and *MTRF1* were changed from 9th intron and 5th intron to 8th intron and 9th intron, respectively, when we used the Sscrofa11.1/susScr11 database.



#### *Metabolites* **2020**, *10*, 201

*Metabolites* **2020**, *10*, 201

The linkage disequilibrium (LD) pattern for all significant SNPs is shown in the Supplementary Figure S2. From the LD results for 58 significant SNPs on chromosome one, we found that 51 significant SNPs associated with isovalerylcarnitine (*p*-value = 2.79 × <sup>10</sup>−8) and propionylcarnitine (*p*-value = 8.32 × <sup>10</sup>−10) from second sampling time were in strong LD (Figure 3). Among the 58 significant SNPs, five of them were not in LD with the other 53 significant SNPs (Supplementary Figure S2), so they were excluded in the haplotype visualization in the Figure 3. In detail, SNPs annotated to *LOC780435* (NM\_001078684), *FHL5* (NM\_001243314), *FBXL4* (NM\_001171752), *CCNC* (NM\_001190160)/*MCHR2* (NM\_001044609), and *SIM1* (NM\_001172585) were in block 2, block 4, block 6, block 8, and block 9/10, respectively. Furthermore, ALGA0004000 in the 6th intron of *FBXL4* was in LD of block 6, together with another five SNPs (INRA0002726, MARC0075306, ALGA0003995, ALGA0004002, and ALGA0004005) that were located in the intergenic regions of *FBXL4*. Especially, three SNPs in strong LD consisted of block 8 with two SNPs (ALGA0004041 and ALGA0004042) located in the second and ninth intron of *CCNC* (Figure 3, Table 3, and Supplementary Table S1). The number of significant SNPs in strong LDs of the other chromosomes was less than the significant SNP number on chromosome one (Supplementary Figure S2). Notably, MARC0110390 in the 7th intron region of *SFXN1* (NM\_001098602) on chromosome two and ALGA0061605 in the 5th intron region of *MTRF1* (NM\_001243580) on chromosome eleven were not in the LD with other SNPs. However, ASGA0093565 in the 8th intron region of *DNAJC6* (NM\_001145378) was in strong LD with WU\_10.2\_6\_135312468 that was annotated to *LEPROT* (NM\_001145388) (Supplementary Table S1 and Supplementary Figure S2).

**Figure 3.** Linkage disequilibrium (LD) pattern for 53 significant SNPs on chromosome one. Note: the solid line triangle indicates LD. One square indicates LD level (r2) between two SNPs, and the squares are colored by the D' & LLOR standard scheme. D' & LLOR standard scheme is that red indicates LLOR > 2, D' = 1; pink indicates LLOR > 2, D < 1; blue indicates LLOR < 2, D' = 1; and white indicates LLOR < 2, D' < 1. LLOR is the logarithm of likelihood odds ratio and the reliable index to measure D'.

#### *2.3. Gene and Metabolite Interaction Network*

The most significantly enriched gene-based pathways were the human papillomavirus infection (ssc05165) with five genes (i.e., *CCND2*, *CTNNB1*, *JAK1*, *LAMC1*, and *NFKB1*), followed by the PI3K-Akt signaling pathway (ssc04151) with five genes (i.e., *CCND2*, *F2R*, *JAK1*, *LAMC1*, and *NFKB1*) and the hepatitis C (ssc05160) with four genes (i.e., *CLDN8*, *CTNNB1*, *JAK1*, and *NFKB1*) (Figure 4A). Based on the gene–gene interaction network analysis, *CCNC* was in good connection with *CDK8*, *CDK3*, and *N6AMT1* whereas *N6AMT1* was linked to *FBXL4* (Figure 4B). Unfortunately, no gene–metabolite interaction network was found in this study. After the clustering of the SNP-related gene component-associated metabolites (Table 3), we found that aspartic acid, 1-hexadecyl-sn-glycero-3-phosphocholine, 1-myristoyl-sn-glycero-3-phosphocholine, and lysoPC(16:0) were clustered in the lower cluster while the upper cluster included the metabolites of isovalerylcarnitine, propionylcarnitine, and pyruvic acid (Figure 4C). Results show that metabolites from Duroc pigs have higher values in the upper cluster than those from lower cluster, but the metabolite values of Landrace pigs are higher in the lower cluster (Figure 4C). Afterwards, we investigated the metabolite values of aspartic acid, isovalerylcarnitine, propionylcarnitine, and pyruvic acid for which the associated significant SNPs were in the introns of *MTRF1*, *FBXL4*/*CCNC*, *SFXN1* (Table 3). Generally, propionylcarnitine from the low RFI group had higher values while other three metabolite values in the high RFI group seemed higher, but they are not significantly different between low and high RFI groups (*p*-value > 0.05) (Figure 4D).

#### *Metabolites* **2020**, *10*, 201

**Figure 4.** Gene pathway, metabolite cluster, and the interaction network: (**A**) Pathway for significant SNP-related genes. (**B**) Network for the genes in which significant SNPs were annotated to gene components. (**C**) Heatmap cluster for the metabolites that were associated with significant SNPs annotated to gene components. (**D**) Metabolites (i.e., aspartic acid, isovalerylcarnitine, propionylcarnitine, and pyruvic acid) values in high and low residual feed intake (RFI) pigs associated with the genes in which significant SNPs annotate to gene components. Note: The high RFI pigs and low RFI pigs were from left and right parts of all the pigs (n = 108) with one SD of actual RFI values.
