Skip to main content

Advertisement

Dynamic transcriptome and DNA methylome analyses on longissimus dorsi to identify genes underlying intramuscular fat content in pigs

BMC Genomicsvolume 18, Article number: 780 (2017) | Download Citation

Article metrics

  • 1356 Accesses

  • 9 Citations

Abstract

Background

The intramuscular fat content (IMF) refers to the amount of fat within muscles, including the sum of phospholipids mainly found in cell membranes, triglycerides and cholesterol, and is determined both by hyperplasia and hypertrophy of adipocyte during the development of pigs. The IMF content is an important economic trait that is genetically controlled by multiple genes. The Laiwu pig is an indigenous fatty pig breed distributed in North China, characterized by excessively higher level of IMF content (9%~12%), therefore, is suitable for the identification of genes controlling IMF variations. To identify genes underlying IMF deposition, we performed genome-wide transcriptome and methylome analyses on longissimus dorsi (LD) muscle in Laiwu pigs across four developmental stages.

Results

A total of 22,524 expressed genes were detected and 1158 differentially expressed genes (DEGs) were hierarchically clustered in the LD muscle over four developmental stages from 60 d to 400 d. These genes were significantly clustered into four temporal expression profiles, and genes participating in fat cell differentiation and lipid biosynthesis processes were identified. From 120 d to 240 d, the period with the maximum IMF deposition rate, the lipid biosynthesis related genes (FOSL1, FAM213B and G0S2), transcription factors (TFs) (EGR1, KLF5, SREBF2, TP53 and TWIST1) and enriched pathways (steroid biosynthesis and fatty acid biosynthesis) were revealed; and fat biosynthesis relevant genes showing differences in DNA methylation in gene body or intergenic region were detected, such as FASN, PVALB, ID2, SH3PXD2B and EGR1.

Conclusions

This study provides a comprehensive landscape of transcriptome of the LD muscle in Laiwu pigs ranging from 60 to 400 days old, and methylome of the LD muscle in 120 d and 240 d Laiwu pigs. A set of candidate genes and TFs involved in fat biosynthesis process were identified, which were probably responsible for IMF deposition. The results from this study would provide a reference for the identification of genes controlling IMF variation, and for exploring molecular mechanisms underlying IMF deposition in pigs.

Background

Pigs are the main sources for human meat consumption and provide 43% of total meat production in the world [1, 2]. Skeletal muscle is the primary meat production tissue of pigs, responsible for 20%~50% of their body mass [3]. The intramuscular fat content (IMF) refers to the amount of fat within muscles, which includes the sum of phospholipids mainly found in cell membranes, triglycerides serving as the main forms of energy reserves and cholesterol [4]. The triglycerides in mammalian and avian muscles are mainly stored within intramuscular adipocytes and myofibres cytoplasm in droplets in close vicinity to mitochondria [4,5,6], which appear as the flecks and streaks of fat within the lean sections of meat [4, 7]. The IMF is considered as late developing fat storage, which is determined both by hyperplasia and hypertrophy of adipocyte during the development of pigs [8]. The fatty acid composition of these fat components in muscle tissue is associated with eating quality [9]; saturated and monounsaturated fatty acids (SFA and MUFA) are positively, while polyunsaturated fatty acids (PUFA) are negatively correlated with meat flavour [10].

The IMF content is a polygenic trait in livestock species that affects meat quality traits such as flavour, drip loss and shear force [4, 11]. Studies on the genetic basis of IMF in pigs have been extensively carried out, and some genes associated with IMF were reported [12, 13], however, the molecular mechanisms of porcine IMF are still unclear. The Laiwu pig is an indigenous fatty pig breed distributed in North China, characterized by excessively higher level of IMF content (9%~12%), far exceeding the commercial lean pigs [14], therefore, is suitable for the identification of genes controlling IMF variations. Several studies have analyzed the genetic basis of IMF deposition in Laiwu pigs using candidate gene approach [15,16,17]. Recently, genome-wide analysis of differentially expressed genes (DEGs) has been proved as an effective approach to elucidating genetic mechanisms of complex traits. The transcriptome of porcine subcutaneous adipose [18, 19], muscle [2, 20] and liver [21, 22] have been investigated, while genome-wide analysis on the molecular mechanisms of IMF deposition in Laiwu pigs is still lacking.

In this study, we present a comprehensive survey of the transcriptome profiles of the longissimus dorsi muscle (LD) in Laiwu pigs across 60, 120, 240 and 400 days of age, which are representative of major morphological and developmental stages in Laiwu pigs. Moreover, DNA methylation at CpG dinucleotide in the promoter region plays an important role in regulating the transcription of genes [21, 23]. Based on the transcriptome results, we further compared the changes in DNA methylation profile during the fastest fat deposition period (from 120 d to 240 d) and analyzed the role of DNA methylation in the expression dynamics of genes related to IMF changes. The results of this study revealed several genes and pathways involved in adipogenesis and lipogenesis in LD muscle that likely underlie IMF trait in pigs.

Results

Dynamics in IMF and fatty acid composition of porcine LD muscle across four developmental stages

From 60 d to 400 d, along with the growth of pigs, the IMF content of porcine LD muscle increased significantly (p < 0.01), and from 120 d to 240 d, it is increased from 3.59% to 9.88%, representing the fastest fat deposition stage of the LD muscle in Laiwu pigs (Fig. 1a, b). The saturated fatty acid proportion of both palmitic acid (C16: 0) and stearic acid (C18: 0) significantly increased from 120 d to 240 d, then decreased at 400 d (Fig. 1c, d). Oleic acid (C18: 1) is the most abundant unsaturated fatty acid, accounting for approximately 40% of total muscle fatty acids and increased slightly over time (Fig. 1e). By contrast, the proportion of linoleic acid (C18: 2n-6) decreased significantly from 60 d to 240 d (Fig. 1f). The SFA and MUFA accounted for relatively larger proportion of fatty acid in LD muscle from 240 d to 400 d, whereas the opposite is for the PUFA (Fig. 1g).

Fig. 1
figure1

The phenotypic characteristics in the longissimus dorsi (LD) muscle of Laiwu pigs. Dynamics in the (a) live weight, (b) intramuscular fat (IMF) content and fatty acid composition of (c) palmitic acid (C16: 0), (d) stearic acid (C18: 0), (e) oleic acid (C18: 1n-9) and (f) linoleic acid (C18: 2n-6) across 60, 120, 240 and 400 days of age. (g) The proportion of saturated fatty acids (SFA), monounsaturated fatty acids (MUFA) and polyunsaturated fatty acids (PUFA). The data show the means ± SD analyzed by one-way ANOVA followed by Duncan’s multiple comparison. a-c: p ≤ 0.05, A-D: p ≤ 0.01

Full size image

Summary of RNA-seq sequencing

High-throughput RNA-seq generated 54.19 Gb clean data for 12 LD muscle samples of Laiwu pigs across four postnatal developmental stages of 60, 120, 240 and 400 d, and 70.15%–73.61% reads were successfully aligned to pig genome. All 12 samples had at least 88.47% reads equal to or exceeding Q30 (Table 1). The pair wise correlation was more than 0.83 among individuals of each group (Fig. 2a). A total of 22,524 expressed genes were detected and 1158 DEGs were hierarchically clustered over four developmental stages (Fold-Change ≥ 2 and FDR ≤ 0.01) (Fig. 2b). Six genes including C/EBPA, C/EBPD, MSTN, PDK4, PPARG and SCD were chosen and quantified using qRT-PCR to validate sequencing data, and a strong correlation (r > 0.76) between qRT-PCR and RNA-seq data was observed, suggesting that the RNA-seq result was reliable (Additional file 1: Figure S1).

Table 1 Summary of RNA-Seq metrics from transcriptomes across four developmental stages of Laiwu pigs
Full size table
Fig. 2
figure2

Hierarchical clustering analysis for all the samples and differentially expressed genes (DEGs). a Heat map matrix of 12 LD muscle samples of Laiwu pigs constructed using Pearson’s correlation. b Hierarchical clustering analysis for DEGs in the longissimus dorsi (LD) muscle of Laiwu pigs across four developmental stages from 60 d to 400 d

Full size image

Classification of the expressed genes in porcine LD muscle

