Molecular Regulation of Differential Lipid Molecule Accumulation in the Intramuscular Fat and Abdominal Fat of Chickens

Reducing abdominal fat (AF) accumulation and increasing the level of intramuscular fat (IMF) simultaneously is a major breeding goal in the poultry industry. To explore the different molecular mechanisms underlying AF and IMF, gene expression profiles in the breast muscle (BM) and AF from three chicken breeds were analyzed. A total of 4737 shared DEGs were identified between BM and AF, of which 2602 DEGs were upregulated and 2135 DEGs were downregulated in the BM groups compared with the AF groups. DEGs involved in glycerophospholipid metabolism and glycerolipid metabolism were potential regulators, resulting in the difference in lipid metabolite accumulation between IMF and AF. The PPAR signaling pathway was the most important pathway involved in tissue-specific lipid deposition. Correlation analysis showed that most representative DEGs enriched in the PPAR signaling pathway, such as FABP5, PPARG, ACOX1, and GK2, were negatively correlated with PUFA-enriched glycerophospholipid molecules. Most DEGs related to glycerophospholipid metabolism, such as GPD2, GPD1, PEMT, CRLS1, and GBGT1, were positively correlated with glycerophospholipid molecules, especially DHA- and arachidonic acid (ARA)-containing glycerophospholipid molecules. This study elucidated the molecular mechanism underlying tissue-specific lipid deposition and poultry meat quality.


Introduction
Lipids are one of the major components of the chicken body. Intensive selection for a rapid growth rate leads to excessive abdominal fat (AF) accumulation [1], which exerts a negative impact on consumer acceptance and health [2]. However, intramuscular fat (IMF) is a major component that produces meat flavor, and the presence of intramuscular fat improves meat tenderness and juiciness [3][4][5]. To meet consumers' needs, reducing AF accumulation and increasing the levels of IMF simultaneously has become a major breeding goal in the poultry industry.
A previous study verified that a desirable broiler line with higher IMF and lower AF could be obtained using genetic selection [6]. AF-derived preadipocytes have higher adipogenic differentiation ability than IMF-derived preadipocytes [7]. Large transcriptomic differences were identified between IMF-and AF-derived preadipocyte differentiation [7]. These results indicated that the deposition of AF and IMF are subject to different regulatory mechanisms in chickens. Although transcriptomics has been applied to reveal the molecular mechanisms underlying AF and IMF deposition between different chicken breeds [8,9], it is more effective to find a balanced selection between AF and IMF, rather than selection for into two separate mixtures, respectively. A total of 18 mixture samples (JYBM1, JYBM2,  JYBM3, JYAF1, JYAF2, JYAF3, GYBM1, GYBM2, GYBM3, GYAF1, GYAF2, GYAF3, TCBM1,  TCBM2, TCBM3, TCAF1, TCAF2, TCAF3) were obtained for further analysis. The library preparations were sequenced using an Illumina Novaseq platform, and 150 bp paired-end reads were generated. Adapter sequences and other low-quality data were removed using Trimmomatic v0.32. Sequencing reads were then aligned to the chicken reference genome of GRCg6a using HISAT2 v2.0.5.

