Differential regulation of intramuscular fat and abdominal fat deposition in chickens

Chicken intramuscular fat (IMF) content is closely related to meat quality and performance, such as tenderness and flavor. Abdominal fat (AF) in chickens is one of the main waste products at slaughter. Excessive AF reduces feed efficiency and carcass quality. To analyze the differential deposition of IMF and AF in chickens, gene expression profiles in the breast muscle (BM) and AF tissues of 18 animals were analyzed by differential expression analysis and weighted co-expression network analysis. The results showed that IMF deposition in BM was associated with pyruvate and citric acid metabolism through GAPDH, LDHA, GPX1, GBE1, and other genes. In contrast, AF deposition was related to acetyl CoA and glycerol metabolism through FABP1, ELOVL6, SCD, ADIPOQ, and other genes. Carbohydrate metabolism plays an essential role in IMF deposition, and fatty acid and glycerol metabolism regulate AF deposition. This study elucidated the molecular mechanism governing IMF and AF deposition through crucial genes and signaling pathways and provided a theoretical basis for producing high-quality broilers.


Introduction
Fat deposition in broilers is affected by genetic factors (breed, sex, and genotype) and nongenetic factors (nutrition, age, and environmental factors). The rate of fat deposition varies across tissues in broilers. In addition, glucose and lipid metabolism and hormone sensitivity vary in adipocytes from different tissues [1][2][3][4]. The deposition of abdominal fat (AF) and intramuscular fat (IMF) is regulated by different mechanisms in chickens [5], and this characteristic allows the genetic selection of broilers with high IMF and low AF [6].
Previous studies showed a significant and positive correlation between IMF and AF [7,8]. In contrast, a study found that the chicken with lower AF had higher IMF. The discrepancy may be caused by breeding characteristics, environmental factors, genetic selection, and sampling methods [6]. A balanced selection population (increasing IMF, decreasing AF percentage (AFP)) was more effective than selecting for IMF alone in Jing-xing yellow chicken [9]. Various genetic mechanisms regulate lipid deposition in tissues. Genes encoding perilipin and long-chain acyl CoA dehydrogenase are highly expressed in subcutaneous adipocytes, whereas genes encoding bone morphogenetic proteins 4 and 7 are abundantly expressed in intramuscular adipocytes [10]. With the development of cost-effective sequencing technologies, a large amount of transcriptome data has been generated under different biological back scenes. Currently, based on these highthroughput sequencing data, Weighted correlation network analysis (WGCNA) has been successfully applied to identify gene expression networks and biomarkers of interest in various biological domains [11][12][13][14][15][16]. Studies have shown that the transcription factor CREB3L1 was associated with low AF, whereas the transcription factor Open Access *Correspondence: wenjie@caas.cn † Na Luo and Jingting Shu contributed equally to this work. 1 State Key Laboratory of Animal Nutrition, Institute of Animal Sciences, Chinese Academy of Agricultural Sciences, Beijing 100193, China Full list of author information is available at the end of the article L3MBTL1 and the cofactor TNIP1 were related to high IMF and low AF in the Jing-xing Huang chicken line [17]. These results underscore the need to study the genetic mechanisms regulating fat deposition in broiler chickens. Wenchang breed is mainly distributed in Hainan, the lowest altitude in south China. It is an excellent local breed with more than 400 years of breeding history. Wenchang chicken is usually marketed in 100-110 days and weighs about 1.4-1.8 kg. Wenchang chicken has the characteristics of thin and tender skin, fragrant meat, and high IMF. However, the AF content of Wenchang chicken is also high, which will cause waste in production and reduce the uniformity of population. Therefore, it is necessary to understand the characteristics of fat deposition in Wenchang chickens, which is beneficial to improve Wenchang chickens' quality and production efficiency. Few studies jointly evaluated IMF and AF in broilers. This study evaluates differential fat deposition in breast muscle (BM) and AF tissues of Wenchang chicken, gene expression profiles these tissues based on transcriptome data, and genes and signaling pathways involved in IMF and AF deposition in chicken to reveal the difference of molecular regulation mechanisms.