We clustered the expressed genes (at least one group FPKM ≥ 0.1) in the LD muscle of Laiwu pigs using Short Time-series Expression Miner (STEM) to obtain statistically significant clustered temporal expression profiles [24]. A total of 10 expression profiles were obtained, and the first four profiles were significantly clustered, including 1423 genes (p < 0.01, Fig. 3a). Among these profiles, profile 1, 3 and 4 showed a robust correlation and an expression pattern similar to IMF variations in the LD muscle of Laiwu pigs from 60 d to 400 d (Fig. 3b). The genes from these three files (Additional file 2: Table S1) were subsequently uploaded into ClueGO, and those genes participating in fat cell differentiation, lipid metabolism processes and skeletal muscle development were identified, including ACACA, SCD, ACLY, ELOVL1 and FASN that participate in long-chain fatty-acyl-CoA, triglyceride and lipid biosynthetic processes; APOE, APOA1, TSPO, FGFR4, FDPS, G6PD, HSD17B7, LSS, ADM, SQLE, DHCR24, CYP51A1 and SCARB1 that are involved in steroid biosynthetic processes (Fig. 3c); NFATC2, TMEM8C, MYOG, MYOD1, MYF5 and MYF6 that are involved in myoblast fusion and ID2, TNNT2, MEGF10 and MYL3 that participate in muscle tissue development (Additional file 3: Figure S2).

Fig. 3
figure3

Analysis of short time-series gene expression cluster across four developmental stages. a Four significant cluster profiles under the order of time points. ‘n’ represents the gene numbers of each profile. b The relationship between these four significant clustered profiles was calculated based on correlation coefficients. Strongly positive or negative correlations between gene pairs are shown in black or red, respectively. c The differentially expressed genes and functions associated with fat metabolism in profile 1, 3 and 4. GO terms and genes are represented as nodes based on their kappa score more than 0.4 and networks with at least three nodes. The node size represents the GO terms enrichment significance

Full size image

Expression analysis on IMF relevant QTL in the LD muscle of Laiwu pigs

