Identification of differentially expressed genes and pathways between intramuscular and abdominal fat-derived preadipocyte differentiation of chickens in vitro

Background The distribution and deposition of fat tissue in different parts of the body are the key factors affecting the carcass quality and meat flavour of chickens. Intramuscular fat (IMF) content is an important factor associated with meat quality, while abdominal fat (AbF) is regarded as one of the main factors affecting poultry slaughter efficiency. To investigate the differentially expressed genes (DEGs) and molecular regulatory mechanisms related to adipogenic differentiation between IMF- and AbF-derived preadipocytes, we analysed the mRNA expression profiles in preadipocytes (0d, Pre-) and adipocytes (10d, Ad-) from IMF and AbF of Gushi chickens. Results AbF-derived preadipocytes exhibited a higher adipogenic differentiation ability (96.4% + 0.6) than IMF-derived preadipocytes (86.0% + 0.4) (p < 0.01). By Ribo-Zero RNA sequencing, we obtained 4403 (2055 upregulated and 2348 downregulated) and 4693 (2797 upregulated and 1896 downregulated) DEGs between preadipocytes and adipocytes in the IMF and Ad groups, respectively. For IMF-derived preadipocyte differentiation, pathways related to the PPAR signalling pathway, ECM-receptor interaction and focal adhesion pathway were significantly enriched. For AbF-derived preadipocyte differentiation, the steroid biosynthesis pathways, calcium signaling pathway and ECM-receptor interaction pathway were significantly enriched. A large number of DEGs related to lipid metabolism, fatty acid metabolism and preadipocyte differentiation, such as PPARG, ACSBG2, FABP4, FASN, APOA1 and INSIG1, were identified in our study. Conclusion This study revealed large transcriptomic differences between IMF- and AbF-derived preadipocyte differentiation. A large number of DEGs and transcription factors that were closely related to fatty acid metabolism, lipid metabolism and preadipocyte differentiation were identified in the present study. Additionally, the microenvironment of IMF- and AbF-derived preadipocyte may play a significant role in adipogenic differentiation. This study provides valuable evidence to understand the molecular mechanisms underlying adipogenesis and fat deposition in chickens.


Background
Chicken is generally accepted as one of the main protein sources worldwide. In the last several decades, meat quality has decreased as a result of genetic selection for growth rate and feed conversion. Intramuscular fat (IMF) content, an important factor influencing meat quality, contributes to multiple meat quality characteristics, such as flavour, tenderness and juiciness [1][2][3]. Abdominal fat (AbF) is an important carcass trait in chickens. A higher growth rate induces larger fiber diameters and lower IMF deposition, which severely deteriorates the quality of meat [4,5]. However, the overemphasis on selection for a rapid growth rate leads to excessive fat accumulation, especially AbF accumulation [6]. Excessive fat is often discarded as waste [6][7][8]. Reducing the levels of AbF and increasing the levels of IMF have therefore become a major breeding goals in the chickens industry [6,9].
Previous studies have indicated that adipose tissues from different locations display unique physiological and biochemical characteristics [10][11][12]. In addition, glucose utilization and lipid metabolism mechanisms and hormone sensitivities are different among tissues from different locations [13][14][15][16]. IMF has specific biological features compared with fat from other locations. Previous studies have suggested that AbF has higher triglyceride (TG) levels than IMF tissue [17][18][19]. Hrdinka C et al. demonstrated that the fatty acid composition of AbF differs significantly from that of IMF [20]. Zhou et al. found that dietary supplementation with 3% conjugated linoleic acid (CLA) decreased AbF accumulation but increased IMF content [21]. Leng et al. indicated that a desirable broiler line with high IMF content but low AbF content could be obtained by genetic selection [6]. This may be because AbF deposition and IMF deposition are subject to different regulatory mechanisms. However, the mechanisms underlying regional differences in chicken adipogenesis remain unknown.
In the current study, Ribo-Zero RNA-Seq was used to systematically identify differentially expressed genes (DEGs), mRNAs (DEMs) and novel genes (DENGs) and different pathways between preadipocytes and adipocytes of IMF and AbF. These data may contribute to a more thorough understanding of tissue-specific adipogenic differentiation and poultry meat quality.