Differential Expression Analysis and KEGG Enrichment Analysis
Feature Counts v1.5.0-p3 was used to count the read numbers mapped to each gene. Then, the FPKM of each gene was calculated based on the length of a gene and the read count mapped to that gene. Differential expression of transcripts was performed using the DESeq2 R package (v 1.20.2). The resulting p-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate (FDR). Genes with an adjusted p (Padj)-value < 0.05 and |log2(Fold Change)| ≥ 1 were considered DEGs. A Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using the KEGG online website (http://www.genome.jp/kegg/ (accessed on 6 November 2022)). We used the clusterProfiler R package to test the statistical enrichment of differentially expressed genes in the KEGG pathway. The pathways with an adjusted p-value less than 0.05 were considered significantly enriched.

cDNA Synthesis and Real-Time qRT-PCR Analysis
The PrimeScript™ RT Reagent Kit was used for the reverse transcription of total RNA. Primers used for quantitative real-time (qRT-PCR) were listed in Table 1. The qRT-PCR was conducted with a CFX-96 (Bio-Rad, Inc., Richmond, CA, USA) qRT-PCR system using a 10 µL reaction volume with 1 µL of cDNA, 0.5 µL of the forward and reverse primers (10 µM) for each gene, 5 µL of TB Green Premix Ex Taq (Tli RNase H Plus) (TaKaRa), and 3 µL of double-distilled H 2 O. All mRNA expression levels were normalized to the GAPDH mRNA level. The relative expressions of related genes were calculated using the 2 −∆∆Ct method, and three biological replicates were performed on each sample. Pearson's correlation coefficients were generated to evaluate the accuracy of RNA-Seq.
The raw data files were processed using Compound Discoverer 3.1 (Thermo Fisher) to perform peak alignment, peak picking, and quantitation of each metabolite. Then, the peak intensities were normalized to total spectral intensity according to the set parameters. The main parameters were set as follows: retention time tolerance, 0.2 min; actual mass tolerance, 5 ppm; signal intensity tolerance, 30%; signal/noise ratio, 3; and minimum intensity, 100,000. The normalized data were used to predict the molecular formula based on additive ions, molecular ion peaks, and fragment ions. Then, the peaks were matched with Lipidmaps and Lipidblast databases to obtain accurate qualitative and relative quantitative results.

Integrated Analysis of Transcriptomics and Lipidomics
In our recently published study, we used the same samples for lipidomic analysis. A total of 209 lipid molecules were commonly determined as the maker candidates for AF and IMF among the three chicken breeds. Based on those results, 3 breast muscle samples and 3 abdominal fat samples were pooled into two separate mixtures, respectively, when we conducted transcriptomics analysis. We also calculated the mean of the relative abundance for the corresponding samples (3 breast muscle or 3 abdominal fat samples) of lipid molecules to guarantee that the absolute same samples were used for integrated analysis of transcriptomics and lipidomics. Pair-wise Pearson's correlation coefficients were generated to evaluate the correlation between differential metabolites and DEGs between AF and BM in Guangyuan grey chickens, Jiuyuan black chickens, and Tibetan chickens, respectively. Correlation coefficients were visualized with heatmap plots using the R package "complex_heat_map".

Overall Description of Sequencing Data
Overall, we obtained 118.9 Gb clean bases and 792,671,934 clean reads from RNA-seq data. These clean reads were uniquely mapped to the chicken genome (GRCg6a), and the mapping frequencies were found to vary from 86.91% to 91.27% (Table S2). Among the mapped reads, an average of 92.95% of the total mapped reads was mapped to exons, 3.30% was mapped to introns, and 3.75% was mapped to the intergenic regions.

KEGG Enrichment Analysis of DEGs Involved in Lipid Metabolism
To eliminate the effect of breed on lipid metabolism, we used 4737 common DEGs simultaneously identified in the three chicken breeds for functional enrichment analysis. Here, we display the top 20 pathways for the GYBM vs. GYAF, JYBM vs. JYAF, and TCBM vs. TCAF comparisons. Based on 4737 shared DEGs, Figure 2A-C shows the results for KEGG enrichment in Guangyuan grey chickens, Jiuyuan black chickens, and Tibetan chickens, respectively. Interestingly, we found that these shared DEGs identified in the GYBM vs. GYAF, JYBM vs. JYAF, and TCBM vs. TCAF comparisons were all enriched in 15 KEGG pathways (p < 0.05), including carbon metabolism; biosynthesis of amino acids; the citrate cycle (TCA cycle); pyruvate metabolism; glycine, serine, and threonine metabolism; 2-oxocarboxylic acid metabolism; glycolysis/gluconeogenesis; propanoate metabolism; the calcium signaling pathway; the pentose phosphate pathway; glyoxylate and dicarboxylate metabolism; fructose and mannose metabolism; the apelin signaling pathway; the peroxisome proliferator-activated receptor (PPAR) signaling pathway; and vascular smooth muscle contraction (p < 0.05). Thirty-one shared DEGs were significantly enriched in the PPAR signaling pathway, which was the sole pathway involved in lipid metabolism ( Figure 2D). In Guangyuan grey chickens and Tibetan chickens, 28 DEGs were significantly enriched in the adipocytokine signaling pathway (p < 0.05), while no significant difference was observed in Jiuyuan black chickens (p > 0.05). the peroxisome proliferator-activated receptor (PPAR) signaling pathway; and vascular smooth muscle contraction (p < 0.05). Thirty-one shared DEGs were significantly enriched in the PPAR signaling pathway, which was the sole pathway involved in lipid metabolism ( Figure 2D). In Guangyuan grey chickens and Tibetan chickens, 28 DEGs were significantly enriched in the adipocytokine signaling pathway (p < 0.05), while no significant difference was observed in Jiuyuan black chickens (p > 0.05).
PEER REVIEW 9 of 17 The node color indicates the expression of genes: (red) up-regulated and (green) down-regulated in the BM groups relative to the AF groups. GeneRatio: the ratio of the number of differential genes annotated to the KEGG pathway number to the total number of differential genes.

Gene Expression Validation of DEGs using qRT-PCR
The results from the RNA-Seq analysis were further confirmed by validating the expression data for six randomly selected DEGs. The log2 (fold change) of RNA-Seq data for each gene was the average value, so we also calculated the average value of the qRT-PCR result. The accuracy of the RNA-Seq was authenticated by the consistency in the qRT-PCR results using the RNA-Seq data with a correlation coefficient of 0.94 (Figure 3). The node color indicates the expression of genes: (red) up-regulated and (green) down-regulated in the BM groups relative to the AF groups. GeneRatio: the ratio of the number of differential genes annotated to the KEGG pathway number to the total number of differential genes.

Gene Expression Validation of DEGs using qRT-PCR
The results from the RNA-Seq analysis were further confirmed by validating the expression data for six randomly selected DEGs. The log2 (fold change) of RNA-Seq data for each gene was the average value, so we also calculated the average value of the qRT-PCR result. The accuracy of the RNA-Seq was authenticated by the consistency in the qRT-PCR results using the RNA-Seq data with a correlation coefficient of 0.94 (Figure 3).
The results from the RNA-Seq analysis were further confirmed by validating the expression data for six randomly selected DEGs. The log2 (fold change) of RNA-Seq data for each gene was the average value, so we also calculated the average value of the qRT-PCR result. The accuracy of the RNA-Seq was authenticated by the consistency in the qRT-PCR results using the RNA-Seq data with a correlation coefficient of 0.94 (Figure 3).

Correlations between Shared DEGs and Differential Lipid Molecules
In our recently published study, the same samples were used for lipidomic analysis. We identified a large number of shared glycerophospholipid lipid molecules that were significantly upregulated in IMF compared with those in AF. These glycerophospholipid lipid molecules included 11 cardiolipin (CL), 1 phosphatidic acid (PA), 33 phosphatidylcholines (PC), 19 phosphatidylethanolamines (PE), 4 phosphatidylinositols (PI), and 10 phosphatidylserines (PS). The difference between these individual glycerophospholipid

Correlations between Shared DEGs and Differential Lipid Molecules
In our recently published study, the same samples were used for lipidomic analysis. We identified a large number of shared glycerophospholipid lipid molecules that were significantly upregulated in IMF compared with those in AF. These glycerophospholipid lipid molecules included 11 cardiolipin (CL), 1 phosphatidic acid (PA), 33 phosphatidylcholines (PC), 19 phosphatidylethanolamines (PE), 4 phosphatidylinositols (PI), and 10 phosphatidylserines (PS). The difference between these individual glycerophospholipid compounds was the root cause of the higher content of phospholipids in intramuscular fat. In the present study, the shared DEGs related to glycerophospholipid metabolism and the PPAR signaling pathway were screened. It is reasonable to speculate that these DEGs were involved in the biosynthesis of glycerophospholipid molecules in BM. Therefore, an integrated analysis of transcriptomics and lipidomics was conducted to evaluate the correlation between differential glycerophospholipid metabolites and DEGs related to glycerophospholipid metabolism and the PPAR signaling pathway. In total, 59 shared DEGs (31 related to the PPAR signaling pathway and 28 related to glycerophospholipid metabolism) and 78 shared differential glycerophospholipid metabolites were used in the correlation analysis for Guangyuan grey chickens, Jiuyuan black chickens, and Tibetan chickens, respectively. The heatmap plots are presented in Figure 4. Each row represents a glycerophospholipid molecule, and each line represents a gene. A total of 1896, 3034, and 1496 significant correlations between DEGs and differential metabolites were identified in Guangyuan grey chickens, Jiuyuan black chickens, and Tibetan chickens, respectively. In order to identify the candidate genes for specific lipid molecule accumulation more accurately, shared significant correlations among the three chicken breeds were screened. A total of 777 significant shared correlations were detected. Considering the important role of PUFA in meat flavor and nutritional value, we concentrated on the PUFA-enriched glycerophospholipid molecules ( Figure 5). We found most representative DEGs enriched in the PPAR signaling pathways were negatively correlated with PUFA-enriched glycerophospholipid molecules (p < 0.01).      In our previous study, we identified that two triacylglycerol metabolites (TAG (16:1-18:1-18:1) and TG (16:1(9Z)/16:1(9Z)/18:1(9Z))[iso3])) were significantly higher in AF tissue, which might be a possible reason for the higher triglyceride content verified in AF. To decrease AF content, we also conducted a correlation analysis between two triacylglycerol molecules and the shared DEGs involved in glycerolipid metabolism and the PPAR signaling pathway. Unfortunately, no shared correlations were identified among the three chicken breeds (p > 0.01).