To validate the DEGs potentially related to pig IMF content trait, we examined the expression changes of IMF content relevant QTL across four development stages of Laiwu pigs based on pig QTL database (http://www.animalgenome.org/cgi-bin/ QTLdb/SS/index). A total of 14 specifically expressed genes in LD muscle were detected, most of them have the highest FPKM value at 240 d, including SCD, FASN, FADS2, ANK1, ACSF3 and SNAI2, and three genes including PIK3R6, PIK3C3 and CA3 showed the peak FPKM value at 400 d (Table 2).

Table 2 The gene expression changes of IMF relevant QTLs across four developmental stages of Laiwu pigs
Full size table

DEGs and the pathway enrichment between adjacent developmental stages in porcine LD muscle

The DEGs between adjacent developmental stages, i.e., 60 d vs 120 d, 120 d vs 240 d and 240 d vs 400 d, were screened (Additional file 4: Table S2), and pathways with these DEGs were performed using DAVID. From 60 d to 120 d, 208 genes were differentially expressed, including 155 up-regulated and 53 down-regulated genes, which are significantly enriched in 11 pathways including cell proliferation (ECM-receptor interaction, focal adhesion) and immune (B cell receptor signalling pathway) (p < 0.01) (Fig. 4a). From 120 d to 240 d, 881 genes were differentially expressed, including 628 up-regulated and 253 down-regulated genes, which are involved in 16 pathways including phagocytosis (phagosome, endocytosis), bone differentiation (osteoclast differentiation), energy metabolism and fat deposition (type I diabetes mellitus, insulin signalling pathway, steroid biosynthesis, fatty acid biosynthesis) and protein metabolism (protein digestion and absorption, arginine and proline metabolism) (p < 0.01) (Fig. 4b). From 240 d to 400 d, 251 genes were differentially expressed, including 83 up-regulated and 168 down-regulated genes, which are involved in 13 pathways including immune (cytokine-cytokine receptor interaction, primary immunodeficiency, chemokine signalling pathway), energy metabolism and fat deposition (insulin signalling pathway, fatty acid biosynthesis) (p < 0.01) (Fig. 4c). Pathways related to fat deposition, including fatty acid biosynthesis, steroid biosynthesis, biosynthesis of unsaturated fatty acid and adipocytokine signalling pathways, were revealed in the comparison between 120 d and 240 d. The fatty acid biosynthesis pathway was also significantly enriched in 400 d, but most of the DEGs were down-regulated.

Fig. 4
figure4

The differentially expressed genes and enriched pathways in the LD muscle of Laiwu pigs. a 60 d vs 120 d, (b) 120 d vs 240 d and (c) 240 d vs 400 d. The number in circle represent the detected expressed gene number, the intersection was the differentially expressed gene number, red and green represented the up-regulated and down-regulated genes, respectively. Grey dashed line in pathway indicates a threshold of p = 0.05

Full size image

Functional enrichment in porcine LD muscle between 120 d and 240 d

From 120 d to 240 d, the fastest IMF deposition stage, the top 10 significantly up- or down-regulated DEGs were found to be involved in biological process such as energy metabolism (ATP5J2, DNAJB1), transcription (EGR1, KLF11), muscle development (CHRNG, TNNT2) and fat metabolism (G0S2, CYP1A1, FAM213B) (Table 3). We uploaded 881 DEGs to IPA for functional analysis and revealed 127 DGEs (93 up-regulated and 34 down-regulated) that are involved in lipid synthesis process (Z-score = 1.25, p < 0.001). Among those genes related to lipid synthesis, the top 5 up-regulated genes were CYP1A1, SERPINA1, LDLR, EGR1 and FOSL1, and the top 5 down-regulated genes were DIAPH1, SORBS1, PDK4, ACSL1 and ASPA (Fig. 5), and 84 genes were included in profiles 1, 3 and 4 (Additional file 5: Table S3), such as FASN, ACACA, ACLY and SCD.

Table 3 The top 10 up- and down-regulated genes of 120 d vs 240 d Laiwu pigs
Full size table
Fig. 5
figure5

Candidate genes associated with fat biosynthesis process between 120 d and 240 d. The DEGs are organized by cellular substructure, up-regulated genes in 240 d are indicated in red, down-regulated genes in 240 d are indicated in green and the colour depth represents the regulated degree. IPA predicted activation Z-score = 1.25, p-value = 3.73E-10

Full size image

Furthermore, the upstream key regulators of IMF deposition were predicted using IPA [2, 8], considering gene expression is transcriptionally controlled by upstream regulators, particularly transcription factors (TFs). Even if TFs do not always exhibit some extent change, the expression of downstream genes would vary significantly [8]. A total of 1620 upstream regulators were predicted and 16 regulators would be expected to have a greater impact on IMF deposition based on their biological function annotations (Table 4). Meanwhile, the predicted downstream genes (Additional file 6: Table S4) and the binding site (Additional file 7: Table S5) of TFs were summarized. The TFs including EGR1, KLF5, SREBF2, TP53 and TWIST1 were activated (Z-score ≥ 2.0) [25], and the expression of EGR1 and TP53 was significantly up-regulated in 240 d compared with 120 d.

Table 4 Upstream regulators in the LD muscle of Laiwu pigs from 120 d to 240 d
Full size table

mRNA comparison of Laiwu pigs with other pig breeds at 120 d

At 120 d, IMF content in the LD muscle of Laiwu pigs was significantly higher than that of DYL pigs (Additional file 8: Figure S3). The qRT-PCR analysis showed that the mRNA expression of FADS2, FABP4, ACSF3 and ATP5B in LD muscle was significantly higher in Laiwu pigs than in DYL pigs, while that of LCLAT1, SCD, LDLR and CA3 was not significantly different (Fig. 6). Meanwhile, Ovilo [26] has reported that the IMF of IB × DU pig was 2.87% at 120 d, which was lower than that of Laiwu pig (3.59%) at the same age. The transcriptome data of IB × DU pig at 120 d were downloaded from GEO of NCBI (GSE86441). A total of 682 DEGs were screened (Additional file 9: Table S6), and 49 IMF candidate genes were obtained based on their GO or KEGG annotations, of which 37 genes were higher and 11 genes were lower in the LD muscle of Laiwu pigs compared with IB × DU pigs (Additional file 10: Table S7).

Fig. 6
figure6

The detection of IMF candidate genes between DYL and Laiwu pigs using qRT-PCR. *p < 0.05; **p < 0.01

Full size image

Landscape of DNA methylomes of porcine LD muscle

To further understand the role of DNA methylation in regulating the transcription of genes related to IMF deposition from 120 d to 240 d, we performed reduced representation bisulfite sequencing (RRBS) on the genomic DNA from the LD muscle of Laiwu pigs. Approximately 12.2G clean data were obtained in each group with 46.27% to 46.64% clean reads uniquely mapped (Table 5). We plotted the genome-wide distribution of cytosine according to sequencing read coverage and depth across chromosomes (Additional file 11: Figure S4) and the methylation distribution of genes within different functional components, including 2 kb region upstream, gene body and downstream region (Fig. 7a). Similar CpG methylation trend was shown between 120 d and 240 d: a relative low methylation levels were found in upstream region, reaching the lowest level around transcription start sites (TSS), being stable in the first exon, followed by a sharp increase in the first intron, and hitting a plateau till transcription termination site (TTS).

Table 5 Data generated by RRBS of LD muscle in 120 d and 240 d Laiwu pigs
Full size table
Fig. 7
figure7

The methylation distribution around gene body of Laiwu pigs between 120 d and 240 d. a Methylation level of different sites around gene body. b Distribution of differentially methylated regions (DMRs) in gene body and intergenic region. Most (52.5%) DMRs located in intergenic region, followed by gene body (34.7%) and promoter (24.3%) regions. TSS: transcription start sites; TTS: transcription termination site

Full size image

A total of 124,693 differentially methylated cytosines (DMCs) between 120 d and 240 d were identified, with sequencing depth ≥ 4, and FDR ≤ 0.05 (Additional file 12: Table S8). Similarly, a total of 3676 differentially methylation regions (DMRs) were generated corresponding to 2303 genes. The DMRs were screened by the sequencing depth ≥ 10×, at least 3 DMCs, with the minimum value ≥0.3, and p < 0.05 using Fisher’s exact test. Then, these DMRs were annotated according to the up- or down-stream 3000 bp of gene using MOABS method [27] (Additional file 13: Table S9). Distribution of the DMRs showed that intergenic region accounted for the most (52.5%) methylated variations, followed by gene body (34.7%), and gene promoter (24.3%) (Fig. 7b). The number of DMRs distributed on genes was shown in Table 6, hyper-methylated regions generally exceeding hypo-methylated regions in 240 d than 120 d. In addition, we performed a GO functional enrichment analysis for genes with DMRs in their promoters [28]. The top 12 significantly enriched GO terms were mainly involved in transcription, myotube differentiation and osteoblast differentiation (Table 7). Notably, the fat metabolism process was also significantly enriched, the involved genes include FASN, LEP, LCN12, GDPD3 and FADS6.

Table 6 Number of differentially methylated genes in DMRs in the LD muscle of Laiwu pigs
Full size table
Table 7 Top 15 Gene Ontology (GO) categories enriched for genes with DMRs in their promoters
Full size table

Integration analysis of transcriptome and DNA methylome in porcine LD muscle between 120 d and 240 d

The gene expression level was divided into four groups according to FPKM (highest, lowest, medium high, medium low) to analyze the relationship between methylation variation and mRNA expression levels of genes in LD muscle (Additional file 14: Table S10). The expression level of the genes among the four groups was significantly different. Analysis showed that the lowest methylation level was found around the TSS in both 120 d and 240 d groups, the highest FPKM group was significantly different from others (Additional file 15: Figure S5), with an obviously sharp decrease in their upstream regions (Fig. 8). Furthermore, we calculated the proportions of DEG with DMR (4.3%), DEG without DMR (22.4%) and DMR without DEG (73.3%) (Additional file 16: Figure S6), and detected a total of 127 DEGs with DMRs, suggesting a likely role of methylation on transcription (Additional file 17: Table S11). The lipid biosynthesis relevant genes showed different DNA methylation level in gene body or intergenic region on 120 d and 240 d, such as FASN, PVALB, ID2, SH3PXD2B and EGR1 (Table 8).

Fig. 8
figure8

CpG methylation of genes with different expression level between 120 d and 240 d. Gene expression level was divided into four equal parts according to FPKM value. The red, blue, green, purple line indicate the highest, lowest, medium high and medium low FPKM value, respectively. a 120 days old Laiwu pigs; and (b) 240 days old Laiwu pigs. TSS: transcription start sites; TTS: transcription termination site

Full size image
Table 8 The fat relevant genes with DMRs located in between 120 d and 240 d
Full size table

RRBS data validation of EGR1 and FASN via BSP

According to transcriptome and RRBS data, the two candidate IMF genes EGR1 and FASN were differentially expressed and with DMRs in 120 d vs 240 d. Herein, they were validated through qRT-PCR and BSP. The mRNA expression level of EGR1and FASN were up-regulated in 240 d (Fig. 9a, e), in consistent with transcriptome results. BSP analysis revealed that the proportion of CpG methylation in the intergenic region of EGR1 was significantly lower in the LD muscle of 240 d pigs than that of 120 d (p < 0.0001), and that of FASN promoter region was not significantly different (p = 0.0575, Fig. 9d, h). The BSP analysis results on EGR1 and FASN were basically consistent with RRBS data.

Fig. 9
figure9

The detection of mRNA expression and methylation of candidate genes. a The relative mRNA expression of EGR1 between 120 d and 240 d Laiwu pigs; The (b)120 d and (c) 240 d EGR1 intergenic CpG methylation patterns; (d) The methylated cytosine proportion of EGR1 intergenic CpG region between 120 d and 240 d. e The relative mRNA expression of FASN between 120 d and 240 d Laiwu pigs; The (f)120 d and (g) 240 d FASN promoter CpG methylation patterns; (h) The methylated cytosine proportion of FASN promoter CpG region between 120 d and 240 d. *p < 0.05; ***p < 0.001

Full size image

Discussion

In this study, by a combinational use of transcriptome and DNA methylome analysis, we screened the genes and pathways related to IMF deposition in the LD muscle from Laiwu pigs across four developmental satges. The Laiwu pigs exhibit strong capacity to deposit fat in skeletal muscles, therefore, are ideal animal model for analyzing mechanisms underlying IMF deposition. To identify genes that are related to IMF deposition in Laiwu pigs, we first obtained transcriptome of LD muscle from Laiwu pigs of 60, 120, 240 and 400 days old, and selected genes with similar dynamics to IMF; then compared the DEGs between adjacent stages, especially focused on the genes and pathways between 120 d and 240 d, the stage of quickly increased IMF deposition; and finally analyzed the DNA methylation changes in the IMF related genes in the LD muscle from 120 d to 240 d.

By STEM analysis on the transcriptome data of four developmental stages, we found that three temporal expression profiles correlated with IMF variation in the LD muscle. From these profiles, 15 DEGs with known function in lipid metabolism were identified. Among these genes, ACACA, SCD, ACLY, ELOVL1 and FASN are lipogenesis genes, all of them involved in long-chain fatty-acyl-CoA biosynthetic process and played regulatory and/or catalysis roles in fatty acid biosynthesis [4, 29,30,31]. Four genes (SCD, FADS1, FADS2, and FADS3) are known to encode desaturases, responding for the biosynthesis of unsaturated fatty acids [32,33,34,35,36]. Six genes (CYP51, DHCR24, EBP, HSD17B7, SOAT1 and SQLE) are involved in steroid biosynthesis pathway, execute their catalytic function during cholesterolopoiesis process [37, 38]. The expression dynamics of these genes from 60 d to 400 d are in accordance with the trend of IMF deposition dynamics in the LD muscle of Laiwu pigs. Moreover, since the growth and development of LD muscle is gradually increased from 60 d to 400 d of Laiwu pigs, reflected in the changes in skeletal muscle fibre area, genes involved in myoblast fusion and muscle tissue development were also identified (Additional file 3: Figure S2), which will be further investigated in future study.

IMF deposition in skeletal muscle results from balance among the uptake, synthesis and oxydrolysis of lipids, and the development of IMF is a complex multi-organ process regulated by coordinated actions of muscle, adipocyte, and connective tissues; angiogenesis and immune system-related genes and pathways are also involved [4, 8, 39]. For instance, animals with high muscularity generally display a reduced developmental of IMF [4], suggesting that IMF candidate genes may negatively regulate muscle development. Blood vessels growth and remodelling not only supply optimal levels nutrients and oxygen to nourish adipocytes, but also provide adipose precursor and stem cells that control adipose tissue mass and function [18, 40, 41]. Accordingly, we observed the enrichment of multi-pathways in the LD muscle of Laiwu pigs from 60 d to 400 d, including cell proliferation and immune (from 60 d to 120 d), phagocytosis and protein metabolism (from 120 d to 240 d) (Fig. 4).

The 120 d vs 240 d group has the most DEGs and mostly up-regulated in 240 d. The fat deposition pathways firstly emerged in this period, concordant with the quick increase of IMF content, representing a critical fat deposition period. The fatty acid biosynthesis relevant DEGs were identified in this group, such as FASN, ACACA, ACLY and SCD. FASN (fatty acid synthase) catalyzes all of the reaction steps in the synthesis of palmitate from acetyl-CoA and malonyl-CoA, which plays a central role in de novo lipogenesis in mammals [29], and was reported to be related with pork IMF content [42, 43]. ACACA (Acetyl-CoA carboxylase alpha) is a key regulator of de novo fatty acid synthesis [30, 44] and the haplotype A 4899 C 5196 can decrease IMF percentage in a Duroc commercial population [45]. ACLY (ATP citrate lyase) is a rate-limiting lipogenic enzyme in pigs associated with fatty acid biosynthesis [31], its effects on IMF deposition has not been studied. Notably, SCD (stearoyl-coenzyme A desaturase) is also a rate-limiting enzyme in lipogenesis, regulating the biosynthesis of C16: 1 and C18: 1 from C16: 0 and C18: 0, respectively, and this process is regulated by SREBF1 [33, 34]. Bessa et al. [35] reported that IMF is positively related to SCD protein expression level in commercial pig breeds. Therefore, it is reasonable to speculate that it is the relative highly expressed SCD that triggers the increase of MUFA percentage [46], particularly between 120 and 240 days old Laiwu pigs.

Besides these well-known IMF relevant genes, we also focused on dramatically up- or down-regulated genes that probably relevant to fat metabolisms, such as LDLR, FOSL1, EGR1, G0S2, FAM213B and SORBS1. LDLR (low density lipoprotein receptor) plays an important role in lipid transportation and shows high positive correlation with IMF content in Piau pigs [47]. FOSL1 (FOS like 1, AP-1 transcription factor subunit) could be activated by the upstream AP1 motif, induced by Liver X receptors (LXRs) in mouse, which involves in regulating lipid synthesis and transport [48]. Boyle et al. [49] showed that EGR1 (early growth response transcription factor 1) is a negative regulator during 3 T3-L1 differentiation [50]. G0S2 (G0/G1 switch gene 2) can attenuate the activity of the intracellular triacylglycerol hydrolase and adipose triglyceride lipase [51]. And the role of the gene FAM213B (family with sequence similarity 213 member B) in fat deposition was rarely reported, which deserves further study. While the down-regulated gene SORBS1 (sorbin and SH3 domain containing 1) plays a key role in adipogenesis [52], and the expression level in the visceral adipose depots correlated negatively with body mass index [53]. Taken together, only LDLR has been reported related with IMF content among the six mentioned genes, the function of the other five genes in IMF deposition still deserves further investigations.

The well studied TFs relevant to lipid metabolism, such as KLF5, SREBF2, TP53 and TWIST1, were also identified. KLFs (Krüppel-like factors) regulate adipocyte differentiation in mammals, and KLF5 acts as a key regulator, by binding to SREBP1, to enhance the SREBP1-mediated increase in FASN promoter activity [54]. SREBFs (sterol regulatory element binding factors) are major regulators of carbohydrate and lipid metabolism [30]. SREBF2 is closely associated with cholesterol biosynthesis [55], showing an activation Z-score of 3.51 in comparison between 120 d and 240 d. TP53 (tumor protein p53) is involved in pig myogenesis [56] and functional p53 can regulate cellular glucose metabolism through several proteins such as TP53-induced glycolysis and phosphoglycerate mutase [57]. TWIST1 (twist family BHLH transcription factor 1) was demonstrated as a regulator during the adipocytes gene expression in 3 T3-L1 [58, 59]. The roles of these TFs in IMF deposition of pigs have not been reported.

It has been well studied that DNA methylation plays an important role in the regulation of gene expression [21, 23]. In this study, we obtained the DNA methylation landscapes of LD muscle in 120 d and 240 d Laiwu pigs. The results showed that most of the DMRs located in intergenic regions, which is consistent with another study [18]. The methylation level of CpG dinucleotide around TSS region showed a ‘V’ shaped curve in the LD muscle of both 120 d and 240 d Laiwu pigs, this pattern was also reported in other species, such as human [60], cattle [23] and rat [61]. Moreover, the DMRs located in promoter region were mainly involved in muscle, adipose and bone development, implying that methylation probably plays an important role during the LD muscle development.

Generally, the genes with high expressed levels are often associated with a relative lower promoter methylation [23, 62], while for the gene body regions, it remains unclear [63,64,65]. Studies assumed that DNA methylation in gene body is positively correlated with gene expression level, because DNA methylation might alter chromatin structure and transcription elongation efficiency [66, 67]. In this study, the expressed genes showed the lowest methylation level around TSS, and for the highly mRNA expressed genes both in the LD muscles of 120 d and 240 d, the methylation level in the promoter region was decreased, which is consistent with previous reports [63, 68, 69]. Integration analysis of transcriptome and methylome revealed a set of candidate IMF deposition genes probably regulated by methylation, such as FASN, ID2, ITGB2, PEMT and SH3PXD2B. These mentioned genes were all hypomethylated in 240 d in the promoter or intergenic regions. ID2 is an adipogenic transcription factor that stimulates PPAR expression and adipocyte differentiation [70], loss of ID2 expression leads to decreased white adipose tissue development and adipogenesis, both in vitro and in vivo [71]. ITGB2 was an up-regulated adipokine associated to obesity [72]. And SH3PXD2B positively regulates the fat cell differentiation [73]. The role of DNA methylation on the expression of these genes requires further investigations.

In this study, we performed transcriptome comparison on LD muscle from four developmental stages, and DNA methylome comparison between two developmental stages, three individuals in each stage. Although a relatively high correlation was observed in each group, ideal analysis should be performed on muscle biopsies collected in the same pig at all developmental stages to eliminate the individual effect. Function of the IMF related genes identified by this study will be further characterized on cell level.

Conclusions

In summary, this study provides a comprehensive landscape of transcriptome of the LD muscle in Laiwu pigs ranging from 60 d to 400 d. We identified 127 DEGs related to lipid biosynthesis in 120 d vs 240 d, including genes FASN, ACACA, ACLY and SCD, and transcription factors EGR1, SREBF2, TP53 and TWIST1. The fat biosynthesis relevant genes FASN, PVALB, ID2, SH3PXD2B and EGR1, showed differences in DNA methylation in the gene body or intergenic regions from 120 d to 240 d in Laiwu pigs. The role of these genes in porcine IMF deposition needs further investigations.

Methods

Animal sampling and meat composition measurements

A total of 12 pure castrated male Laiwu pigs reared under similar environmental and feeding conditions in Laiwu Pig Conservation Center (Laiwu, Shandong) were randomly selected at four developmental stages: 60, 120, 240 and 400 days of age (n = 3). Each stage includes three individuals of similar body weights. These pigs were weighed and electrically stunned to ameliorate the suffering of the animals before death, then blood drawn and sampled from the LD muscle at the third lumbar vertebra, and vacuum-packed in a low-oxygen environment until subsequent IMF content measurement by Soxhlet petroleum ether extraction after sampling 24-h [74, 75]. Similarly, four DYL pigs of 120 days of age were randomly selected from Xingyue Pig Breeding Co. Ltd. (Taian, Shandong). The DYL pigs were weighed, electrically stunned, blood drawn and slaughtered, all efforts were made to minimize suffering. Immediately after slaughter, a sample from the LD muscle at the third lumbar vertebra was collected and frozen in liquid nitrogen for further study.

Total lipids were extracted from the LD samples using a benzene-petroleum ether (1:1) mixture. The lipids were directly methylated using 0.4 mol/L KOH methyl alcohol solution according to Demirel et al. [76]. Fatty acid methyl esters (FAME) were quantified using a gas chromatograph (GC) equipped with a DB-17 column from Agilent Technologies (30 m × 0.25 mm × 0.25 μm; Palo Alto, CA, USA) on a Shimadzu GC-2010 (Kyoto, Japan). The fatty acid composition was analyzed with an Automatic Amino Acid Analyzer, and the operating steps were performed as previously described [77].

Total RNA extraction and RNA-sequencing

Total RNA was extracted from Laiwu pig muscle samples (stored at −80 °C) according to the EZNA Tissue RNA Kit instruction manual (Omega-Biotech, Doraville, USA). RNA concentration and purity were measured with a NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE). RNA integrity was assessed using the RNA Nano 6000 Assay Kit with the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