Results
Ribo-zero RNA-Seq of different chicken different adipose tissue-derived preadipocytes and adipocytes Intramuscular and abdominal preadipocytes were isolated and cultured in growth medium until they reached 80-90% confluence (Fig. 1a). Microscopy showed that the IMF preadipocytes shared the same fibroblast-like morphology as the AbF preadipocytes (Fig. 1b). To construct intramuscular and abdominal adipogenic differentiation models, MDI medium supplemented with oleic acid was used for adipogenic differentiation. After induction with adipogenic agents for 10 days, chicken preadipocytes readily differentiated into mature adipocytes, and lipid droplets were visible under a microscope after 10 days of induction (Fig. 1b). AbF-derived preadipocytes exhibited a higher adipogenic differentiation ability (96.4% + 0.6) than IMFderived preadipocytes (86.0% + 0.4) (p < 0.01). The expression level of the adipogenic marker genes peroxisome proliferator-activated receptor γ (PPARγ, PPARG) and fatty acid binding protein 4 (FABP4, ap2) significantly increased with adipogenic differentiation (p < 0.01) (Fig. 1c). Pearson correlation analysis showed that the gene expression correlation coefficient within each group was noticeably higher than that between the groups, reflecting a good linear correlation between the independent samples of preadipocytes or adipocytes in the IMF and AbF groups (Pearson correlation coefficient, r = 0.98) (Fig. 1d). Principal component analysis (PCA) also showed global differences among the preadipocyte, adipocyte, IMF and AbF groups (Fig. 1e). All evidence suggested that our data was repeatability and reproducibility.

Global analysis of gene expression patterns in chicken adipocytes
As shown in Table 1, 963,374,122 raw reads were produced from 8 cDNA libraries. We identified a considerable number of genes in preadipocytes and adipocytes derived from chicken breast muscle and abdominal tissues. The percentage of clean reads in each library ranged from 94.21 to 96.09%. The mapping percentage of the 8 samples ranged from 88.82 to 93.42%. We found that the genomic loci of the genes were widely distributed across chromosomes (Fig. 2a). The sequencing depth was saturated at 68 M reads for each library (Fig. 2b). The mapping percentages of the different samples on different regions of the genome are displayed in Fig. 2c. More than 75% of the reads were mapped to gene regions. In preadipocytes and adipocytes derived from breast muscle and abdominal tissues, less than 0.1% of the reads were mapped to splicing sites of the genome.

DEG, DEM and DENG profiles between different chicken adipose tissue-derived preadipocytes and adipocytes
To identify potential candidate genes related to adipogenic differentiation, we examined the expression level of genes in preadipocytes and adipocytes. A total of 2039 genes were found to be differentially expressed between preadipocytes and adipocytes in both the IMF and AbF groups (Fig. 3a). As expected, we noticed that large numbers of genes or transcription factors related to adipogenic differentiation and lipid metabolism were differentially expressed, including peroxisome proliferator-activated receptor γ (PPARG), bone morphogenetic protein 4 (BMP4), fatty acid synthase (FASN), adiponectin, C1Q and collagen domain-containing (ADIPOQ), perilipin 2 (PLIN2), lipin 1 (LPIN1), and carnitine palmitoyltransferase 1A (CPT1A) (Fig. 3b). The ten most abundant DEGs between preadipocytes and adipocytes in the IMF and AbF groups are presented in Table 2. Furthermore, we determined the global gene expression profiles in preadipocytes and adipocytes between different groups. The number of differentially expressed genes (DEGs) (Additional file 1: Table S1), mRNAs (DEMs) (Additional file 2: Table S2) and novel genes (DENGs) (Additional file 3: Table S3) between the different groups are shown in Fig. 4a. As illustrated in Fig. 4b,  We found that 2742 DEGs were upregulated and 1705 DEGs were downregulated in the AbAd group compared with the IMAd group, that 3437 DEGs were upregulated and 1310 DEGs were downregulated in the AbPre group compared with the IMPre group, that 2797 DEGs were upregulated and 1896 DEGs were downregulated in the AbPre group compared with the AbAd group, and that 2055 DEGs were upregulated and 2348 DEGs were downregulated in the IMPre group compared with the IMAd group (Fig. 4c).

