Integrative Analysis of Transcriptomic and Lipidomic Profiles Reveal Differential Mechanism of Subcutaneous Adipose Tissue among Ningxiang, Berkshire and Their Offspring Pigs

1. Introduction

Pork is a major source of amino acids and fatty acids, especially in China. Healthy, safe and nutritious high-quality and special flavor pork is more and more popular for consumer, and has huge demand. Adipose tissue is a dynamic metabolic and endocrine organ, which also stores fat and secrete biological activity substances, such as leptin and lipid mediators [1]. It was considered that body fat accumulation may be the result of the balance among the absorption of fat, de novo fat synthesis, fat metabolism and catabolism [2]. Fatty acids composition of adipose tissue contributed greatly to the flavor, quality and nutritional value of meat [3]. The lipid content and the fatty acids composition in adipose tissue were strongly affected by genotype and nutrient supply [4].

Indigenous pigs, with excellent meat quality and special flavor, is more and more popular with consumers and researchers compared to modern pig breeds [5,6]. While low growth rate and lean meat rate limit the large-scale breeding of indigenous pigs. Fat deposition exhibited huge differences between indigenous pig breeds and modern pig breeds [7]. The selective breeding for high lean meat rate and low backfat thickness resulted in a strongly reduced lipogenic potential in the modern pig breeds, while indigenous pig breeds preserved high fat deposition capacity [4]. New varieties and breeds with indigenous pig lineage, which had high growth rate, high lean meat rate, low backfat thickness, excellent flavor, were strongly requirement for market. These new varieties and breeds were cultivated by hybrid breeding between indigenous pig breeds and modern pig breeds, such as Xiangcun black pigs, Sutai pigs. Nowadays, how to apply multi-omics data to accelerate pig breeding may be a key problem. Recently, multi-omics association analysis had brought opportunities for exploring the molecular mechanisms of complex traits and for more efficient mining of candidate genes for complex traits [8,9]. The clearly molecular mechanisms of given trait and candidate gene resources may contribute to pig breeding. The molecular mechanisms of differences in adipose tissue between indigenous pig breeds and modern pig breeds was complex. Several studies compared the lipogenic enzyme activities of subcutaneous adipose tissue in multiple groups (Basque pigs vs. Large White pigs, Alentejano pigs vs. Large White pigs, Meishan pigs vs. Large White pigs), and the result showed that Basque pigs and Alentejano pigs exhibit higher activities of ACACA, G6PDH and Meishan pigs exhibit lower SCD activities compared to Large White pigs [10,11,12]. Thus, the expression and activities difference of these key genes among different varieties may contribute to the phenotypic difference.

The Ningxiang pigs, one of four well-known indigenous breeds in China, has abundant polyunsaturated fatty acid content and superior meat quality [13]. The Ningxiang pigs is a typical fatty pig breed, with a slow growth rate and strong fat deposition ability. Berkshire pigs is a typical lean pig breed, famous for its high lean meat rate and fast growth rate. Ningxiang pigs, acted as male parent, crossed to Berkshire pigs, the meat quality of F₁ pigs was superior. Compared with Ningxiang pigs, F₁ pigs had narrower backfat thickness. The obvious phenotype differences existed in the subcutaneous adipose tissue among Ningxiang, Berkshire and F₁ pigs. multi-omics association analysis had proved to be an effective method for investigating molecular mechanisms fat deposition [14]. However, there is little multi-omics data on the difference among Ningxiang, Berkshire and their offspring pigs, and molecular mechanisms of fat deposition difference among Ningxiang, Berkshire and their offspring pigs still unknown. Therefore, this study aimed to research the molecular mechanisms of fat deposition difference among Ningxiang, Berkshire and F₁ pigs and identify the key genes and lipids. This study also provided insight into selection of backfat thickness and fat composition of adipose tissue for pig breeding.

2. Materials and Methods

2.1. Animals

A total forty-two pigs (fourteen purebred Ningxiang pigs, fourteen Berkshire pigs and fourteen hybrid progenies of Ningxiang pig × Berkshire pig) from the same farm were used in this study. All pigs were raised on the same farm and housed under standard management conditions at the Dalong livestock Co. Ltd., in Hunan province, China, with the same diet and free access to water. All pigs were healthy. The subcutaneous adipose tissue sample of pigs were collected from Chu Weixiang Slaughtering and Cutting Plant in Ningxiang City, Hunan Province. We randomly selected six Ningxiang pigs with a left carcass weight 32–38 kg, six Berkshire pigs with a left carcass weight 41–47 kg, and six hybrid progenies of Ningxiang pig X Berkshire pig (F₁ pigs) with a left carcass weight 37–41 kg, and collected the subcutaneous adipose between the 6th and 11th ribs samples from each pig. The samples were crushed and packed, and then stored in a refrigerator at -80 ℃ with constant room temperature and humidity. The experiment was approved by the Animal Care Committee of the Institute of Subtropical Agriculture.

2.2. RNA Extraction

The total RNA was extracted using the TRIzol reagent (Life technologies, California, USA). RNA concentration and purity was measured using NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, DE). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). A total amount of 1 μg RNA per sample was used as input material for the RNA sample preparations.