Discussion
In poultry, excessive AF accumulation causes low slaughter efficiency and exerts a negative impact on consumer acceptance and health [2]. Increased IMF content contributes to better meat quality, including tenderness, flavor, and nutritional value [4].
Lowering the AF content and increasing the IMF content can effectively increase the economic value of broilers, which has become a major breeding goal in the poultry industry. Most previous studies concentrated on exploring the potential gene markers regulating AF or IMF deposition unilaterally based on a comparison of two types of chicken breeds whose AF deposition or IMF deposition was exceptionally varied. For example, after obtaining gene expression profiles of breast muscles for Beijing-You and AA chickens, researchers identified candidate genes related to IMF deposition [8]. A previous study also used extremely high-and low-abdominal fat chicken groups to investigate crucial genes related to AF deposition [16]. Although the genetic mechanisms underlying chicken fat deposition have been widely studied, few studies were conducted to determine the mechanism that leads to tissue-specific lipid molecule accumulation. Nowadays, indigenous chicken breeds are preferred by customers due to their better meat quality. Previous studies indicated that the intramuscular fat content in these breeds is higher than that in commercial chicken breeds [8]. In this study, we selected three indigenous chicken breeds for integrative analysis of transcriptome and previous lipidomics data.
Unlike the marbling distribution of IMF in domestic animals, the IMF of chickens cannot be separated visually. To elucidate the regulatory mechanism underlying tissuespecific lipid metabolism, DEGs were identified between AF and BM. Those genes related to glycerolipid metabolism, glycerophospholipid metabolism, and sphingolipid metabolism were displayed. To eliminate the effects of breed, three chicken breeds (the Guangyuan grey chicken, Jiuyuan black chicken, and Tibetan chicken) were selected in this study. Only the shared DEGs identified in GYBM vs. GYAF, JYBM vs. JYAF, and TCBM vs. TCAF simultaneously were used for subsequent analysis. A total of 4737 shared DEGs showed a completely consistent trend between BM and AF, of which 2602 DEGs were upregulated and 2135 DEGs were downregulated in the BM group compared with the AF group. We found that more shared DEGs related to glycerophospholipid metabolism were upregulated in the BM groups, while the expressions of more DEGs related to glycerolipid metabolism were higher in the AF groups. Interestingly, a previous study indicated that phospholipids were found in a considerable proportion of BM [11], while the triglyceride content in AF was significantly higher than that in BM [17]. Thus, these DEGs were potential regulators causing the difference in triglyceride and phospholipid content between AF and BM. It is well known that the AF content is far greater than the IMF content in chickens. At the cellular level, researchers have indicated that the lipogenesis of preadipocytes derived from AF was significantly increased compared with that of preadipocytes derived from BM [18]. Consistent with previous results in vitro, the downregulated DEGs related to glycerolipid metabolism were responsible for the decrease in total lipids in BM.
Based on 4737 shared DEGs, the results of the KEGG signaling pathway analysis showed that DEGs between BM and AF were jointly enriched in 15 pathways. Thirty-one shared DEGs were significantly enriched in the PPAR signaling pathway with the classic mediation of lipid metabolism. Among these 31 DEGs, the RNA-seq data showed that the PLIN1 , FABP4, SCD, PPARG, GK2, CD36, APOC3, FABP7, LPL, RXRG, SLC27A6, FABP3,  ACSL1, DBI, ACOX1, PLTP, APOA1, ANGPTL4, FABP5, ACSL4, ILK, ACAA1, PCK2,  and SLC27A2 genes were downregulated, while only 7 genes (CYP27A1, HMGCS1, ME1,  LOC121106914, RXRA, ACOX2, and FABP1) were upregulated in the BM group. Many studies have confirmed the role of these genes in lipid metabolism [19][20][21]. For example, PLIN1, which was reported to protect lipid droplets from the hydrolytic activity of hormonesensitive lipase [18], had significantly higher expression levels in AF. PPARG, the most adipocyte-specific member of the PPAR family, is mainly expressed in adipose tissue [20]. Fatty acid-binding proteins (FABPs), also known as intracellular lipid chaperones, are a family that regulates lipid trafficking and responses in cells [19]. Fatty acids cross the cell membrane via a protein-mediated mechanism, and FABPs are responsible for fatty acid uptake [22]. Different isoforms of the FABP family are uniquely expressed in tissues involved in lipid metabolism [23]. The mRNA expressions of heart FABP (H-FABP/FABP3), adipocyte FABP (A-FABP/FABP4), and epidermal FABP (E-FABP/FABP5), brain FABP (B-FABP/FABP7) were significantly higher in AF, while the expression of liver FABP (L-FABP/FABP1) was significantly higher in IMF. This indicated that FABPs may provide tissue-specific control of lipid metabolism. Consistent with previous results in vitro [18], it was considered that tissue-specific lipid deposition mostly depended on the changed expression of related genes through the PPAR pathway.
Glycerophospholipids are dominant in cell membranes and play a vital role in cellular functions such as signal transduction, protein function, and regulation of transport processes [24]. In our recently published study, lipidomic analysis showed that a significantly higher abundance of 11 CL, 1 PA, 33 PC, 19 PE, 4 PI, and 10 PS was observed in IMF compared with that in AF. PUFAs are predominantly deposited in the glycerophospholipid of muscle rather than free fatty acids [25,26]. The higher abundance of PUFA-containing glycerophospholipids in IMF could elucidate its beneficial roles in imparting a characteristic flavor and nutritional value to meat [12]. Therefore, exploring the important candidate gene markers for specific PUFA-containing glycerophospholipid molecules would contribute to a more thorough understanding of IMF deposition and poultry meat quality.
The limited amounts of n-3 PUFAs present in meats are nutritionally important [27]. It is worthwhile to mention the beneficial function of n-3 PUFAs, notably EPA and DHA. They play a vital role in the prevention and treatment of certain diseases, such as Alzheimer's disease, cardiovascular disease (CVD), type 2 diabetes, hypertension, psychiatric diseases, and several cancers [28,29]. Furthermore, DHA is the core nutrient for human brain function, neuronal development, and visual acuity [30]. In the present study, significantly negative correlations were observed between four genes (FABP5, PPARG, ACOX1, and GK2) and seven glycerophospholipid molecules carrying DHA and EPA (PC . However, previous studies have shown that DHA and EPA could be natural PPAR agonists to activate PPARG signaling [31,32]. In golden pompano hepatocytes, the overexpression of PPARG increased the DHA content, whereas the suppression of PPARG expression diminished this effect, suggesting that PPARG plays an active role in regulating DHA content [33]. Here, we conducted a correlation analysis between PPARG and individual DHA-and EPA-containing glycerophospholipid rather than total DHA or EPA content. The molecular mechanism underlying the PPARG regulation of these seven glycerophospholipid molecules carrying DHA and EPA needs to be verified in chicken intramuscular preadipocytes. Similar to DHA, ARA also has a significant physiological function during development and growth, especially for the central nervous system and retina [34]. In this study, GPD2, GPD1, PEMT, CRLS1, and GBGT1 were collectively involved in the positive regulation of ARA-enriched metabolites. Whether these DEGs positively regulated ARA content or these ARA-enriched molecules also needs to be verified at the cellular level.

Conclusions
In this study, an analysis of transcriptome data identified 4737 shared differentially expressed genes between IMF and AF deposition among three chicken breeds. Differentially expressed genes involved in glycerophospholipid metabolism and glycerolipid metabolism were potential regulators resulting in the difference between lipid metabolite accumulation in IMF and AF. The PPAR signaling pathway was the most important pathway involved in tissue-specific lipid deposition. Most representative DEGs enriched in the PPAR signaling pathways, such as FABP5, PPARG, ACOX1, and GK2, were negatively correlated with PUFA-enriched glycerophospholipid molecules. Most DEGs related to glycerophospholipid metabolism, such as GPD2, GPD1, PEMT, CRLS1, and GBGT1, were positively correlated with glycerophospholipid molecules, especially DHA-and ARA-containing glycerophospholipid molecules. These results contributed to a more thorough understanding of lipid deposition and poultry meat quality.