Functional enrichment analysis of DEGs involved in the adipogenic differentiation of chicken preadipocytes
To investigate the functions of the DEGs in chicken adipogenic differentiation, GO (Additional file 4: Table S4) and KEGG pathway (Additional file 5: Table S5) analyses were performed in the present study. Our results suggested that the DEGs in the AbPre vs AbAd and IMPre vs IMAd comparisons were significantly enriched in ECM-receptor interaction, the PPAR signalling pathway, and focal adhesion (Fig. 5). Interestingly, we found that ABC transporters, glutathione metabolism and fatty acid biosynthesis were significantly enriched for the IMF groups, while steroid biosynthesis and the p53 signalling pathway were significantly enriched for the AbF groups. The DEGs in the IMPre vs AbPre and IMAd vs AbAd comparisons were significantly enriched in the focal adhesion, ECM-receptor interaction, and fatty acid metabolism pathways, among others ( Fig. 5). To identify the gene expression patterns associated with adipogenic differentiation in both IMF-and AbF-derived adipocytes, the DEGs shared between the IMF group and the AbF group were analysed in the present study. Our results suggested that the shared DEGs were enriched for the metabolism, cellular processes and translation terms (Fig. 6a) and for pathways including ECM-receptor interaction, DNA replication, the cell cycle and the PPAR signalling pathway (Fig. 6b). Furthermore, we found that the shared downregulated DEGs (number: 750) were significantly enriched in the cell cycle, focal adhesion and purine metabolism pathways, while the upregulated DEGs (number: 792) were significantly enriched in the PPAR signalling pathway, the adipocytokine signalling pathway and the retinol metabolism pathway (Fig. 7).

Differentially expressed transcription factors and their potential interacting genes involved in adipogenic differentiation
To identify the differential gene expression patterns of transcription factors in the context of adipogenic differentiation, we compared the expression levels of transcription factors (TFs). Interestingly, our results suggested that most TFs shared the same gene expression patterns in both the IMF and AbF groups, while KLF9 and MYOG showed significantly different expression patterns between the groups (Fig. 8). Furthermore, as shown in Fig. 9, the TFs and genes related to adipogenic differentiation were significantly positively correlated (r > 0.86, p < 0.01).

Genes expressed in a fat depot-specific manner and the validation of the DEGs by qRT-PCR
To identify whether gene expression is fat depot-specific ( Fig. 10), the expression levels of genes related to adipogenic differentiation at different differentiation stages and in different groups (IMF and AbF) were compared. Our results suggested that most genes have different gene expression patterns during the adipogenic differentiation process. Most genes reached the highest expression levels after 6 days of induced differentiation in the AbF group and after 8 days in the IMF group (Fig. 11). To confirm the results from RNA-Seq, qRT-PCR was performed in the present study. Our qRT-PCR results were consistent with the RNA-Seq results (Additional file 7: Figure S1).

Integrated analysis of DEG-pathway network between different fat-derived chicken preadipocytes and adipocytes
To further understand the adipogenic differentiationregulated mechanisms of different fat-derived chicken preadipocytes, we visualized the integrated DEGpathway networks for the IMF (Fig. 12a) and AbF (Fig. 12b) groups. Our results showed that a large number of DEGs in the IMF group were mainly enriched in the PPAR signaling pathway, fatty acid biosynthesis, focal adhesion and ECM-receptor interaction (Fig. 12a).
For the AbF group, we noticed that the DEGs were mainly enriched in the PPAR signaling pathway, focal adhesion, ECM-receptor interaction, steroid biosynthesis and the p53 signaling pathway (Fig. 12b).