RNA sequencing libraries were generated using the NEBNext Ultra RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA) with multiplexing primers, according to the manufacturer’s instructions. Sequencing was performed using a paired-end 125-cycle rapid run on the Illumina HiSeq2500 (Illumina Inc., San Diego, CA, USA). Low quality reads were removed by perl script, and the clean reads were filtered from the raw reads and mapped to the Sus scrofa genome (Sscrofa10.2) using Tophat2 software [78]. Gene expression levels were estimated based on the FPKM values obtained using Cufflinks software [79]. Here, only genes with an absolute value of log2 (Fold-Change) ≥ 2 and an FDR < 0.01 were used for subsequent analysis. The transcriptome data have been deposited with the NCBI Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) under accession number GSE90135. Additionally, the transcriptome data of IB × DU pigs at 120 d were downloaded from GEO of NCBI (GSE86441) to compared with the Laiwu pig at the same age [16]. The DEGs were also screened based on log2FC ≥ 2 and an FDR < 0.01.

DNA preparation and RRBS sequencing

The genomic DNAs were extracted using the TIANamp genomic DNA kit (Tiangen, Beijing, China) following the manufacturer’s instructions of the 120 and 240 days old pig muscle tissues, based on the phenotypic and transcriptome data showed that this stage was the IMF deposition fastest period. The construction of DNA methylation library was performed at Biomarker (Beijing, China) using RRBS. The RRBS library was constructed using 2 μg of high-quality genomic DNA, pooling with an equal amount from three 120 and 240 days old pigs, respectively. Briefly, DNA was restriction digested using MspI enzyme, which cut the DNA at sites CCGG, then the fragment were end-repaired and dA-tailing to blunt-end products, followed by adaptor-ligation with T overhang. The ligation products were purified by 2% agarose gel electrophoresis and size-selected of DNA fragments 150–400 bp long (including 100 bp adaptor). Size-selected DNA was bisulfite-conversion with the NEXTflex Bisulfite-Seq Kit (Bioo Scientific, Austin, TX, USA). The final library was generated by PCR-amplified, enriching for fragments with adapters on both ends and the RRBS was performed by Illumina HiSeq X Ten (Illumina Inc., San Diego, CA, USA). The clean reads were aligned to the pig reference genome (Sscrofa10.2) (ftp://ftp.ensembl.org/pub/release-75/fasta/sus_scrofa) and produced by the BS-seeker2 v.2.0.8 using Bowtie2 v.2.1.0 [80] in local alignment mode and no more than 2 mismatches per read. Methylation status was determined using weighted methylation level [81]. The DMRs were produced, which sites coverage depth more than 10×, and at least 3 different methylation sites and Fisher’s exact test p < 0.05 using MOABS [27]. The RRBS sequences were submitted to Gene Expression Omnibus (GEO) of NCBI under the study accession number GSE93563.

Bioinformatics analysis

The DEGs in 60 d vs 120 d, 120 d vs 240 d, and 240 d vs 400 d comparisons were uploaded to DAVID Bioinformatics Resources (version 6.8, https://david.ncifcrf.gov/home.jsp) to pathway enrichment analysis. STEM software (version 1.3.9, http://www.cs.cmu.edu/~jernst/stem) was used to classify all of the genes into profiles with similar expression patterns [24], and the correlation coefficients of significant clustered profiles were calculated using the MultiExperiment Viewer (version 4.9.0, http://www.tm4.org). All of the DEGs were analyzed using Cytoscape ClueGO plug-in (version 2.3.2; http://apps.cytoscape.org/apps/cluego) as a complementary analysis method. Only Benjamini-corrected values of p < 0.05 were considered statistically significant. The obtained IMF candidate genes were compared with pig IMF relevant QTL dataset to validate fat metabolism genes. Finally, the upstream regulators prediction was predicted by Ingenuity Pathway Analysis (Ingenuity Systems, Redwood City, CA, USA; http://www.ingenuity.com), a bioinformatics tool for biological functions, canonical pathways and gene networks. Meanwhile, the TFs binding sites were predicted using JASPAR (http://jaspar.genereg.net). The DMRs were blasted [82] and pathway enriched also using online DAVID analysis.

qRT-PCR

Quantitative PCR was performed to validate the DEGs according to the Pearson’s correlation between the FPKM from the RNA-Seq data and the relative expression data obtained using qRT-PCR. First-strand cDNA was synthesized using Primescript RT reagent (Takara Bio Inc., Otsu, Japan) in a 20 μl total volume, containing 1 μL total RNA, 1 μL gDNA Eraser, 2 μL 5 × gDNA Eraser Buffer, 4 μL 5 × Prime Script Buffer 2, 1 μL Prime Script RT Enzyme Mix, 1 μL RT Primer Mix4 and RNase-Free dH2O. The expression of six genes was quantified using SYBR Premix Ex Taq (Takara Bio Inc., Otsu, Japan) in Agilent Mx3000P system (Agilent Co., Wilmington, Delaware, USA) in a total volume of 25 μL containing 12.5 μL 2× Pre Ex Taq, 0.5 μL ROX II, 2 μL cDNA, 0.5 μL each of the forward and reverse primer (10 μM) and dH2O. The standard curve for each gene was used to confirm amplification specificity and efficiency, and an appropriate 3-fold dilution ratio was used. The data were normalized to GAPDH and RPL7, and the experiments were run in triplicate. The 2-∆∆CT method was used to calculate relative gene expression. Concordance correlation coefficient [83] was calculated between FPKM and qRT-PCR experiments.

Validation of RRBS data using BSP

We performed BSP to validate the RRBS results. The methylation primers were designed using Methyl Primer Express v1.0, which were provided in Additional file 18: Table S12. The bisulfite conversion of pig genome DNA was executed according to protocol of EZ DNA Methylation-Gold Kit (Zymo Research, Irvine, CA, USA). Then amplification converted DNA using Ex Taq Hot Start Version (Takara Bio Inc., Otsu, Japan). And the PCR product was cloned into the pMD18-T vector (Takara Bio Inc., Otsu, Japan). Twenty subclones were selected for each fragment and the positive clones were sequenced using ABI3730XL DNA Analyzer (ABI, CA, USA). All the sequences were analyzed using BiQ Analyzer v2.0.

Statistical analysis

The live weight, IMF content and fatty acid composition over four developmental stages were performed with one-way ANOVA, followed by Duncan’s multiple range test (p < 0.05) by using SAS 8.2 (SAS Institute, Cary, NC, USA). Student t-test was employed to analyze the mRNA over four developmental stages of pig muscle. The GO and KEGG analysis were performed by Fisher’s t-test. p < 0.05 was considered significantly different.

Abbreviations

BP:

Biological process

CC:

Cellular component

DEGs:

Differentially expressed genes

DMCs:

Differentially methylated cytosines

DMRs:

Differentially methylated regions

DYL:

Duroc × Yorshire × Landrace

FC:

Fold-change

FDR:

False discovery rate

FPKM:

Fragments Per Kilobase of transcript per Million mapped reads

GO:

Gene ontology

IMF:

Intramuscular fat

IPA:

Ingenuity Pathway Analysis

KEGG:

Kyoto encyclopedia of genes and genomes

LD:

longissimus dorsi

MF:

molecular function

MUFA:

Monounsaturated fatty acids

PUFA:

Polyunsaturated fatty acids

RRBS:

Reduced representation bisulfite sequencing

SFA:

Saturated fatty acids

STEM:

Short Time-series Expression Miner

TFs:

Transcription factors

TSS:

Transcription start sites

TTS:

Transcription termination site

References

  1. 1.

    Puig-Oliveras A, Ballester M, Corominas J, Revilla M, Estellé J, Fernández AI, Ramayo-Caldas Y, Folch JM. A co-association network analysis of the genetic determination of pig conformation, growth and fatness. PLoS One. 2014;9:e114862.

  2. 2.

    Ayuso M, Fernández A, Núñez Y, Benítez R, Isabel B, Barragán C, Fernández AI, Rey AI, Medrano JF, Cánovas Á, et al. Comparative analysis of muscle transcriptome between pig genotypes identifies genes and regulatory mechanisms associated to growth, fatness and metabolism. PLoS One. 2015;10:e145162.

  3. 3.

    Wu T, Zhang Z, Yuan Z, Lo LJ, Chen J, Wang Y, Peng J. Distinctive genes determine different intramuscular fat and muscle fiber ratios of the longissimus dorsi muscles in Jinhua and landrace pigs. PLoS One. 2013;8:e53181.

  4. 4.

    Hocquette JF, Gondret F, Baéza E, Médale F, Jurie C, Pethick DW. Intramuscular fat content in meat-producing animals: development, genetic and nutritional control, and identification of putative markers. Animal. 2010;4:303–19.

  5. 5.

    Gerbens F, Jansen A, van Erp AJ, Harders F, Meuwissen TH, Rettenberger G, Veerkamp JH, Te PM. The adipocyte fatty acid-binding protein locus: characterization and association with intramuscular fat content in pigs. Mamm Genome. 1998;9:1022–6.

  6. 6.

    Essen-Gustavsson B, Karlsson A, Lundstrom K, Enfalt AC. Intramuscular fat and muscle fibre lipid contents in halothane-gene-free pigs fed high or low protein diets and its relation to meat quality. Meat Sci. 1994;38:269–77.

  7. 7.

    Guo B, Kongsuwan K, Greenwood PL, Zhou G, Zhang W, Dalrymple BP. A gene expression estimator of intramuscular fat percentage for use in both cattle and sheep. J Anim Sci Biotechnol. 2014;5:35.

  8. 8.

    Ovilo C, Benitez R, Fernandez A, Nunez Y, Ayuso M, Fernandez AI, Rodriguez C, Isabel B, Rey AI, Lopez-Bote C, et al. Longissimus dorsi transcriptome analysis of purebred and crossbred Iberian pigs differing in muscle characteristics. BMC Genomics. 2014;15:413.

  9. 9.

    Wood JD, Enser M, Fisher AV, Nute GR, Sheard PR, Richardson RI, Hughes SI, Whittington FM. Fat deposition, fatty acid composition and meat quality: a review. Meat Sci. 2008;78:343–58.

  10. 10.

    Cameron ND, Enser M, Nute GR, Whittington FM, Penman JC, Fisken AC, Perry AM, Wood JD. Genotype with nutrition interaction on fatty acid composition of intramuscular fat and the relationship with flavour of pig meat. Meat Sci. 2000;55:187–95.

  11. 11.

    Blanchard PJ, Willis MB, Warkup CC, Ellis M. The influence of carcass backfat and intramuscular fat level on pork eating quality. J Sci Food Agr. 2000;80:145–51.

  12. 12.

    Fan B, ZQ D, Rothschild MF. The fat mass and obesity-associated (FTO) gene is associated with intramuscular fat content and growth rate in the pig. Anim Biotechnol. 2009;20:58–70.

  13. 13.

    Cho KH, Kim MJ, Jeon GJ, Chung HY. Association of genetic variants for FABP3 gene with back fat thickness and intramuscular fat content in pig. Mol Biol Rep. 2011;38:2161–6.

  14. 14.

    Yang H, Huang X, Zeng Z, Zhang W, Liu C, Fang S, Huang L, Chen C. Genome-wide association analysis for blood lipid traits measured in three pig populations reveals a substantial level of genetic heterogeneity. PLoS One. 2015;10:e131667.

  15. 15.

    Chen Q, Zeng Y, Wang H, Yang L, Yang Y, Zhu H, Shi Y, Chen W, Hu Y. Molecular characterization and expression analysis of NDUFS4 gene in m. Longissimus dorsi of Laiwu pig (Sus Scrofa). Mol Biol Rep. 2013;40:1599–608.

  16. 16.

    Cui JX, Chen W, Zeng YQ. Development of FQ-PCR method to determine the level of ADD1 expression in fatty and lean pigs. Genet Mol Res. 2015;14:13924–31.

  17. 17.

    Cui J, Chen W, Liu J, Xu T, Zeng Y. Study on quantitative expression of PPARgamma and ADRP in muscle and its association with intramuscular fat deposition of pig. Spring. 2016;5:1501.

  18. 18.

    Zhang S, Shen L, Xia Y, Yang Q, Li X, Tang G, Jiang Y, Wang J, Li M, Zhu L. DNA methylation landscape of fat deposits and fatty acid composition in obese and lean pigs. Sci Rep-UK. 2016;6:35063.

  19. 19.

    Li M, Wu H, Luo Z, Xia Y, Guan J, Wang T, Gu Y, Chen L, Zhang K, Ma J, et al. An atlas of DNA methylomes in porcine adipose and muscle tissues. Nat Commun. 2012;3:850.

  20. 20.

    Zhao Y, Li J, Liu H, Xi Y, Xue M, Liu W, Zhuang Z, Lei M. Dynamic transcriptome profiles of skeletal muscle tissue across 11 developmental stages for both Tongcheng and Yorkshire pigs. BMC Genomics. 2015;16:377.

  21. 21.

    Choi M, Lee J, Le MT, Nguyen DT, Park S, Soundrarajan N, Schachtschneider KM, Kim J, Park J, Kim J, et al. Genome-wide analysis of DNA methylation in pigs using reduced representation bisulfite sequencing. DNA Res. 2015;22:343–55.

  22. 22.

    Ramayo-Caldas Y, Mach N, Esteve-Codina A, Corominas J, Castello A, Ballester M, Estelle J, Ibanez-Escriche N, Fernandez AI, Perez-Enciso M, et al. Liver transcriptome profile in pigs with extreme phenotypes of intramuscular fatty acid composition. BMC Genomics. 2012;13:547.

  23. 23.

    Huang Y, Sun J, Zhang L, Li C, Womack JE, Li Z, Lan X, Lei C, Zhang C, Zhao X, et al. Genome-wide DNA methylation profiles and their relationships with mRNA and the microRNA transcriptome in bovine muscle tissue (Bos taurine). Sci Rep-UK. 2014;4:6546.

  24. 24.

    Ernst J, Bar-Joseph Z. STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics. 2006;7:191.

  25. 25.

    Kramer A, Green J, Pollard JJ, Tugendreich S. Causal analysis approaches in ingenuity pathway analysis. Bioinformatics. 2014;30:523–30.

  26. 26.

    Ayuso M, Fernandez A, Nunez Y, Benitez R, Isabel B, Fernandez AI, Rey AI, Gonzalez-Bulnes A, Medrano JF, Canovas A, et al. Developmental stage, muscle and genetic type modify muscle transcriptome in pigs: effects on gene expression and regulatory factors involved in growth and metabolism. PLoS One. 2016;11:e167858.

  27. 27.

    Sun D, Xi Y, Rodriguez B, Park HJ, Tong P, Meong M, Goodell MA, Li W. MOABS: model based analysis of bisulfite sequencing data. Genome Biol. 2014;15:R38.

  28. 28.

    Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57.

  29. 29.

    Munoz G, Ovilo C, Noguera JL, Sanchez A, Rodriguez C, Silio L. Assignment of the fatty acid synthase (FASN) gene to pig chromosome 12 by physical and linkage mapping. Anim Genet. 2003;34:234–5.

  30. 30.

    Sevane N, Armstrong E, Cortés O, Wiener P, Wong RP, Dunner S. Association of bovine meat quality traits with genes included in the PPARG and PPARGC1A networks. Meat Sci. 2013;94:328–35.

  31. 31.

    Davoli R, Braglia S, Zappaterra M, Redeghieri C, Comella M, Zambonelli P. Association and expression analysis of porcine ACLY gene related to growth and carcass quality traits in Italian large white and Italian Duroc breeds. Livest Sci. 2014;165:1–7.

  32. 32.

    Lee H, Park WJ. Unsaturated fatty acids, desaturases, and human health. J Med Food. 2014;17:189–97.

  33. 33.

    Xing K, Zhu F, Zhai L, Chen S, Tan Z, Sun Y, Hou Z, Wang C. Identification of genes for controlling swine adipose deposition by integrating transcriptome, whole-genome resequencing, and quantitative trait loci data. Sci Rep-UK. 2016;6:23219.

  34. 34.

    Zulkifli RM, Parr T, Salter AM, Brameld JM. Regulation of ovine and porcine stearoyl coenzyme a desaturase gene promoters by fatty acids and sterols. J Anim Sci. 2010;88:2565–75.

  35. 35.

    Bessa RJB, Hughes RA, Jeronimo E, Moreira OC, Prates JAM, Doran O. Effect of pig breed and dietary protein level on selected fatty acids and stearoyl-coenzyme a desaturase protein expression in longissimus muscle and subcutaneous fat. J Anim Sci. 2013;91:4540–6.

  36. 36.

    Malerba G, Schaeffer L, Xumerle L, Klopp N, Trabetti E, Biscuola M, Cavallari U, Galavotti R, Martinelli N, Guarini P, et al. SNPs of the FADS gene cluster are associated with polyunsaturated fatty acids in a cohort of patients with cardiovascular disease. Lipids. 2008;43:289–99.

  37. 37.

    Ačimovič J, Rozman D. Steroidal triterpenes of cholesterol synthesis. Molecules. 2013;18:4002–17.

  38. 38.

    Mazein A, Watterson S, Hsieh W, Griffiths WJ, Ghazal P. A comprehensive machine-readable view of the mammalian cholesterol biosynthesis pathway. Biochem Pharmacol. 2013;86:56–66.

  39. 39.

    Xing K, Zhu F, Zhai L, Liu H, Wang Y, Wang Z, Chen S, Hou Z, Wang C. Integration of transcriptome and whole genomic resequencing data to identify key genes affecting swine fat deposition. PLoS One. 2015;10:e122396.

  40. 40.

    Stupka N, Kintakas C, White JD, Fraser FW, Hanciu M, Aramaki-Hattori N, Martin S, Coles C, Collier F, Ward AC, et al. Versican processing by a disintegrin-like and metalloproteinase domain with thrombospondin-1 repeats proteinases-5 and -15 facilitates myoblast fusion. J Biol Chem. 2013;288(3):1907–17.

  41. 41.

    Cao Y. Angiogenesis and vascular functions in modulation of obesity, adipose metabolism, and insulin sensitivity. Cell Metab. 2013;18:478–89.

  42. 42.

    Lim KS, Lee KT, Park JE, Chung WH, Jang GW, Choi BH, Hong KC, Kim TH. Identification of differentially expressed genes in longissimus muscle of pigs with high and low intramuscular fat content using RNA sequencing. Anim Genet. 2016; doi:10.1111/age.12518.

  43. 43.

    Zappaterra M, Deserti M, Mazza R, Braglia S, Zambonelli P, Davoli R. A gene and protein expression study on four porcine genes related to intramuscular fat deposition. Meat Sci. 2016;121:27–32.

  44. 44.

    Canovas A, Estany J, Tor M, Pena RN, Doran O. Acetyl-CoA carboxylase and stearoyl-CoA desaturase protein expression in subcutaneous adipose tissue is reduced in pigs selected for decreased backfat thickness at constant intramuscular fat content. J Anim Sci. 2009;87:3905–14.

  45. 45.

    Gallardo D, Quintanilla R, Varona L, Diaz I, Ramirez O, Pena RN, Amills M. Polymorphism of the pig acetyl-coenzyme a carboxylase alpha gene is associated with fatty acid composition in a Duroc commercial line. Anim Genet. 2009;40:410–7.

  46. 46.

    Henriquez-Rodriguez E, Tor M, Pena RN, Estany J. A polymorphism in the stearoyl-CoA desaturase gene promoter increases monounsaturated fatty acid content in dry-cured ham. Meat Sci. 2015;106:38–43.

  47. 47.

    Serão NVL, Veroneze R, Ribeiro AMF, Verardo LL, Braccini Neto J, Gasparino E, Campos CF, Lopes PS, Guimarães SEF. Candidate gene expression and intramuscular fat content in pigs. J Anim Breed Genet. 2011;128:28–34.

  48. 48.

    Shen Q, Bai Y, Chang KC, Wang Y, Burris TP, Freedman LP, Thompson CC, Nagpal S. Liver X receptor-retinoid X receptor (LXR-RXR) heterodimer cistrome reveals coordination of LXR and AP1 signaling in keratinocytes. J Biol Chem. 2011;286:14554–63.

  49. 49.

    Boyle KB, Hadaschik D, Virtue S, Cawthorn WP, Ridley SH, O'Rahilly S, Siddle K. The transcription factors Egr1 and Egr2 have opposing influences on adipocyte differentiation. Cell Death Differ. 2009;16:782–9.

  50. 50.

    Yang HW, Xia T, Chen ZL, Feng SQ, Peng Y, Zhou L, Gan L, Yang ZQ. Cloning, chromosomal localization and expression patterns of porcine Kruppel-like factors 4, −5, −7 and the early growth response factor 2. Biotechnol Lett. 2006;29:157–63.

  51. 51.

    Heckmann BL, Zhang X, Xie X, Saarinen A, Lu X, Yang X, Liu J. Defective adipose lipolysis and altered global energy metabolism in mice with adipose overexpression of the lipolytic inhibitor G0/G1 switch gene 2 (G0S2). J Biol Chem. 2014;289:1905–16.

  52. 52.

    Lin WH, Chiu KC, Chang HM, Lee KC, Tai TY, Chuang LM. Molecular scanning of the human sorbin and SH3-domain-containing-1 (SORBS1) gene: positive association of the T228A polymorphism with obesity and type 2 diabetes. Hum Mol Genet. 2001;10:1753–60.

  53. 53.

    Yang WS, Lee WJ, Huang KC, Lee KC, Chao CL, Chen CL, Tai TY, Chuang LM. MRNA levels of the insulin-signaling molecule SORBS1 in the adipose depots of nondiabetic women. Obes Res. 2003;11:586–90.

  54. 54.

    Brey CW, Nelder MP, Hailemariam T, Gaugler R, Hashmi S. Kruppel-like family of transcription factors: an emerging new frontier in fat biology. Int J Biol Sci. 2009;5:622–36.

  55. 55.

    Chen Z, Ding Z, Ma G, Liu N, Qian Q. Sterol regulatory element-binding transcription factor (SREBF)-2, SREBF cleavage-activating protein (SCAP), and premature coronary artery disease in a Chinese population. Mol Biol Rep. 2011;38:2895–901.

  56. 56.

    Verardo LL, Nascimento CS, Silva FF, Gasparino E, Martins MF, Toriyama E, Faria VR, Botelho ME, Costa KA, Lopes PS, et al. Identification and validation of differentially expressed genes from pig skeletal muscle. J Anim Breed Genet. 2013;130:372–81.

  57. 57.

    Kim SS, Kim JR, Moon JK, Choi BH, Kim TH, Kim KS, Kim JJ, Lee CK. Transcriptional alteration of p53 related processes as a key factor for skeletal muscle characteristics in Sus Scrofa. Mol Cells. 2009;28:565–73.

  58. 58.

    Ma W, Lu S, Sun T, Wang X, Ma Y, Zhang X, Zhao R, Wang Y. Twist 1 regulates the expression of PPARgamma during hormone-induced 3T3-L1 preadipocyte differentiation: a possible role in obesity and associated diseases. Lipids Health Dis. 2014;13:132.

  59. 59.

    Ren R, Chen Z, Zhao X, Sun T, Zhang Y, Chen J, Lu S, Ma W. A possible regulatory link between twist 1 and PPARgamma gene regulation in 3T3-L1 adipocytes. Lipids Health Dis. 2016;15:189.

  60. 60.

    Zhou J, Sears RL, Xing X, Zhang B, Li D, Rockweiler NB, Jang HS, Choudhary MNK, Lee HJ, Lowdon RF, et al. Tissue-specific DNA methylation is conserved across human, mouse, and rat, and driven by primary sequence conservation. BMC Genomics. 2017;18:724.

  61. 61.

    Sati S, Tanwar VS, Kumar KA, Patowary A, Jain V, Ghosh S, Ahmad S, Singh M, Reddy SU, Chandak GR, et al. High resolution methylome map of rat indicates role of intragenic DNA methylation in identification of coding region. PLoS One. 2012;7:e31621.

  62. 62.

    Kass SU, Landsberger N, Wolffe AP. DNA methylation directs a time-dependent repression of transcription initiation. Curr Biol. 1997;7:157–65.

  63. 63.

    Ball MP, Li JB, Gao Y, Lee JH, LeProust EM, Park IH, Xie B, Daley GQ, Church GM. Targeted and genome-scale strategies reveal gene-body methylation signatures in human cells. Nat Biotechnol. 2009;27:361–8.

  64. 64.

    Laurent L, Wong E, Li G, Huynh T, Tsirigos A, Ong CT, Low HM, Kin SK, Rigoutsos I, Loring J, et al. Dynamic changes in the human methylome during differentiation. Genome Res. 2010;20:320–31.

  65. 65.

    Cokus SJ, Feng S, Zhang X, Chen Z, Merriman B, Haudenschild CD, Pradhan S, Nelson SF, Pellegrini M, Jacobsen SE. Shotgun bisulphite sequencing of the Arabidopsis genome reveals DNA methylation patterning. Nature. 2008;452:215–9.

  66. 66.

    Klose RJ, Bird AP. Genomic DNA methylation: the mark and its mediators. Trends Biochem Sci. 2006;31:89–97.

  67. 67.

    Lorincz MC, Dickerson DR, Schmitt M, Groudine M. Intragenic DNA methylation alters chromatin structure and elongation efficiency in mammalian cells. Nat Struct Mol Biol. 2004;11:1068–75.

  68. 68.

    Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat Rev Genet. 2012;13:484–92.

  69. 69.

    Xu Y, Wu F, Tan L, Kong L, Xiong L, Deng J, Barbera AJ, Zheng L, Zhang H, Huang S, et al. Genome-wide regulation of 5hmC, 5mC, and gene expression by Tet1 hydroxylase in mouse embryonic stem cells. Mol Cell. 2011;42:451–64.

  70. 70.

    Murad JM, Place CS, Ran C, Hekmatyar SKN, Watson NP, Kauppinen RA, Israel MA. Inhibitor of DNA binding 4 (ID4) regulation of adipocyte differentiation and adipose tissue formation in mice. J Biol Chem. 2010;285:24164–73.

  71. 71.

    Park KW, Waki H, Villanueva CJ, Monticelli LA, Hong C, Kang S, MacDougald OA, Goldrath AW, Tontonoz P. Inhibitor of DNA binding 2 is a small molecule-inducible modulator of peroxisome proliferator-activated receptor-gamma expression and adipocyte differentiation. Mol Endocrinol. 2008;22:2038–48.

  72. 72.

    Gil A, María Aguilera C, Gil-Campos M, Cañete R. Altered signalling and gene expression associated with the immune system and the inflammatory response in obesity. Brit. J Nutr. 2007;98:S121–6.

  73. 73.

    Fan C, Dong H, Yan K, Shen W, Wang C, Xia L, Zhan D, Qi K. Genome-wide screen of promoter methylation identifies novel markers in diet-induced obese mice. Nutr Hosp. 2014;30:42–52.

  74. 74.

    de Koning DJ, Janss LL, Rattink AP, van Oers PA, de Vries BJ, Groenen MA, van der Poel JJ, de Groot PN, Brascamp EW, van Arendonk JA. Detection of quantitative trait loci for backfat thickness and intramuscular fat content in pigs (Sus Scrofa). Genetics. 1999;152:1679–90.

  75. 75.

    Gerbens F, Verburg FJ, Van Moerkerk HT, Engel B, Buist W, Veerkamp JH, Te PM. Associations of heart and adipocyte fatty acid-binding protein gene expression with intramuscular fat content in pigs. J Anim Sci. 2001;79:347–54.

  76. 76.

    Demirel G, Wachira AM, Sinclair LA, Wilkinson RG, Wood JD, Enser M. Effects of dietary n-3 polyunsaturated fatty acids, breed and dietary vitamin E on the fatty acids of lamb muscle, liver and adipose tissue. Br J Nutr. 2004;91:551–65.

  77. 77.

    Jiang YZ, Zhu L, Li XW, Si T. Evaluation of the Chinese indigenous pig breed Dahe and crossbred Dawu for growth and carcass characteristics, organ weight, meat quality and intramuscular fatty acid and amino acid composition. Animal. 2011;5:1485–92.

  78. 78.

    Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14:R36.

  79. 79.

    Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28:511–5.

  80. 80.

    Schachtschneider KM, Madsen O, Park C, Rund LA, Groenen MAM, Schook LB. Adult porcine genome-wide DNA methylation patterns support pigs as a biomedical model. BMC Genomics. 2015;16:743.

  81. 81.

    Schultz MD, Schmitz RJ, Ecker JR. Leveling’ the playing field for analyses of single-base resolution DNA methylomes. Trends Genet. 2012;28:583–5.

  82. 82.

    Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25:3389–402.

  83. 83.

    Miron M, Woody OZ, Marcil A, Murie C, Sladek R, Nadon R. A methodology for global validation of microarray experiments. BMC Bioinformatics. 2006;7:333.

Download references

Acknowledgements

The authors would like to thank Mr. Yanxiao Sun and Dr. Yujun Gao for providing the Laiwu pigs and DYL pigs. The authors would also like to thank Biomarker of Beijing for their bioinformatics support contributions.

Funding

This research was financially supported by grants from the Natural Science Foundation of Shandong Province (ZR2012CZ002), the Agricultural Breed Project of Shandong Province (2016), the Application of Agricultural Technology of Shandong Province (2013), and Shandong “Double Tops” Program (SYL2017YSTD12). The funding body had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Availability of data and materials

The datasets generated and/or analysed during the current study are available in the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) repository, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE90135,

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE93563.

Author information

Affiliations

  1. Shandong Provincial Key Laboratory of Animal Biotechnology and Disease Control and Prevention, College of Animal Science and Veterinary Medicine, Shandong Agricultural University, No. 61 Daizong Street, Taian, 271018, People’s Republic of China

    • Yuding Wang
    • , Cai Ma
    • , Yi Sun
    • , Li Kang
    •  & Yunliang Jiang

  2. Shandong Provincial Key Laboratory of Animal Disease Control and Breeding, Institute of Animal Science and Veterinary Medicine, Shandong Academy of Agricultural Sciences, Jinan, 250100, People’s Republic of China

    • Yuding Wang

  3. Central Hospital of Taian, Taian, 271018, People’s Republic of China

    • Yi Li

Authors

  1. Search for Yuding Wang in:

  2. Search for Cai Ma in:

  3. Search for Yi Sun in:

  4. Search for Yi Li in:

  5. Search for Li Kang in:

  6. Search for Yunliang Jiang in:

Contributions

YW and YJ designed and drafted the manuscript. YW and CM carried out animal care, prepared samples and performed the experiments. YW and LK performed the data processing and biological information analysis. YJ, YS and YL conceived the study and the experimental design and helped draft the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Yunliang Jiang.

Ethics declarations

Ethics approval and consent to participate

The animal experiments were carried out in accordance with the protocols of the ‘Guidelines for Experimental Animals’ of the Ministry of Science and Technology (Beijing, China) and all efforts were made to minimize suffering. The animal experiments were approved by the Institutional Animal Care and Use Ethics Committee of Shandong Agricultural University (Permit Number: 2,007,005). We have obtained the informed consents of a written form from Mr. Yanxiao Sun and Dr. Yujun Gao, who was the owner of the Laiwu pigs and DYL pigs, respectively.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1: Figure S1.

Validation of the expression of candidate genes using qRT-PCR. The Y-axis on the left side of the histogram represents the gene expression level according to qRT-PCR (marked as *), and the right Y-axis of line represents the standard value of FPKM based on transcription (indicated as #). The Pearson’s correlation coefficient (r) and Gene Symbol are shown above the figure. #/*p < 0.05; ##/**p < 0.01; ###/***p < 0.001. (PNG 2842 kb)

Additional file 2: Table S1.

The four significantly clustered profiles and the corresponding genes across four developmental stages in Laiwu pigs. A total of 10 profiles were generated and only four profiles significantly clustered across 60, 120, 240 and 400 days old in Laiwu pigs. (XLSX 76 kb)

Additional file 3: Figure S2.

The differentially expressed genes and functions associated with muscle metabolism in profile 1, 3 and 4. GO terms and genes are represented as nodes based on their kappa score more than 0.4 and networks with at least three nodes. The node size represents the GO terms enrichment significance. (PNG 1770 kb)

Additional file 4: Table S2.

The differentially expressed genes (DEGs) in 60 d vs 120 d, 120 d vs 240 d and 240 d vs 400d. We showed all DEGs only with absolute value of log2 (Fold-Change) ≥ 1 and FDR ≤ 0.01 in adjacent development stages of Laiwu pigs. (XLSX 81 kb)

Additional file 5: Table S3.

The fat relevant differentially expressed genes (DEGs) with similar expression trend in the comparison 120 d vs 240 d of Laiwu pigs. (XLSX 12 kb)

Additional file 6: Table S4.

The upstream regulators and their target genes in 120 d vs 240 d Laiwu pigs. (XLSX 11 kb)

Additional file 7: Table S5.

The TFs binding sites predicted by JASPAR database. (DOCX 280 kb)

Additional file 8: Figure S3.

The IMF content between DYL and Laiwu pig breeds. *p < 0.05. (PNG 5 kb)

Additional file 9: Table S6.

All the DEGs between IB × DU and Laiwu pig breeds. The DEGs were screened with absolute value of log2 (Fold-Change) ≥ 1 and FDR ≤ 0.01. (XLSX 281 kb)

Additional file 10: Table S7.

The IMF candidate genes between IB × DU and Laiwu pig breeds. (DOCX 17 kb)

Additional file 11: Figure S4.

Distribution of CpG methylation in the LD muscle of 120 d (a) and 240 d (b) Laiwu pigs across pig chromosomes. (PNG 1168 kb)

Additional file 12: Table S8.

The differentially methylated cytosines (DMCs) and corresponding genes in 120 and 240 days old Laiwu pigs. The DMCs were produced by sequencing depth ≥ 4×, and FDR ≤ 0.05. (XLSX 6864 kb)

Additional file 13: Table S9.

The differentially methylated regions (DMRs) in 120 and 240 days old Laiwu pigs. The DMRs were generated by sequencing depth ≥ 10×, at least 3 DMCs, meanwhile the minimum value ≥ 0.3, and p < 0.05 using Fisher’s exact test. (XLSX 461 kb)

Additional file 14: Table S10.

The classification of transcriptome according to their FPKM value at 120 d and 240 d of Laiwu pigs. (XLSX 2057 kb)

Additional file 15: Figure S5.

The ANOVA test among the highest, medium high, medium low and lowest FPKM groups in 120 d vs 240 d Laiwu pigs. ***p < 0.001. (PNG 10 kb)

Additional file 16: Figure S6.

The proportion of DEGs with DMRs, DEGs without DMRs and DMRs without DEGs in 120 d vs 240 d Laiwu pigs. (JPEG 75 kb)

Additional file 17: Table S11.

Description of DEGs with DMRs, DEGs without DMRs and DMRs without DEGs in 120 d vs 240 d Laiwu pigs. (XLSX 1094 kb)

Additional file 18: Table S12.

The primers for methylation analysis on intergenic region of EGR1 and promoter region of FASN by BSP. (DOCX 14 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Keywords

  • Pig
  • longissimus dorsi muscle
  • Intramuscular fat content
  • Transcriptome
  • DNA methylome
  • Gene

BMC Genomics

ISSN: 1471-2164

Contact us