2.3. mRNA Sequencing and Transcriptome Data Analysis

RNA-seq was performed by Biomarker Technology Co., Ltd. (Beijing, China). Sequencing libraries were generated using NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, USA), and then library quality was assessed on the Agilent Bioanalyzer 2100 system. After cluster generation, the library preparations were sequenced on an Illumina platform and paired-end reads were generated. The raw data were filtered using the FastQC program (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), clean reads were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated to ensure high quality of clean reads. These clean reads were then mapped to the pig reference genome Sus scrofa11.1 (http://www.ensembl.org/info/data/ftp/index.html) by Hisat2 software (http://www.ccb.jhu.edu/software/hisat). Only reads with a perfect match or one mismatch were further analyzed and annotated based on the pig reference genome. Differential expression analysis in the subcutaneous adipose tissue group (four subcutaneous adipose tissue of Ningxiang samples, four subcutaneous adipose tissue of Berkshire samples, and four subcutaneous adipose tissue of hybrid progenies) was performed using the DESeq2 [40]. Differential gene screening conditions were: |log2Foldchange| > 1, and false discovery rate (FDR) < 0.05. Gene Ontology (GO) enrichment analysis of DEGs was implemented by the GOseq R package [41]. KOBAS [42] software was used to test the statistical enrichment of differential expression genes in KEGG pathways. P-values were adjusted with FDR method and adjusted P-values < 0.05 were considered significant.

2.4. Lipids Extration and LC-MS/MS Analysis

The 100 μL sample was mixed with 1mL extraction solvent (MTBE: MeOH =3:1, v/v) containing internal standard mixture, and then vortexed for 15 min and sonicated in an ice bath for 10 min, 200 μL of ultrapure water was added. 200 μL of the extration solution was collected and evaporated using a vacuum concentrator.The dry extract was reconstituted using 400 μL mobile phase B prior to LC-MS/MS analysis. The sample extracts were analyzed using an LC-ESI-MS/MS system (UPLC, ExionLC AD，https://sciex.com.cn/; MS, QTRAP^® System, https://sciex.com/). The analytical conditions were as follows, UPLC: column, Thermo Accucore™ C30 (2.6 μm, 2.1 mm*100 mm i.d.); solvent system, A: acetonitrile/water (60/40,V/V, 0.1% formic acid, 10 mmol/L ammonium formate), B: acetonitrile/isopropanol (10/90 V/V, 0.1% formic acid, 10 mmol/L ammonium formate); gradient program, A/B（80:20, V/V）at 0 min, 70:30 v/v at 2.0 min, 40:60 v/v at 4 min, 15:85 v/v at 9 min, 10:90 v/v at 14 min, 5:95 v/v at 15.5 min, 5:95 v/v at 17.3 min, 80:20 v/v at 17.3 min, 80:20 v/v at 20 min; flow rate, 0.35 ml/min; temperature, 45 °C; Injection volume: 2μl. The effluent was alternatively connected to an ESI-triple quadrupole-linear ion trap (QTRAP)-MS. LIT and triple quadrupole (QQQ) scans were acquired on a triple quadrupole-linear ion trap mass spectrometer (QTRAP). The ESI source operation parameters were as follows: ion source, turbo spray; source temperature 500 °C; ion spray voltage (IS) 5500 V（Positive）,-4500 V (Neagtive); ion source gas 1 (GS1), gas 2 (GS2), curtain gas (CUR) were set at 45, 55, and 35 psi, respectively; the collision gas (CAD) was medium. Instrument tuning and mass calibration were performed with 10 and 100 μmol/L polypropylene glycol solutions in QQQ and LIT modes, respectively. QQQ scans were acquired as MRM experiments with collision gas (nitrogen) set to 5 psi. DP and CE for individual MRM transitions was done with further DP and CE optimization. A specific set of MRM transitions were monitored for each period according to the lipids eluted within this period.

2.5. Statistical Analysis

Unsupervised PCA (principal component analysis) was performed by statistics function prcomp within R (Version 3.5.1, www.r-project.org). The data was unit variance scaled before unsupervised PCA. Variable important in the projection (VIP) were calculated in the OPLS-DA model, and generated using R package MetaboAnalystR (Version 1.01). Differential lipids were determined by VIP (VIP ≥ 1) and absolute Log2FC (|Log2FC| ≥ 1.0). In order to avoid overfitting, a permutation test (200 permutations) was performed. Then the differential lipids were annotated using KEGG Compound database (http://www.kegg.jp/kegg/compound/), and mapped into biochemical pathways through KEGG enrichment analysis. Pathways with significantly regulated metabolites mapped to were then fed into MSEA (metabolite sets enrichment analysis), their significance was determined by hypergeometric test’s p-values. Heat maps of DEGs was performed using the Metware Cloud, a free online platform for data analysis (https://cloud.metware.cn).

3. Results

3.1. The Backfat Thickness of Ningxiang pigs, Berkshire pigs and F₁ Pigs

We measured the backfat thickness of pigs of the same age and found that the backfat thickness of Ningxiang pigs was significantly higher than that of and F₁ pigs, the backfat thickness of F₁ pigs was significantly higher than that of Berkshire pigs (Figure 1).

3.2. Lipidomic Data Analysis

To investigate the difference in lipid metabolite compositions of the subcutaneous adipose tissue among Ningxiang pigs, Berkshire pigs and F₁ pigs. LC-MS/MS analysis was performed. A total of 1293 lipids were detected, including 506 glycerphospholipids (GP), 421 glycerolipids (GL), 226 sphingolipids (SP), 88 triglycerides (TG), 40 fatty acyls (FA), 10 sterol lipids (ST), 2 prenol lipid (PR). In the detected lipid subclasses, triglycerides (TG), diglycerides (DG), Phosphatidylethanolamine (PE), Phosphatidylcholine (PC), Phosphatidylserine (PS), ceramide (Cer-AS) were abundant in the subcutaneous adipose tissue (Figure S1). Principle component analysis (PCA) indicated that three different groups were partly separated in PC1 × PC2 score plots (Figure 2A). Further orthogonal partial least squares discriminant analysis (OPLS-DA) showed that there were apparent distinctions among three groups, and small distinctions within group (Figure 2B-D). Besides, the overall sample cluster heat map, which was drawed by the Z-sore for the expression level of all samples, and K-Means analysis suggested that the relative lipid content of the subcutaneous adipose tissue of Berkshire pigs was higher than that of Ningxiang pigs and F₁ pigs for the detected lipids in this study (Figure 2E, Figure S2). Then, SCLs were filtered and determined by Variable Importance in Projection (VIP) value ≥ 1.0 and |Log2(fold change)| ≥ 1.0. 423 SCLs were screened in the subcutaneous adipose tissue between Berkshire pigs and Ningxiang pigs. Compared with Berkshire pigs, Ningxiang pigs had 412 downregulated, primarily GPs, GLs, SPs, DGs, FAs, STs, PR (238, 93, 44, 24, 10, 2 and 1, respectively), only 11 upregulated including 8 TGs, 2 SPs and 1 GP, in the subcutaneous adipose tissue (Figure 3A); 252 SCLs were markedly changed in the subcutaneous adipose tissue of F₁ pigs compared to that of Berkshire pigs. F₁ pigs had 251 downregulated, containing 3 FAs, 13 GLs, 219 GPs, 12 SPs and 4 TGs, only 1 SP upregulated in the subcutaneous adipose tissue (Figure 3B); For the groups of Ningxiang pigs vs. F₁ pigs, 50 SCLs were screened. Compared with Ningxiang pigs, F₁ pigs had just 7 downregulated including 4 GPs, 1 SP, 1 TG, 1 GL, 43 upregulated containing 12 GLs, 11 GPs, 11SPs, 4STs, 3FAs, and 2TGs in the subcutaneous adipose tissue (Figure 3C).

Subsequently, KEGG enrichment analysis was performed to uncover biological mechanisms of the subcutaneous adipose tissue among Ningxiang pigs, Berkshire pigs and F₁ pigs. With respect to the groups of Berkshire pigs vs. Ningxiang pigs, 423 SCLs were significantly annotated into 108 pathways, and the majority of SCLs were involved in fatty acid and amino acid metabolism related pathway such as glycerophospholipid metabolism, linoleic acid metabolism, arachidonic acid metabolism, alpha-linolenic acid metabolism, inositol phosphate metabolism, and glycine, serine and threonine metabolism (Figure 3D). For the group of Berkshire pigs vs. F₁ pigs, 252 SCLs were significantly annotated into 38 pathways, including glycerophospholipid metabolism, linoleic acid metabolism, arachidonic acid metabolism, alpha-linolenic acid metabolism, glycine, serine and threonine metabolism, and cAMP signaling pathway (Figure 3E). For the groups of Ningxiang pigs vs. F₁ pigs, 50 SCLs were significantly annotated into 31 pathways, including teichoic acid biosynthesis, inositol phosphate metabolism, steroid biosynthesis, ovarian steroidogenesis, glycerolipid metabolism, and phosphatidylinositol signaling system (3F).Thus, there were 423, 252, and 50 SCLs up- or downregulated in the groups of Berkshire pigs vs. Ningxiang pigs, Berkshire pigs vs. F₁ pigs, and Ningxiang pigs vs. F₁ pigs, respectively (Figure 4A). A total of 13 SCLs were commonly identified in three pairwise comparisons (Figure 4A). Additionally, there were 19, 16, and 7 KEGG pathways (P < 0.05) enriched in the groups of Berkshire pigs vs. Ningxiang pigs, Berkshire pigs vs. F₁ pigs, and Ningxiang pigs vs. F₁ pigs, respectively (Figure 4B). Interestingly, fourteen pathways were commonly enriched in Berkshire pigs vs. Ningxiang pigs and Berkshire pigs vs. F₁ pigs, including seven (kO05231, kO00563, kO00591, kO00592, kO00590, and kO00564) involving in lipid metabolism, one (kO00260) involving glycine, serine and threonine metabolism, one (kO04723) involving retrograde endocannabinoid signaling (Figure 4B). Four pathways were commonly enriched in Berkshire pigs vs. Ningxiang pigs and Ningxiang pigs vs. F₁ pigs, including inositol phosphate metabolism (kO00562), teichoic acid biosynthesis (kO00552), and phosphatidylinositol signaling system (kO04070) (Figure 4B). Besides, three pathways were uniquely enriched in Ningxiang pigs vs. F₁ pigs, including steroid biosynthesis (kO00100), ovarian steroidogenesis (kO04913) and bile secretion (kO04976).

3.3. Transcriptome Data Analysis

Transcriptome sequencing of the subcutaneous adipose tissue of Ningxiang pigs, Berkshire pigs and F₁ pigs were performed, in order to verified the reliability of the results, their correlations were analyzed. The correlation analysis indicated that strong correlation within groups, and that the correlation between the groups was slightly lower (Figure S3A). As the heatmap was showed in Figure S3B, C and D, the results showed that the expression patterns of DEGs consistent with the sample groups. We counted the DEGs among three difference pigs. We found 395 up-regulated genes and 285 down-regulated genes, and a total of 680 DEGs in the subcutaneous adipose tissue of the group of Berkshire pigs vs. Ningxiang pigs (Figure 5A). There were 317 up-regulated genes, 275 down-regulated genes and a total of 592 DEGs in the subcutaneous adipose tissue of F₁ pigs compared to that of Ningxiang pigs (Figure 5B). There were 230 up-regulated genes, 154 down-regulated genes, and a total of 384 DEGs in the subcutaneous adipose tissue of F₁ pigs compared to that of Berkshire pigs (Figure 5C). Interestingly, the number of DEGs form the group of Berkshire pigs vs. Ningxiang pigs were largest among three difference comparisons.

Based on the identified DEGs, functional annotation and classification were performed using GO analysis. Genes and DEGs were classified into three primary categories, such as biological process (BP), cellular component (CC), and molecular function (MF). For the groups of Berkshire pigs vs. Ningxiang pigs, after GO enrichment of DEGs, the biological processes were mainly enriched in the process related to xenobiotic metabolic, organic acid metabolic and L-serine biosynthetic (Figure S4A), the cellular component were mainly enriched in collagen, lipid droplet and cytoplasm (Figure S4B), the molecular function mainly enriched in lipid binding, transmembrane receptor protein serine/threonine kinase activity (Figure S4C). For the groups of Ningxiang pigs vs. F₁ pigs, after GO enrichment of DEGs, the biological processes were mainly enriched in the process related to translation, brown fat cell differentiation and ATP synthesis coupled proton transport (Figure S4D), the cellular component were mainly enriched in mitochondrial respiratory chain complex I, collagen trimer, extracellular matrix (Figure S4E), the molecular function mainly enriched in cytochrome-c oxidase activity, structural constituent of ribosome and NADH dehydrogenase activity (Figure S4F). For the groups of Berkshire pigs vs. F₁ pigs, after GO enrichment of DEGs, the biological processes were mainly enriched in the process related to response to xenobiotic stimulus, immune response and stem cell fate specification (Figure S4G), the cellular component were mainly enriched in extracellular region, cytoplasm and collagen trimer (Figure S4H), the molecular function mainly enriched in chemokine activity, phosphoglycerate dehydrogenase activity and sterol esterase activity (Figure S4I). Then, the DEGs were annotated with the KEGG pathway. There were 52, 30, and 35 KEGG pathways (P < 0.05) enriched in the groups of Berkshire pigs vs. Ningxiang pigs, Berkshire pigs vs. F₁ pigs, and Ningxiang pigs vs. F₁ pigs, respectively (Figure 6). Interestingly, twelve pathways were commonly enriched in the groups of Berkshire pigs vs. Ningxiang pigs and Berkshire pigs vs. F₁ pigs, seven pathways were commonly enriched in the groups of F₁ pigs vs. Ningxiang pigs and Berkshire pigs vs. F₁ pigs, twenty pathways were commonly enriched in the groups of Berkshire pigs vs. Ningxiang pigs and Ningxiang pigs vs. F₁ pigs (Figure 6). For the groups of Berkshire pigs vs. Ningxiang pigs, KEGG analysis of the subcutaneous adipose tissue yielded significant pathways, such as carbon metabolism, fatty acid metabolism, fatty acid biosynthesis, and PI3k-Akt signaling pathway (Figure 7A, Table S1). For the groups of Ningxiang pigs vs. F₁ pigs, KEGG analysis revealed that significant pathways were closely related to oxidative phosphorylation, thermogenesis and protein digestion and absorption (Figure 7B). For the groups of Berkshire pigs vs. F₁ pigs, KEGG analysis indicated that significant pathways were closely related to TNF signaling pathway, cysteline and methionine metabolismand and Glycine, serine and threonine metabolism (Figure 7C).

3.4. Joint Analysis of the Transcriptome and Lipidome Data

A joint analysis of the transcriptome and lipidome were carried out to identify the origin of the transcriptomic and lipidomic differences among Ningxiang pigs, Berkshire pigs and F₁ pigs. A total of 680 DEGs (395 upregulated, 285 downregulated) and 423 SCLs (11 upregulated, 412 downregulated) were identified in the group of Berkshire pigs vs. Ningxiang pigs, 384 DEGs (230 upregulated, 154 downregulated) and 252 SCLs (1 upregulated, 251 downregulated) were identified in the group of Berkshire pigs vs. F₁ pigs, 592 DEGs (317 upregulated, 275 downregulated) and 50 SCLs (43 upregulated, 7 downregulated) were identified in the group of Ningxiang pigs vs. F₁ pigs. Furthermore, pathway annotation was performed, the significant pathways were annotated with the DEGs and SCLs in the groups of Berkshire pigs vs. Ningxiang pigs (Figure 8A), Berkshire pigs vs. F₁ pigs (Figure 8B), and Ningxiang pigs vs. F₁ pigs (Figure 8C), respectively. Interestingly, the most significant pathways enriched by DEGs and SCLs in both the groups of Berkshire pigs vs. Ningxiang pigs and the groups of Berkshire pigs vs. F₁ pigs are the lycine, serine and threonine metabolism (Figure 8A,B).

3.5. Candidate genes and Lipids Screening with the Combined Transcriptome and Lipidome Data

KO00260 ‘the lycine, serine and threonine metabolism’ pathways containing both DEGs and SCLs in the group of Berkshire pigs vs. Ningxiang pigs (Figure S5A) and Berkshire pigs vs. F₁ pigs (Figure S5B), respectively. For the group of Berkshire pigs vs. Ningxiang pigs, KO00260 ‘the lycine, serine and threonine metabolism’ pathways related DEGs phosphoglycerate dehydrogenase (PHGDH), phosphoglycerate dehydrogenase like protein (LOC100156167), and glycine hydroxymethyltransferase (SHMT1) were significantly downregulated, phosphoserine phosphatase (PSPH) and membrane primary amine oxidase-like (LOC110256000) were significantly upregulated (Figure 9A). In addition, the content of the KO00260 ‘the lycine, serine and threonine metabolism’ pathways related SCLs (32 Phosphatidylserines) of the subcutaneous adipose tissue were reduced in Berkshire pigs vs. Ningxiang pigs (Figure S5C). For the group of Berkshire pigs vs. F₁ pigs, PHGDH was significantly downregulated, LOC110256000 and serine dehydratase (SDS) were significantly upregulated (Figure 9B). Moreover, the content of the KO00260 ‘the lycine, serine and threonine metabolism’ pathways related SCLs (31 phosphatidylserines) of the subcutaneous adipose tissue were reduced in the group of Berkshire pigs vs. F₁ pigs (Figure S5D). PHGDH is a key enzyme of the serine synthesis pathway [15]. In this study, PHGDH was significantly downregulated in the subcutaneous adipose tissue of Ningxiang pigs and F₁ pigs compared to that of Berkshire pigs, and the content of phosphatidylserines in the subcutaneous adipose tissue of Ningxiang pigs and F₁ pigs was reduced compared to that of Berkshire pigs. The SCLs and DEGs (pearson correlation coefficient > 0.8 and pvalue < 0.05) in the KEGG pathway were selected to study the correlation between SCLs and DEGs. The result indicated that the expression of PHGDH and LOC100156167 were positive correlated with the content of phosphatidylserines, while the expression of LOC110256000 was negative correlated with the content of phosphatidylserines in the group of Berkshire pigs vs. Ningxiang pigs (Figure 9C) and Berkshire pigs vs. F₁ pigs (Figure 9D).

For the group of Berkshire pigs vs. Ningxiang pigs, the KO00061 ‘fatty acid biosynthesis’ pathways related DEGs, such as, fatty acid synthase (FASN), acetyl-CoA carboxylase (ACACA), carboxyl reductase 4 (CBR4), and medium chain acyl hydrolase (MCH) were upregulated (Figure 10A, C). Besides, the content of the KO00061 ‘fatty acid biosynthesis’ pathways related SCLs (palmitoleic acid) of the subcutaneous adipose tissue was reduced in the group of Berkshire pigs vs. Ningxiang pigs (Figure 10B). While, KO00061 was not significantly enriched in the groups of Ningxiang pigs vs. F₁ pigs and Berkshire pigs vs. F₁ pigs, indicating that after crossing Berkshire pigs with Ningxiang pigs, KO00061 ‘fatty acid biosynthesis’ pathways related genes (FASN, ACACA, CBR4, MCH) of F₁ pigs were upregulated and the relative content of palmitoleic acid was reduced in the subcutaneous adipose tissue of F₁ pigs compared to that of Berkshire pigs.

The kO01040 ‘biosynthesis of unsaturated fatty acids’ pathway was enrichment in the group of Berkshire pigs vs. Ningxiang pigs. The kO01040 ‘biosynthesis of unsaturated fatty acids’ related DEGs, such as elongation of very long chain fatty acids protein 6 (ELOV6) and stearoyl-CoA desaturase (SCD) were upregulated and very long-chain (3R)-3-hydroxyacyl-CoA dehydratase (HACD) was downregulated in the subcutaneous adipose tissue of Ningxiang pigs compared to that of Berkshire pigs (Figure 11A). In addition, the content of the kO01040 ‘biosynthesis of unsaturated fatty acids’ related SCLs, such as palmitoleic acid, linoleic acid, icosadienoic acid and arachidonic acid were reduced in the subcutaneous adipose tissue of the the group of Berkshire pigs vs. Ningxiang pigs (Figure 11A). The KO00591 ‘linoleic acid metabolism’ pathway was significantly enrichment in the group of Berkshire pigs vs. Ningxiang pigs. There were one DEGs, cytochrome P450 family 3 subfamily A46 (CYP3A46) and three SCLs, linoleic acid, arachidonic acid and 9,10-DHOME were involved in linoleic acid metabolism pathway, CYP3A46 was upregulated in the subcutaneous adipose tissue of Ningxiang pigs compared to that of Berkshire pigs, the content of linoleic acid, arachidonic acid and 9,10-DHOME were reduced in the subcutaneous adipose tissue of the the group of Berkshire pigs vs. Ningxiang pigs (Figure 11B).

In kO00590 ‘arachidonic acid metabolism’, the DEGs, gamma-glutamyltranspeptidase (GGT5), cytochrome P450 4F6-like (LOC110255237), microsomal prostaglandin-E synthase 1 (PTGES), salivary lipocalin (SAL1) and glutathione peroxidase 1 (GPX1) were upregulated in the subcutaneous adipose tissue of Ningxiang pigs compared to that of Berkshire pigs, while cytochrome P450 family 2 subfamily B22 (CYP2B22) and glutathione peroxidase 3 (GPX3) were downregulated in the subcutaneous adipose tissue of Ningxiang pigs compared to that of Berkshire pigs (Figure 12A, C). The relative content of arachidonic acid in the subcutaneous adipose tissue of Ningxiang pigs was lower compared to that of Berkshire pigs. In the group of Berkshire pigs vs. F₁ pigs, kO00590 ‘arachidonic acid metabolism’ pathway related DEGs, glutathione peroxidase 1 (GPX1), prostaglandin-endoperoxide synthase 2 (PTGS2), phospholipase A2 group IVA (PLA2G4A) were upregulated (Figure 12B, D). The correlation analysis of SCLs and DEGs in KO00590 pathway suggested that the expression of PTGES and LOC110255237 were positive correlated with the content of arachidonic acid and some phosphatidylcholines, and the expression of GPX1 was negative correlated with the content of arachidonic acid and some phosphatidylcholines in the group of Berkshire pigs vs. Ningxiang pigs (Figure 12E). the expression of GPX1 was negative correlated with the content of some phosphatidylcholines and the expression of PLA2G4A was negative correlated with the content of one phosphatidylcholine in the group of Berkshire pigs vs. F₁ pigs (Figure 12F).

The result indicated that the content of palmitoleic acid, arachidonic acid, linoleic acid, and icosadienoic acid in the subcutaneous adipose tissue of Ningxiang pig was lower than that of Berkshire pigs and the key genes of the KEGG pathways belonging to fatty acids biosynthesis and metabolism exhibit difference expression in the subcutaneous adipose tissue across the groups of Berkshire pigs vs. Ningxiang pigs and Berkshire pigs vs. F₁ pigs. The identified these key genes will contribute to the improvement the content of palmitoleic acid, arachidonic acid, linoleic acid, icosadienoic acid for Ningxiang pig through molecular design breeding.

4. Discussion

4.1. Huge Difference Existed in the Lipid Composition of the Subcutaneous Adipose Tissue among Ningxiang pigs, Berkshire pigs and F₁ Pigs

A total of 1293 lipids were detected. 32.48% of detected lipids were TG, 6.03% of detected lipids were DG, 6.19% of detected lipids were Cer-AS, 6.11% of detected lipids were PE, 5.88% of detected lipids were PC, 5.34% of detected lipids were PS. Compared with Berkshire pigs, Ningxiang pigs had 412 downregulated and 11 upregulated lipids, interestingly, Cer(d18:0/18:0) was highly significantly regulated, while linolenic acid, pentadecanoic acid were highly significantly downregulated (Figure S6). Correlation analysis was performed on the SCLs, mapping the top 50 VIP values, in the subcutaneous adipose tissue of the group Berkshire pigs vs. Ningxiang pigs, and the majority lipids were positively correlated, while Cer(d18:0/18:0) and TG(18:0_16:1_24:6) were negatively correlated (Figure S7). F₁ pigs had 251 downregulated and 1 SP upregulated lipids in the subcutaneous adipose tissue compared to that of Berkshire pigs, Cer(d18:0/18:0) was highly significantly upregulated, while linolenic acid, pentadecanoic acid were not highly significantly downregulated (Figure S8A). For the groups of Ningxiang pigs vs. F₁ pigs, 50 SCLs were screened. Compared with Ningxiang pigs, F₁ pigs had just 7 downregulated and 43 upregulated lipids in the subcutaneous adipose tissue, and linolenic acid, pentadecanoic acid were highly significantly downregulated, while Cer(d18:0/18:0) not highly significantly upregulated (Figure S8B). The relative content of linolenic acid, pentadecanoic acid were highly significantly regulated in the subcutaneous adipose tissue of F₁ pigs compared to that of Ningxiang pigs.

4.2. The KEGG Pathways Belonging to Fatty Acids Biosynthesis and Metabolism Related DEGs and SCLs Contribute to the Difference of Fatty Acids Composition of Adipose Tissue and Fat Deposition among Ningxiang, Berkshire and F₁ Pigs

The Ningxiang pig is a famous indigenous pig breed in China with excellent meat quality. Berkshire pigs is a typical lean pig breed, famous for its high lean meat rate and fast growth rate. Besides, Berkshire pigs has been under intensive artificial selective and genetic improvement of important traits, such as growth rate, lean meat rate and backfat thickness [16]. Healthy diet attracts attention, the consumer is prefer to increase their PUFAs intake and reduce their SFAs [17]. Thus, we focus on backfat thickness and fatty acids composition of subcutaneous adipose tissue. Muti-omics association analysis for local pig breeds and modern pig breeds may screen the candidate genes for backfat thickness and fatty acids composition of subcutaneous adipose tissue, and these genes was probably selected by mankind in pigs breeding. In this study, a total of 40 fatty acyls (14 SFAs, 13 MUFAs, 11 PUFAs) were detected. Two SFAs including C(14:0), pentadecanoic acid, 3 MUFAs including 9,10-DiHOME, palmitoleic acid, heptadenoic acid, and 5 PUFAs including linoleic acid, linolenic acid, icosadienoic acid, arachidonic acid, docosapentaenoic acid were downregulated in the subcutaneous adipose tissue of the group Berkshire pigs vs. Ningxiang pigs. Besides, 9,10-DiHOME, linoleic acid and docosapentaenoic acid were downregulated in the subcutaneous adipose tissue of the group Berkshire pigs vs. F₁ pigs. In addition, C(14:0), pentadecanoic acid and linolenic acid was upregulated in the subcutaneous adipose tissue of the group Ningxiang pigs vs. F₁ pigs. Which suggested that the relative content of MUFAs and PUFAs in subcutaneous fat of Berkshire pigs was higher compared to that of Ningxiang pigs or F₁ pigs, and the relative content of C(14:0), pentadecanoic acid and linolenic acid of F₁ pigs was higher compared to that of Ningxiang pigs. It was reported that the adipose tissue of local pig breed contained lower proportion of PUFAs compared to that of modern pigs [4,18].Which was consistent with our study. Interestingly, these MUFAs (palmitoleic acid) and PUFAs (linolenic acid, linoleic acid, arachidonic acid, and icosadienoic acid) involved in multiple KEGG pathways belonging to fatty acids biosynthesis and metabolism, such as KO00061 ‘fatty acid biosynthesis’ pathway, kO01040 ‘biosynthesis of unsaturated fatty acids’ pathway, KO00590 ‘arachidonic acid metabolism’ pathway, KO00591‘linoleic acid metabolism’ pathway. Moreover, KO00061 ‘fatty acid biosynthesis’ pathway related DEGs, FASN, ACACA, CBR4, and MCH, kO01040 ‘biosynthesis of unsaturated fatty acids’ pathway related DEGs, ELOV6, SCD, and HACD, KO00590 ‘arachidonic acid metabolism’ pathway related DEGs, GGT5, PTGES, SAL1, GPX1, GPX3, CYP2B22, PTGS2, PLA2G4A and KO00591 ‘linoleic acid metabolism’ pathway related DEGs, CYP3A46, may play vital role in fatty acid composition of adipose tissue and the content of MUFAs and PUFAs. FASN is a key enzyme of fatty acids de novo synthesis pathway, ACACA catalyze irreversible formation of malonyl-CoA from acetyl-CoA and acetyl-CoA was the precursor of de novo synthesized fatty acids, which played important role in fatty acid composition of meat [19]. ELOV6, FASN, SCD were identified as candidate genes for fatty acid composition in backfat by genome-wide association studies (GWAS) of multiple pig populations[20,21]. In the Alentejano fat breed, the acetyl-CoA carboxylase activity in the dorsal subcutaneous fat was three- and ninefold higher than that in the Large White lean breed [22], which was consistent with the result of this study. Besides, several comparative transcriptomic studies indicated that FASN and ACACA were upregulated in local pig breeds compared to modern pig breeds [23,24,25]. FASN and ACACA play vital role in fatty acid biosynthesis and lipid metabolism. Moreover, it was reported that the content of palmitoleic acid was closely associated with intramuscular fat, meat color, PH24h, and palmitoleic acid might be promising indicators for superior meat quality [26]. In this study, the relative content of palmitoleic acid in the subcutaneous adipose tissue of Ningxiang pigs was lower compared to that of Berkshire pigs. ELOV6 catalyzes fatty acids elongation [27], and SCD is essential for fatty acid desaturation [28], both ELOV6 and SCD plays important role in MUFAs synthesis. It was reported that Meishan pigs contain lower activity of SCD in the subcutaneous adipose tissue compared to that of Large White pigs. Which was also consistent with our study. In this study, SCD were upregulated in the subcutaneous adipose tissue of the group Berkshire pigs vs. Ningxiang pigs. These DEGs and SCLs contribute to the difference of fatty acids composition of adipose tissue and fat deposition among Ningxiang, Berkshire and F₁ pigs.

4.3. Some of the KEGG Pathways Related genes screened by Muti-Omics Association Analysis were Probably Overlapped with Artificial Selected Genes

Muti-omics association analysis was performed to screen the key genes for backfat thickness and fatty acids composition of subcutaneous adipose tissue, and these key genes was probably selected by mankind in pigs breeding. The analysis revealed 19, 16, and 7 KEGG pathways (P < 0.05) enriched in the groups of Berkshire pigs vs. Ningxiang pigs, Berkshire pigs vs. F₁ pigs, and Ningxiang pigs vs. F₁ pigs, respectively, based on transcriptomics data (Figure 7B). Additionally, 52, 30, and 35 KEGG pathways (P < 0.05) were enriched in the groups of Berkshire pigs vs. Ningxiang pigs, Berkshire pigs vs. F₁ pigs, and Ningxiang pigs vs. F₁ pigs, respectively, based on transcriptomics data (Figure 4). The comparative genomic analysis between local pig breeds and modern pig breeds uncover the selective regions and KEGG pathways related to the candidate genes in those regions [29,30,31,32].A comparative genomic analysis indicated that the KEGG pathways ‘phosphatidylinositol signaling system’, ‘inositol phosphate metabolism’ related to the candidate genes in selective regions from Fst statistics were significantly enriched in the group of Duroc vs. Xiang pigs, and the KEGG pathway ‘thermogenesis’ related to the candidate genes in selective regions from Fst statistics was enriched in the group of Yorkshire pigs vs. Xiang pigs [33]. Another study found that the KEGG pathways ‘carbohydrate digestion and absorption’ related to the candidate genes in selective regions from Fst statistics were significantly enriched in the group of Anqing six-end-white pig (ASP) vs. Asian wild boars [34].Interestingly, the result of our study suggested that KO04070 ‘phosphatidylinositol signaling system’ pathway, KO00562 ‘inositol phosphate metabolism’ pathway and KO04714 ‘Thermogenesis’ pathway was significantly enriched in the groups of Berkshire pigs vs. Ningxiang pigs and Berkshire pigs vs. F₁ pigs. KO04973 ‘carbohydrate digestion and absorption’ pathway was significantly enriched in the groups of Berkshire pigs vs. Ningxiang pig. Which suggested that these KEGG pathways related genes were probably under intensive artificial selective. For instance, FASN and ACACA were the DEGs of KO00061 ‘fatty acid biosynthesis’ pathway, which was significantly enriched in the groups of Berkshire pigs vs. Ningxiang pigs in this study. Asian allele of FASN was thought to have an important effect on backfat phenotypes and under strong selection in modern European pig breeds [35]. ACACA located selected regions detected by CLR test in Chinese Landrace and Yorkshire Pigs [36].

4.4. Application of Multi-Omics Data to Pig Breeding for Improvement the Backfat Thickness and Fatty Acids Composition of Adipose Tissue

With the dramatically growth of multi-omics data, its integration analysis can reveal candidate genes for complex trait and gene regulatory network. However, it was still a challenge to integrate the multi-omics data [37]. Identification functional variants in candidate genes, which was screened through integrate multi-omics data will provided more useful new functional molecular markers for pigs breeding [19,38,39]. Then, intelligent algorithms were used to optimize the prediction effect of these functional molecular markers. Thus, mining candidate genes using multi-omics data was an important step. In this study, multiple KEGG pathways related genes were associated with backfat thickness and fatty acids composition of adipose tissue, some genes such as FASN, SCD, ACACA, ELOV6, which was considered as important functional genes for backfat thickness and fatty acids composition of adipose tissue. We had validated the credibility of our research results in multiple dimensions. These candidate genes will be used to identify functional variants and develop functional molecular markers. We provided a possible method to apply to multi-omics data to pig breeding.

5. Conclusions

Ningxiang pig, acted as male parent, crossed to Berkshire pigs, the meat quality of their offspring (F₁) pigs were superior. In order to study the molecular mechanism of fat deposition difference among Ningxiang, Berkshire and F₁ pigs. Transcriptomic and lipidomic analysis were conducted. The functional enrichment analysis indicated that the pathway of the lycine, serine and threonine metabolism, arachidonic acid metabolism, fatty acid biosynthesis, linoleic acid metabolism, alpha-linolenic acid metabolism, glycerophospholipid metabolism were significantly enrichment by comparing the groups of Ningxiang pigs vs. Berkshire pigs and F₁ pigs vs. Berkshire pigs. The fat deposition difference among Ningxiang, Berkshire and F₁ pigs may be caused by the difference of the expression pattern of key genes in multiple enrichment KEGG pathways. We found that the DEGs (PHGDH, AOC2) and the SCLs (Phosphatidylserines) were significantly correlated and play vital role to lycine, serine and threonine metabolism. FASN, ACACA, CBR4, SCD, ELOV6, HACD, CYP3A, CYP2B, gpx belonging to fatty acid biosynthesis and metabolism were identified as key genes to affect the content of palmitoleic acid, linoleic acid, arachidonic acid, icosadienoic acid in the subcutaneous adipose tissue. This study screened key genes and lipids, and investigated the molecular mechanism of fat deposition difference among Ningxiang, Berkshire and F₁ pigs.