Evaluation of sequencing data
Transcriptomic data were obtained from the BM and AF of 18 Wenchang chickens. Thirty-six RNA-sequencing (RNA-seq) libraries were constructed and sequenced (Additional file: Table S1). A total of 23,910 genes were detected, of which 17,179 genes were expressed in BM, and 18,912 genes were expressed in AF. In addition, 72.24% (16,753/23,910) of the genes were expressed in both BM and AF (Fig. 1). Raw sequence reads from each library are shown in Additional files: Tables S2 and S3. The most genes abundantly expressed in BM and AF were glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and fatty acid-binding protein (FABP4), respectively.
To explore differences in IMF and AF deposition (Table 1), Principal component analysis (PCA) of fatty acid composition was performed to assess differences in fat deposition in BM and AF tissue ( Fig. 2A) [18].
PCA analysis was performed on 16,753 genes coexpressed in BM and AF (Fig. 2B). PCA analysis was performed on all expressed genes in BM and AF, respectively ( Fig. 2C and 2D).

WGCNA identified module genes significantly related to phenotypes
Genes with fragments per kilobase of transcript per million fragments mapped values of less than one were excluded. Gene co-expression networks were constructed using 18,295 expressed genes from BM and 19,778 expressed genes from AF, and the coordination matrix was converted to a topology matrix. The similarity in gene co-expression between modules was quantified using a dynamic tree-cutting algorithm, and modules with height values of less than 0.25 were merged (the minimum number of genes in the module was set to 30) (Fig. 3).
Twenty and twenty-seven expression modules were detected in BM (Additional file: Table S4) and AF (Additional file: Table S5), respectively (Fig. 4), and the number of genes in each network module was counted. The largest module in BM (blue) and AF (turquoise) contained 6972 and 5755 genes, respectively, and the smallest module in BM (dark brown) and AF (white) contained 113 and 47 genes, respectively.
Gene modules significantly associated with genetic traits were identified. Absolute values of 0 and 1 indicated weak and strong correlation between the module and the trait, respectively. Among the evaluated traits, IMF was significantly and positively correlated with the blue module (r = 0.51, P = 0.04), triglycerides (TG) was significantly and positively associated with the darkred module (r = 0.72, P = 0.002), total cholesterol (TCHO) was significantly and positively related to the darkorange module (r = 0.63, P = 0.008), phospholipids (PLIP) was significantly and positively associated with the orange and blue modules (r = -0.64, P = 0.007; r = -0.51, P = 0.04) and traits such as fatty acids composition showed significant correlations with related modules. AF weight showed a significant and positive correlation with the darkgreen and lightgreen modules (r = 0.65, P = 0.009; r = 0.71, P = 0.003), and AFP showed a significant and positive association with the lightgreen and yellow modules (r = 0.64, P = 0.01; r = 0.54, P = 0.04). In AF tissue, TG was significantly and positively correlated with the darkred and turquoise modules (r = 0.55,

Identification of DEGs in BM and AF
We selected three individuals with high IMF content and three individuals with low IMF content for BM tissue's difference analysis. And we selected three individuals with high AFP and three individuals with low AFP for AF tissue's difference analysis (Additional file: Table S6). Analysis using Deseq2 detected 461 DEGs (163 upregulated and 298 downregulated) in the high and low IMF groups (Additional file: Table S7). Pathway enrichment analysis revealed that these genes were enriched in 29 pathways, mainly including metabolic pathways, MAPK signaling pathways, arachidonic acid metabolism, glycolysis/gluconeogenesis, starch, and sucrose metabolism, fructose and mannose metabolism, and other signaling pathways (Additional file: Table S8). GO analysis showed that DEGs in this group were enriched in 27 biological processes, 16 cellular components, and 12 molecular functions (Additional file: Table S9). A total of 2010 genes (1301 upregulated and 700-downregulated) were detected in the high and low AF groups (Additional file: Table S10). These genes were enriched in 59 KEGG pathways, mainly including metabolic pathways, PPAR signaling pathway, MAPK signaling pathway, tyrosine metabolism, glycerolipid metabolism, phenylalanine metabolism, glycerophospholipid metabolism, pyruvate metabolism, unsaturated fatty acid biosynthesis, fatty acid metabolism, and other signaling pathways (Additional file: Table S11). GO analysis showed that DEGs in this group were enriched in 66 biological processes, 21 cellular components, and 18 molecular functions (Additional file: Table S12).

Identification of pathways in BM and AF
The genes expressed in the modules significantly associated with all phenotypes were selected from BM and AF tissues, respectively. The significantly expressed genes in the modules significantly associated with the traits were pooled, and a total of 3818 expressed genes were identified in BM tissue, while a total of 5826 expressed genes were identified in AF tissue. These genes were analyzed separately in association with differentially expressed genes, in which 114 coexpressed genes were obtained in BM tissue (Fig. 5A), while a total of 1229 co-expressed genes were obtained in AF tissue (Fig. 5B).
The differential module genes obtained from the two tissues were analyzed separately for KEGG pathway enrichment. In BM tissue, 114 genes were assigned to 5 KEGG pathways (P < 0.05) (Fig. 6A), and genes were enriched in MAPK signaling pathway, ECMreceptor interaction, Gap junction, Tight junction, Vascular smooth muscle contraction, GnRH signaling pathway, and Focal adhesion. In AF tissue, 1229 genes were assigned to 41 KEGG pathways (P < 0.05) (Fig. 6B), and genes were mainly enriched in biological processes, including Oxidative phosphorylation, ABC transporters, C-type lectin receptor signaling pathway, Phosphatidylinositol signaling system.

Discussion
Yellow feather broiler chickens are widely consumed in China, and chicken meat production has increased in recent years. IMF and AF contents strongly affect the quality of chicken meat. Fasting before slaughter is essential in the experimental sampling process [19]. The fasting method can reduce the volume of stomach contents and the risk of visceral microbial contamination during slaughter [20]. In general, moderate fasting (no more than 24 h) is suitable for meat quality because it accelerates glycogen consumption [21,22]. Especially in poultry, fasting can decrease breast muscle glycogen and increase initial pH [23][24][25]. Studies have shown that fasting for a certain period before slaughter can control the muscle glycogen content of beef cattle without affecting the marbling score [26][27][28]. Studies have shown that fasting for 12 h has less influence on the structural changes and plasma corticosterone level and less on the energy metabolism level of chest muscle tissue [29]. Therefore, fasting for 12 h before slaughter would not significantly impact subsequent meat quality and energy metabolism levels, and the effect of glucose metabolism in food on the chicken could be ruled out.
Although the mechanisms underlying fat deposition in broiler chickens are well known, differences in the rate of fat deposition in BM and AF tissues have not been determined.
PCA results based on phenotypes showed a strong separation of breast muscle tissue from abdominal fat tissue, as expected. This suggests that BM tissue and AF tissue have different regulatory roles in phenotype. Many DEGs can be identified by RNA-seq analysis, however, the characterization of gene expression patterns and correlation with phenotypes are challenging, underscoring the need to perform weighted co-expression network analyses. In the study, DEGs were identified in BM and AF by RNA-seq analysis and joint analysis with WGCNA to identify specifically expressed genes and signaling pathways for IMF and AF deposition. The molecular regulatory mechanisms of fat deposition in different tissues at the same age were elucidated in chicken.
Based on the differentially expressed genes, the results of KEGG signaling pathway analysis showed that DEGs from BM and AF tissues were jointly enriched in the AGE-RAGE signaling pathway, amino acid biosynthesis, cytokine-cytokine receptor interactions, local adhesion, GnRH signaling pathway, MAPK signaling pathway, metabolic pathway, in diabetic complications. Tight junctions and other signaling pathways. Studies have shown that MAPK signaling pathway is involved in lipid deposition [30,31], which can regulate PPAR pathway [32]. Moreover, PPAR signaling pathway also has the function of regulating lipid metabolism [33,34].
In contrast, DEGs in BM tissue were specifically enriched in pathways related to gluconeogenesis (glycolysis/gluconeogenesis, starch, sucrose metabolism, fructose, mannose metabolism, pentose, and glucuronide). LDHA and LDHB in the glycolysis/gluconeogenesis signaling pathway are involved in pyruvate metabolism and tricarboxylic acid cycle, facilitating the glycolytic process by converting pyruvate to lactate [35]. In addition, GBE1 is essential in starch and glycogen formation metabolism, and there is a strong link between the hexokinase family genes HK1, HK2, HK3 in this pathway and the metabolic process of glucose. These data suggest that IMF deposition in BM tissue mainly depends on glucose metabolism, with some energy metabolism involved in the process.
WGCNA was performed using RNA-seq data to identify expressed genes in different modules and predict the role of genes in lipid deposition. Several modules were significantly associated with IMF/AFW/AFP, lipid composition, and fatty acids metabolism in BM and AF, and the expressed genes significantly associated with the phenotypes in the modules were aggregated together and analyzed jointly with DEGs, which could further identify the functional genes that play essential roles in the IMF and AF deposition. Among all expressed genes, GAPDH has the highest expression level in breast muscle tissue. GAPDH is a key enzyme in glycolysis. It plays an important role in glycolysis by catalyzing the first step of this pathway by converting D-glyceraldehyde 3-phosphate (G3P) to 3-phospho-d-glyceroylphosphate [47,48]. FABP4 gene in the lipid-binding protein family regulates fatty acid uptake and transport [41].
We identified 114 expressed genes associated with Metabolic pathways, Starch and sucrose metabolism, Intestinal immune network for IgA production, Adrenergic signaling in cardiomyocytes, and Glycolysis/Gluconeogenesis in BM. These pathways are regulated by genes such as ALDH5A1, LDHA, GAPDH, GBE1, and GPX1 and may play roles in IMF deposition. A total of 1229 expressed genes associated with 41 signaling pathways were identified in AF, of which 20 pathways that may play roles in AF deposition via PPAR signaling pathway, amino acid biosynthesis, oxidative phosphorylation, may be controlled by several genes, including ADIPOQ, ELOVL6, HMGCS2, ME3, DGKE, AOX1, UBB, SCD, APOC3, and FABP1. The results revealed that IMF deposition in BM tissue was regulated by gluconeogenesisrelated pathways (glycolysis/gluconeogenesis, starch and sucrose metabolism signaling pathways), and by several genes, including GAPDH, LDHA, GPX1, ALDH5A1, and GBE1, among which, the key differential gene ALDH5A1 can catalyze a step in the degradation of the inhibitory neurotransmitter γ-aminobutyric acid [49]. The protein encoded by LDHA is involved in pyruvate metabolism by catalyzing the conversion of lactate to pyruvate in anaerobic glycolysis. GAPDH plays an important role in glycerol metabolism. These data further suggest that IMF deposition in BM is dependent on gluconeogenesis and energy metabolism.
While expressed genes in AF tissue are enriched to PPAR signaling pathway, oxidative phosphorylation, amino acid biosynthesis, and other signaling pathways, acting through genes such as FABP1, ELOVL6, SCD, and ADIPOQ. A PPAR signaling pathway is involved in lipid droplets formation and mitochondrial metabolism, and fatty acid oxidation and lipid synthesis are necessary for cellular signaling [50]. FABP1 plays a crucial role in the PPAR pathway, PPARG expression, and fatty acid uptake, transport, and metabolism in vivo [51]. In addition, PPARϒ can be regulated by modulating SCD1 expression to control fatty acid synthesis in adipocytes [52]. ELOVLs, which encode long-chain and extra-long-chain fatty acid elongases, play an important role in synthesizing fatty acids and can limit elongation. ELOVL6 is involved in synthesizing fatty acid enzymes in vivo, promoting fatty acid elongation [53]. ACAT2 encodes acetyl-coenzyme A acetyltransferase 2, which is involved in acetyl-coenzyme A metabolism [54]. These results suggest that AF deposition may result from changes in fatty acid synthesis through mitochondrial activity.
The results suggest that different signaling pathways regulate fat deposition in BM and AF tissues. IMF deposition in BM was associated with pyruvate and citrate metabolism through GAPDH, LDHA, GPX1, GBE1, and other genes, whereas AF deposition was related to acetyl-coenzyme A and glycerophospholipid metabolism through FABP1, ELOVL6, SCD, and ADIPOQ. The transcriptional regulation of genes in network modules associated with traits and metabolic pathway analysis can provide new insights into the genetic mechanisms governing fat deposition in broiler chickens.

Conclusion
The analysis of transcriptome data from BM and AF tissues, combined with differential expression analysis and WGCNA, showed that the genes GAPDH, LDHA, GPX1, and GBE1 were involved in regulating IMF deposition, and the genes FABP1, ELOVL6, SCD, and ADIPOQ determined AF deposition. Glycolysis/gluconeogenesis and other signaling pathways play an important role in IMF deposition, whereas PPAR metabolism controls AF deposition. Therefore, we hypothesized that IMF deposition in BM tissue may affect energy metabolism during myocyte gluconeogenesis while fatty acid synthesis pathways may affect AF deposition.

Animals and Sample Collection
Wenchang chickens (18 females) of 98 days were obtained from the Institute of Poultry Research of the Chinese Academy of Agricultural Sciences (Yangzhou, China). The diet (Additional file: Table S13) was formulated according to the feeding standard (NY/T33-2004). Under the same standard conditions of light (20 lx), temperature (35℃ ~ 37℃), humidity (not less than 50%), and immunization schedule, after a 12-h fast, the chickens were stunned by electric shock and killed by cervical dislocation at 98 days. After slaughter, breast muscle (BM) and AF were collected and immediately stored at − 80 °C until use.

Measurement of Biochemical Indices [18]
The determination of fatty acids in BM and AF tissues refers to the national standard《GB/T 5413.  Food Safety National Standard for infant food and dairy products determination of fatty acids》in the first method: acetyl chloride-methanol methyl esterification method. The determination of intramuscular fat in BM tissue refers to the national standard《GB/T 5009.6-2016 National Standard for Food Safety determination of fat in food》in the first method: Soxhlet extraction. The TG, TCHO, and PLIP contents in BM and AF tissues were measured using assay kits (Nanjing Jiancheng Bioengineering Institute, Nanjing, China, and Beijing Leichuang Biotechnology Co., LTD, China). We weighed 2.0 ± 0.01 g breast muscle and abdominal fat samples and mixed them with 18 mL anhydrous ethanol in 50 ml centrifuge tubes. Each sample was homogenized by a handheld high-speed homogenizer (35,000 r/min, 15 s/ time). The samples were crushed, centrifuged at 2500 r/min for 10 min. After centrifugation, the supernatant was used for measurement. Refer to the instruction for specific sample loading quantity. After incubation, the absorbance value (510 nm) was measured using a microplate reader (Tecan Infinite 200 Pro, Switzerland).

Transcriptome Profiling
Transcriptome profiling of 18 breast muscle and 18 abdominal fat tissues was performed on the Illumina PE150 platform (Berry Genomics; Beijing, China) (https:// www. berry genom ics. com). RNA was extracted as described previously [55]. Adapter sequences and other low-quality data were removed using Cutadapt. Reads were aligned to the reference genome using HISAT, while clean reads obtained by filtering were compared to the reference genome according to the gene position information specified in the genome annotation file gtf, and the total comparison rate was approximately 90%.

Differential Expression Analysis
Based on phenotypic data, DEGs were identified in BM tissue with high IMF and low IMF (from three animals each) and AF tissue with high AF and low AF (from three animals each) using DESeq2 in R [56]. Adjusted p-values (q-values) were calculated using Benjamin and Hochberg's approach for controlling the false discovery rate. Genes with | log2(fold change) |≥ 0.58 and P < 0.05 were considered differentially expressed.

Weighted Gene Co-expression Network Analysis (WGCNA)
WGCNA was performed using the WGCNA package [57] with default settings and minor modifications. For the analysis of genes expressed in BM tissue (n = 16), the minModuleSize was set to 50, and mergeCutHeight was set to 0.25 for tissue-or stage-specific module detection (soft threshold = 4). For the analysis of genes expressed in AF tissue (n = 15), the minModuleSize was set to 30, and mergeCutHeight was set to 0.25 for tissue-or stagespecific module detection (soft threshold = 10). Dendrograms were created using topology overlap matrix, module detection, and similar module merging functions. Module-trait relationships were calculated, and correlations were adjusted using Benjamin-Hochberg correction.

Pathway Enrichment Analyses
Gene ontology (GO) enrichment analysis was performed with the DAVID database (https:// david. ncifc rf. gov/) to identify gene classes and categories. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using KOBAS version 3.0 (http:// kobas. cbi. pku. edu. cn). The significance level was set at P < 0.05. The gene-pathway interaction networks for the candidate genes were visualized with Cytoscape 3.8.0 (http:// www. cytos cape. org/) [58].

Statistical Analyses
The significance of differences between means was assessed using Student's t-test in SPSS 22.0 (IBM Corp, Armonk, NY, USA). P-values less than 0.05 were considered statistically significant. Data are presented as the mean ± standard error of the mean.