Discussion
Fat deposition is mainly dependent on the proliferation, differentiation and maturation of preadipocytes [22]. IMF content is an important factor that contributes to the tenderness, juiciness and flavour of meat and thus affects meat quality. High AbF content causes low slaughter efficiency [5,7]. It is known that IMF has high genetic correlations with abdominal fat weight (AFW) and moderate correlations with AF percentage (AFP) in chickens [23]. Previous studies have suggested that the proliferation and differentiation abilities of genes involved in lipid metabolism are dramatically lower in intramuscular adipocytes (IMAs) than in subcutaneous adipocytes (SAs) [11,14,[24][25][26]. In the present study, lipids accumulated in chicken IMF and AbF preadipocytes at the late stage of differentiation. The oil red O staining results and the expression levels of two well-known adipogenic markers demonstrated that the model of adipocyte differentiation was successfully established. We noticed that the accumulation of lipids in chicken AbF adipocytes was higher than that in IMF adipocytes, consistent with the findings of previous studies in pigs [11,27,28]. Ribo-Zero RNA-Seq has been applied as an efficient method to explore transcriptional characteristics because it can capture both poly(A) + and poly(−) transcripts [9,29,30]. However, to our knowledge, none of the previous studies were carried out to compare gene expression profiles between IMF and AbF adipocytes, especially in chickens. The process by which preadipocytes differentiate into mature adipocytes is complex and is regulated by various transcription factors [31]. Previous studies have identified a large number of transcription factors, including CCAAT/enhancer-binding protein (C/EBP), peroxisome proliferator-activated receptors (PPARs) and sterol regulatory element-binding protein (SREBP) [32][33][34][35].
In the present study, functional annotation analysis of the DEGs revealed that these genes play important roles in some lipid metabolism-and adipogenic differentiation-related pathways, such as the PPAR signalling pathway, ECM-receptor interaction and fatty acid metabolism [34,43,46]. PPARG, the most adipocyte-specific and adipogenic member of the PPAR family, is mainly expressed in adipose tissue and plays an important role in lipid metabolism and adipocyte differentiation [47,48]. The extracellular matrix (ECM) plays an important role in the regulation of proliferation, adipogenic differentiation and migration of preadipocytes [43]. Adipocyte differentiation can also be affected by fatty acid metabolism via the regulation of transcription factors [34,49].
We noticed that APOA1, CD36, LAMA2, CHAD, MMP1, RRM2, VWF and PIK3R5 were the most upregulated genes during the adipogenic differentiation processes between the   IMF and AbF groups. Therefore, these genes may be involved in the positively regulating position specificity fat deposition. In addition, the differences in secretory functions and hormone sensitivities between IMF and AbF might be caused by the position-specific regulation of adipose tissue in poultry. Further studies are necessary for to elucidate the

Conclusions
In conclusion, our current study showed that abdominal fat (AbF) preadipocytes accumulate more lipids than intramuscular fat (IMF) preadipocytes. This study presents the first analysis of gene expression during the differentiation of intramuscular and abdominal preadipocytes in chickens. A total of 2039 DEGs were identified by a pairwise comparison of preadipocytes at different stages of differentiation. The DEGs were found to be involved in the PPAR signalling pathway, fatty acid biosynthesis, ECMreceptor interaction and focal adhesion, consistent with previous reports on preadipocyte differentiation in chickens. These DEGs and pathways might play significant roles in intramuscular preadipocyte differentiation in chickens. Our findings provide a solid foundation for future studies on the molecular mechanisms underlying tissue-specific fat deposition and on strategies for the improvement of meat quality in poultry.

Ethics statement
All

Primary preadipocyte isolation and culture in vitro
Primary preadipocytes were isolated from the breast muscles and abdominal fat of two-week-old Gushi chickens according our previously described method [50]. In brief, breast muscle and abdominal adipose tissues were separated from the body of each chicken under sterile conditions. The tissues were washed using phosphate-buffered saline (PBS) supplemented with penicillin (100 units/mL) and streptomycin (100 μg/mL). The washed tissue was cut into 1-mm 3 pieces and then digested with collagenase type II (1 mg/mL, Solarbio, Beijing, China) at 37°C for 90 min. The digested cell suspension was filtered using 200-and 500-mesh screens to separate the stromal-vascular fraction from undigested tissue and mature adipocytes, and the fraction was then centrifuged at 1000 x g for 5 min.

Read mapping and transcriptome assembly
The splice-mapping algorithm of HISAT2 (2.0.4) [51] was used to perform genome mapping of the preprocessed reads. The clean data were mapped to the Gallus gallus reference genome (GGA5) with HISAT2 (2.0.4), and the default parameters were used.

Oil red O staining
An Oil Red O staining assay was performed according to the methods in our previous study [50]. Briefly, adipocytes were gently washed with cold PBS three times and fixed with 4% paraformaldehyde for 30 min. Then, the fixed cells were gently washed with cold PBS three times, incubated with 60% filtered Oil Red O solution for 40 min and then observed under a phase-contrast microscope to check for positive cells appearing red. The cells were then washed three times with deionized cold PBS and photographed using an Olympus CKX41-F32FL microscope (Olympus, Tokyo, Japan). Subsequently, Oil Red O was eluted from the stained cells with 100% isopropanol (v/v) and quantified with a microplate reader (Thermo Fisher Scientific) at 500 nm. ImageJ software (National Institutes of Health, Bethesda, MD) was used to estimate the adipogenic differentiation ability. The cells with lipid droplets were regarded as differentiated cells: Adipogenic differentiation ratio = (differentiated cells count / total cells count) × 100%.

RNA isolation and real-time quantitative PCR (qRT-PCR)
Primers for the DEGs were designed using Primer3Plus online software (http://www.primer3plus.com/cgi-bin/ dev/primer3plus.cgi) (Additional file 6: Table S6). qPCR was performed using SYBR® Green PCR Master Mix (TaKaRa, Dalian, China). The PCR mixture contained 5 μL of SYBR® Premix Ex Taq II (2×), 0.5 μL of forward primer (10 μM), 0.5 μL of reverse primer (10 μM), and 200 ng of cDNA, with RNA-free water added to 10 μL. The qRT-PCR was conducted in a LightCycler 96 system (Roche) by the SYBR Green method. The program included an initial step of 95°C for 5 min; 40 cycles of 95°C for 30 s, 60°C for 30 s, and 72°C for 20 s; and melting curve generation, which was performed as follows: 95°C for 10 s, annealing at 65°C for 20 s, and heating through a continuous temperature gradient to 97°C with 5 acquisitions/s. Chicken GAPDH was selected as an internal control gene. All samples were examined in triplicate. All data were analysed using the 2 -ΔΔCt method. All data are shown as fold changes in gene expression compared with the gene expression in the 0d group.

Statistical analysis
Statistical significance between two experimental groups was evaluated with a T-test for comparisons in SPSS 20.0 statistical software (IBM, Chicago, IL, USA). Statistical significance among three or more experimental groups was evaluated by one-way ANOVA followed by Dunnett's test for multiple comparisons in SPSS 22.0. Graphics were drawn using GraphPad Prism 7 (Graph-Pad Software, San Diego, CA, USA) and RStudio 1.1.453 software. All data are expressed as the mean ± standard error (SE). A P-value ≤0.05 was considered statistically significant, and a P-value ≤0.01 was considered extreme significant.