Genome-Wide Analysis of mRNAs and lncRNAs of Intramuscular Fat Related to Lipid Metabolism in Two Pig Breeds

Intramuscular fat deposition is closely related to pork meat quality and appears to vary in different pig breeds. By selection of pigs with higher proportion of lean meat and higher growth rate, porcine backfat thickness is becoming thinner, but this has a negative effect on the intramuscular fat (IMF) levels in pig muscle [1]. In addition, human health problems caused by obesity are attracting more attention. It has been demonstrated that intramuscular adipocytes prefer to use glucose and obesity can lead to intramuscular lipid accumulation that is associated with several metabolic disorders including insulin resistance and type II diabetes [2, 3]. Importantly, pig is anatomically and physiologically similar to human, thus it can be used as a good medical model for studying lipid metabolic disease [4-6]. Clarifying the molecular basis of fat deposition in pigs will deepen our understanding of diseases associated with fat metabolism.

Previous studies have shown that several genes might regulate intramuscular adipogenesis and lipid metabolism in pigs. They include stearoyl-CoA desaturase (SCD) [7-9], fatty acid synthase (FASN) [10, 11], sterol regulatory element binding transcription factor 1 (SREBF1) [10], peroxisome proliferator-activated receptor γ (PPARγ) [12]. Some molecular events in the development of intramuscular fat have been recently revealed using longissimus dorsi muscle or adipocytes in pig [11, 13-15]. These studies indicate that adipocyte proliferation and differentiation is higher in subcutaneous fat than in intramuscular fat as well as in fat-type pig than in lean-type pig. Furthermore, the number and size of intramuscular adipocytes in fat-type pig are larger than those in lean-type pig. However, the regulatory mechanisms have not yet been fully elucidated. LncRNAs, non-protein coding transcripts longer than 200 nucleotides, have important regulatory functions in adipogenesis and lipid metabolism. With the development of high throughput sequencing technology, RNA-seq has been widely applied in the identification and functional analysis of lncRNAs and mRNAs of adipose tissues in pig [7, 13, 16, 17], chicken [18], sheep [19, 20] and cattle [21]. Ramayo et al. [22] identified 270 and 186 lncRNAs in the liver tissues of female pigs with high and low intramuscular fatty acid contents, respectively. Zhou et al [14]. analyzed the methylation of lncRNAs and found the expression level of linc-sscg3623 varied with pig breeds and developmental stages, affecting fat synthesis. By identifying the lncRNA and mRNA expression profiles of preadipocytes at different stages in abdominal adipose tissue of Jinghai Yellow chicken, Zhang et al. found that the target genes of lncRNAs were significantly enriched in MAPK and PPAR signaling pathway associated with adipocyte differentiation [18]. These studies indicate that lncRNAs might determine fat deposition and fatty acid composition, and regulate adipogenic differentiation and lipid metabolism in livestock and poultry. However, there is limited research on the expression profiles and the functions of lncRNAs in intramuscular adipose tissue of pig breeds with significant differences in fat deposition.

Laiwu (LW) pig, a native breed prevalent in north China, is characterized by high quality meat with bright color, high water-holding capacity and especially high intramuscular fat content (10.32%) [23]. On the contrary, Large White (LY) pigs, the most widely distributed lean-type pig breed, are faster-growing pigs with relatively low subcutaneous and intramuscular fat contents in comparison with Erhualian, Laiwu and Lulai Black pigs [15, 24]. Thus, LW and LY pigs provide good resources for studying intramuscular fat deposition, adipogenic differentiation and diseases associated with fat metabolism.

In this study, high-throughput RNA-seq was used to comparatively analyze the gene expression profiles of intramuscular adipose tissues in LW and LY pigs to identify key lncRNAs and genes associated with adipogenesis and lipid metabolism. Furthermore, by Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and Co-expression analyses, we investigated the molecular mechanism of differentially expressed lncRNAs and genes regulating fat deposition. The results can provide useful information for studying fat deposition regulated by lncRNAs in porcine intramuscular adipose tissue.

All the procedures involving animals were approved by the animal care and use committee at the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences, where the experiment was conducted. All the experiments were performed in accordance with the relevant guidelines and regulations set by the Ministry of Agriculture of the People’s Republic of China.

The experimental animals used here included 3 female LW pigs and 3 female LY pigs. All pigs were raised in Daqian Farming Co. Ltd. of Laiwu (Laiwu city, Shangdong Province, China) and fed a diet formulated to meet current nutritional requirements (National Research Council, NRC, 1998) [25] when these pigs were weaned at about 1 month of age with an average weight of about 5.5 kg. At slaughter age (180-day-old), the mean body weights of LW pigs and LY pigs were 40 kg and 108 kg, respectively. There are some differences in growth and carcass quality among different pig breeds. As a result, different pig breeds possess a different slaughter weight at slaughter age. On the other hand, the metabolic reaction is at equilibrium in the slaughter period. Therefore, we sampled at slaughter age, which is more representative. The intramuscular fat tissues were obtained from the longissimus dorsi muscle immediately after the pigs were slaughtered and then deeply frozen at -80 °C until use.

Total RNAs were extracted from the intramuscular adipose tissues not separated from longissimus dorsi samples of individual LW and LY pigs using mirVanaTM RNA Isolation Kits (#AM1561, Ambion, USA) according to the manufacturer’s instructions. The quality and quantity of the RNA samples were analyzed by Bioanalyzer 2100 system using an RNA 6000 Nano kit (Agilent technologies, Inc., Santa Clara, USA). The isolated RNAs were treated with RNase–free DNase I (Ambion, Inc., Austin, USA) to avoid any potential genomic DNA contamination. Then, approximately 2 μg of total RNAs from each sample was used for the construction of cDNA libraries (including LW1, LW2 and LW3; LY1, LY2 and LY3) according to Illumina® TruSeq^TM RNA Sample Preparation protocol. Agilent DNA 1000 kit on a Bioanalyzer 2100 (Agilent technologies, Inc., Santa Clara, USA) was applied to check the size and purity of each cDNA library. These libraries were sequenced on the platform of Illumina HiSeq 2500 with 100 paired-end sequencing (Illumina, CA, USA).

After removing adaptor contaminants and low-quality reads from the raw reads using Fastx_toolkit software (v0.0.14), the remaining sequences were used for further analysis. Clean reads of each sample were then aligned to reference genome Sscrofa10.2 (ftp://ftp.ensembl.org/pub/release-87/fasta/sus_scrofa/dna/) with Tophat2 software (v2.0.12) [26, 27]. The resulting alignment data from Tophat2 were then fed to Cufflinks (v2.1.1) [28, 29] for transcriptome assembly and reconstruction to obtain a non-redundant unified set of transcripts. Subsequently, these transcripts were aligned to reference annotation files Sscrofa10.2.87.chr.gtf (ftp://ftp.ensembl.org/pub/release-87/gtf/sus_scrofa) for identification of their types and positions. Transcript abundances were evaluated in Fragments Per Kilobase of exon per Million fragments mapped (FPKM), which eliminated the influence of sequencing depth, gene length and sample difference on gene expression level.

SNP detection was conducted with the SAMtools pileup program [30] to find variations in the SAM files and then VarScan program [31] to identify and filter variants based on read counts, base quality and allele frequency.

Potential lncRNA candidates were identified as follows. First, transcripts blasted to known mRNAs and belonging to other types of non-coding RNAs (including pseudogenes, pre-microRNA, tRNA, etc.) were removed. Then, transcripts with single exon and length less than 200 bp were excluded. PLEK [32], CNCI [33], CPC [34] and Pfam [35] were applied to assess the protein-coding potential of the above-obtained transcripts and the intersection of the results from these software tools was taken as the final result. In order to identify known lncRNAs, BLASTN tool was used to align lncRNA candidates to ALDB (A Domestic-Animal Long Noncoding RNA Database) [36], a database with a focus on the domestic-animal lncRNAs, with the settings of identity=100%, mismatch=0, E-value< 1e-10 and gap_opening=0.

Three biological replicates were set in the experiment. Differentially expressed lncRNAs and mRNAs were identified based on negative binomial distribution using R package DESeq2 [37]. Multiple hypothesis testing correction of P value was performed with Benjamini-Hochberg algorithm. If |log2FoldChange| ≥ 1 (LW vs LY) and padj (adjusted P value) ≤ 0.05, the genes and lncRNAs were considered to be differentially expressed.

GO (http://www.geneontology.org/) is the international standard classification of gene functions. It classifies gene functions along three aspects including molecular function, biological process and cellular component. KEGG (http://www.genome.jp/kegg) database [38] is a main public database for metabolic analysis and regulatory network research. To explore the main biological functions of differentially expressed genes, GO and KEGG signaling pathway enrichment analyses of mRNAs were performed using ClueGO [39] plugin of Cytoscape, based on hypergeometric distribution. Benjamini-Hochberg algorithm was used for P value correction. GO terms and pathways with Q value (corrected P value) ≤ 0.05 were considered as significantly enriched.

LncRNAs can cis-regulate neighboring target genes and trans-regulate distant target genes. Besides, genes with the same expression pattern might exhibit high correlation between biological functions. On the basis of above theories, the target genes of lncRNAs and key lncRNAs were identified as follows. For expression pattern analysis, Pearson correlation coefficient (PCC) was calculated for the expression values of each lncRNA and each mRNA, and co-expressed lncRNA-mRNA was selected with |PCC| > 0.8 and P value < 0.05. The genes transcribed within a 300-kb window upstream or downstream of lncRNAs were considered as cis target genes if the |PCC| of lncRNA-mRNA > 0.9. RNAplex [40] was applied to select trans-acting target genes according to the interaction between lncRNA and mRNA sequences with |PCC| > 0.9 and Energy value < -20 [41, 42]. GO annotation and KEGG pathway enrichment analysis were performed for these target genes to investigate the biological processes and signaling pathways that lncRNAs were mainly involved in and the functions of lncRNAs.

According to the results of correlation coefficient analysis as well as GO and KEGG analyses, the differentially expressed genes related to lipid metabolism, fat-associated diseases and adipocyte differentiation were selected (|PCCs| > 0.9) for co-expression construction to identify the key lncRNAs regulating fat deposition. Combined with cis- and trans-regulation analysis of lncRNAs, the total key lncRNAs were identified to serve as a basis for further research on lncRNA functions. GO annotation and KEGG enrichment analysis of the target genes of each lncRNA were performed, and the lncRNA-mRNA-signaling pathway network was constructed using Cytoscape [43] to further explore the functions of each key lncRNA.

qRT-PCR was performed to verify the expression level of genes. Eleven differentially expressed genes (five lncRNAs and six mRNAs) were randomly selected. Three biological replicates were set for each gene. GeneAmp® PCR System 9700 (Applied Biosystems, USA) was used to synthesize cDNA templates via reverse transcription. Approximately 0.5 μg of each RNA sample was used. QuantiFast® SYBR® Green PCR Kit (Qiagen, Germany) and LightCycler® 480 Ⅱ Real-time PCR Instrument (Roche, Swiss) were used for qRT-PCR analysis. The reaction system consisted of 1 μL of cDNA, 5 μL of 2× QuantiFast® SYBR® Green PCR Premix (Qiagen, Germany), 0.2 μL of forward primer, 0.2 μL of reverse primer and 3.6 μL of nuclease-free water. Reaction conditions were as follows: pre-denaturation at 95°C for 5 min, followed by 40 cycles; denaturation at 95°C for 10 s; renaturation at 60°C and extension for 30 s. Pig β-actin (ACTB) gene was used as the internal control [44] and 2^-△△Ct method was used to calculate the relative expression levels of genes between samples.

All the data were presented as “means ± SDs”. When comparisons were made, a Student’s t-test was performed and P < 0.05 was considered statistically significant.

By RNA-seq, approximately 11 Gb of sequence data per library were obtained. After pre-processing and low-quality trimming of the sequence data, about 88.9–91.6 million clean reads were obtained for each sample. Approximately 73.5%–76.3% of the reads were successfully aligned to the scrofa reference genome, and more than 60% reads were uniquely mapped (Table 1). A total of 92, 172 transcripts were obtained (FPKM ≥ 0.01), including 22, 712 known transcripts according to transcriptome assembly and reconstruction.

A total of 1, 048, 575 variations due to SNPs in each unigene were identified. Putative SNPs detected in uniquely mapped reference sequence were plotted along pig chromosomes. We found that several SNPs were mapped to some lipid metabolism-related genes, including PER1 Period Circadian Clock 1 (PER1), Microtubule Associated Protein Tau (MAPT) and Eukaryotic Translation Initiation Factor 4E Binding Protein 1 (EIF4EBP1).

A total of 1484 long intergenic noncoding RNAs (lincRNAs), 149 anti-sense lncRNAs, 136 sense lncRNAs and 97 intronic lncRNAs were identified in the six libraries (Fig. 1A). LincRNA was the most common type of lncRNA and its proportion was also the largest, which were consistent with current research (Fig. 1A). Among all the lncRNAs, 754 (about 40.41%) were mapped to ALDB database. The majority of lncRNAs had two exons (up to 31 exons, an average of 2.99 exons), which were significantly less than that of protein-coding genes (up to 118 exons, an average of 9.64 exons) (Fig. 1B). In addition, the amount of exons in lncRNAs in intramuscular adipose tissue was consistent with that in thyroid [45] and testicular tissue [46] of pig. Overall, the distribution of lncRNAs and mRNAs lengths were consistent, and the proportion of relatively long mRNA transcripts was higher than that of lncRNAs (Fig. 1C). The average lengths of lncRNAs and mRNAs were 2317 and 2229 nt, respectively. Violin plot of FPKM values for the transcripts indicated that the expression level of mRNAs was relatively higher than that of lncRNAs (Fig. 1D), which was consistent with the research results on pig [46], mouse [47] and human [48].

The Venn Diagram of the expressed lncRNAs and mRNAs showed that a total of 132 lncRNAs and 1161 mRNAs were expressed specifically in the intramuscular adipose tissue of LW pig, and 137 lncRNAs and 1088 mRNAs were expressed specifically in that of LY pig (Fig. 1E and F). A total of 55 differentially expressed lncRNAs (33 up-regulation and 22 down-regulation) (Table S1 - for all supplemental material see www.karger. com/10.1159/000495101/) and 715 differentially expressed mRNAs (371 up-regulation and 344 down-regulation), including 513 known mRNAs (Table S2), were identified.

GO enrichment analysis results showed that 210 genes were enriched in the biological process category, 144 genes were enriched in the molecular function category, and 62 genes were enriched in the cellular component category. Importantly, many of the significantly enriched GO terms were closely associated with lipid metabolism and adipogenesis (Table S3). As shown in Fig. 2, In the biological process category, many genes (≥15) were enriched in “lipid biosynthetic process”, “lipid metabolic process”, “cellular lipid metabolic process”, “response to lipid”, “MAPK cascade”, “positive regulation of MAPK cascade” and “regulation of MAPK cascade”. In terms of molecular function, only “enzyme inhibitor activity” was a significantly enriched GO term. In the biological process category, the genes were mainly enriched in GO terms associated with extracellular matrix and axon. GO annotation suggested that the regulatory mechanism and regulatory genes of intramuscular fat deposition in LW and LY pigs might be significantly different and closely related to different genes.

KEGG pathway enrichment analysis revealed that 143 genes were enriched in 94 signaling pathways, among which 31 signaling pathways were significantly enriched and involved 108 genes (Table S4). Specifically, many genes were significantly enriched in signaling pathways closely related to lipid metabolism and adipogenesis, including “MAPK signaling pathway”, “PPAR signaling pathway”, “Regulation of lipolysis in adipocytes”, “Glycerophospholipid metabolism”, “Insulin signaling pathway”, etc. (Fig. 3). These results were consistent with the GO analysis results. MAPK and PPAR signaling pathways were the key pathways regulating adipocyte differentiation and lipid accumulation. A total of 14 genes and 9 genes were enriched in the above two pathways respectively, implying that these genes might play critical roles in regulating intramuscular adipogenesis and lipid metabolism in LW and LY pigs.

According to the correlation coefficient analysis, 462 differentially expressed mRNAs were found to co-express with differentially expressed lncRNAs (Table S5). LncRNA-mRNA co-expression network indicated that 26 lncRNAs which were at the central positions of the network with Degree (Degree Centrality) ≥ 30 might play key roles in regulating intramuscular lipid accumulation (Table S6, Fig. 4). Cis-regulation analysis showed that few mRNAs was the target gene of differentially expressed lncRNAs (Table S7). Trans-regulation analysis showed that ten differentially expressed mRNAs related to lipid metabolism were trans-regulated by eight lncRNAs (Table S8, Table 2). According to the above three analyses, a total of 31 key lncRNAs were selected to investigate their functions in adipogenesis and lipid metabolism.

The target genes of lncRNAs were also significantly enriched in MAPK, PPAR signaling pathways and biological processes including lipid biosynthetic process, lipid metabolic process, fatty acid metabolism. This indicated that the differentially expressed lncRNAs might regulate adipogenic differentiation and lipid metabolism by targeting key genes involved in the above signaling pathways and biological processes. LncRNAs whose target genes were significantly enriched in MAPK and PPAR signaling pathways were selected for the construction of lncRNA-mRNA-signaling pathway network (Fig. 5). According to the network, the regulatory patterns of key lncRNAmRNA pairs were explored.

To validate the reliability of RNA-seq results, differentially expressed genes (five lncRNAs, six mRNAs) were randomly selected for further qRT-PCR verification (Table S9). qRT-PCR is considered as the golden standard for quantitative analysis of genes [49]. The results showed that FASN, TRIB3, MAPKAPK2, CD36, XLOC_046142, XLOC_053194 and XLOC_064871 were highly expressed in LY pig, while NR1D2, SCD, XLOC_027632 and XLOC_004938 were highly expressed in LW pig (Fig. 6). These data were consistent with the sequencing results, implying the reliability of the sequencing results.

Many researchers have used genome-wide RNA seq to study gene expression in pig intramuscular adipose tissue [50, 51]. In light of their works, we compared the expression profile of genes and lncRNAs in intramuscular adipose tissue between LW and LY pigs by RNA-seq technology and comprehensively investigated intramuscular fat deposition and metabolism. KEGG enrichment analysis indicated that the differentially expressed genes in intramuscular adipose tissues between LW and LY pigs were mainly enriched in MAPK and PPAR signaling pathways. MAPK pathway, one of the important signal systems mediating cells’ response to external stimuli, plays a key role in cell proliferation and differentiation. In mammalian cells, there are four subgroups of MAPK: ERK1/2, JNK, p38 MAPK and ERK5/ BMK1, which are important for regulating adipocyte proliferation and differentiation [52-54]. The precise understanding of the cascade of these molecular events and the way to regulate them will be crucial for efficiently treating obesity. PPAR signaling pathway is a key pathway closely associated with metabolism of fatty acid and sterols, and adipogenic differentiation [21, 55, 56]. Lim et al. used qRT-PCR to analyze the genes associated with PPAR signaling pathway in the longissimus dorsi muscle of Hanwoo cattle [57].

In this study, AFGF/FGF1, AKT2, GADD45G, HSP27, IL1R2, LOC100624812, MAP2K6, MAP3K12, MAP3K4, MAP3K6, MAPK10, MAPKAPK2, MAPT and PAK1 were significantly enriched in MAPK signaling pathway. Except for AFGF and MAPT, which were associated with ERK1/2 signaling pathway, all other genes were associated with p38 MAPK signaling pathway. Aouadi et al [54]. comprehensively studied the role of p38 MAKP signaling pathway in regulating fat deposition in mice. They suggested that p38 MAPK could inhibit adipogenic differentiation by inhibiting the activity and expression of C/EBPβ and PPARγ during the whole process of lipogenesis in mice. The activity of p38 MAPK was remarkably low in the white adipose tissue of obese mice (ob/ob mice, high-fat diet, hereditary obesity) and preadipocytes separated from obese mice. In addition, other research showed that p38 MAPK signal could negatively regulate preadipocyte differentiation as well as subcutaneous and intramuscular fat deposition in broiler chicken [58-60], and inhibit lipogenesis in 3T3-L1 adipocyte [61]. However, p38 MAKP could promote the differentiation of human preadipocytes [62]. Taken together, the role of p38 MAPK in regulating adipocyte differentiation and lipid deposition might be related to species, cell type and developmental stage.

Note that these genes involved in p38 MAPK signaling pathway were mainly up-regulated in the intramuscular adipose tissue of LY pig in this study. This implied that p38 MAPK signaling pathway might be important for intramuscular fat deposition in pig, which was consistent with our previous study showing that MAPK signaling pathway may play a critical role in porcine adipogenesis [17]. Fat deposition is a complex process related to lipogenesis, adipocyte proliferation and differentiation. Previous research showed that the proliferation, differentiation and lipogenesis capability of preadipocytes was stronger in Bamei (fat-type) pig than in Landrace (lean-type) pig [13]. The intramuscular fat content in Erhualian pig (fat-type) was remarkably different at different growth stages, whereas that in Large White pig changed little during growth [15]. Above research suggested that preadipocytes in obese pig breeds might have higher ability to proliferate and differentiate, thus more intramuscular fat was deposited in these pig breeds. This also implied that the intramuscular adipocytes in LY pig had weaker ability than those in LW pig to differentiate, proliferate and synthesize lipids. This might be a part of the reason that significantly less intramuscular fat was deposited in LY pig than in LW pig. On above basis, these genes enriched in p38 MAPK signaling pathway might mainly negatively regulate intramuscular adipogenic differentiation in LY pig and thus lead to deposition of less intramuscular fat. It was inferred that the 11 key lncRNAs such as XLOC_004398, XLOC_0053194, and XLOC_046142 could regulate the proliferation and differentiation of intramuscular adipocytes as well as lipogenesis in LW and LY pigs by target regulating these genes involved in p38 MAPK signaling pathway. The results were similar to those of Zhang et al., who found that the target genes of lncRNAs in abdominal fat tissue of chicken were significantly enriched in MAPK signaling pathway [18]. Note that XLOC_046142 might play a critical role in MAPK signaling pathway by trans-acting on MAPKAPK2.

In addition, the expression of AFGF (FGF1) and MAPT genes were highly expressed in the intramuscular adipose tissue of LW pig. In humans, FGF1 is a kind of lipogenesis factor secreted by adipose microvascular endothelial cells [63] and can stimulate and induce the proliferation and differentiation of human preadipocytes through FGF1/FGFR axis [63, 64]. In mice, FGF1 can act as a key signal transduction protein to mediate the reconstruction of adipose tissue and regulate metabolic homeostasis through PPARg-FGF1 axis [65]. Therefore, FGF1 might potentially promote the proliferation and differentiation of intramuscular adipocytes in LW pig and play an important role in intramuscular fat deposition.

In this study, SCP2, PLIN2, GK, CD36, ANGPTL4, SCD and CYP7A1 were significantly enriched in PPAR signaling pathway. Among them, SCD and CYP7A1 were highly expressed in LW pig.

SCD is the rate-limiting enzyme in the biosynthesis of unsaturated fatty acid. Besides, SCD gene is a candidate key gene regulating intramuscular fat deposition [8, 9]. It also plays an important role in regulating fatty acid metabolism and determining fatty acid composition in muscle [66, 67]. In this study, SCD was highly expressed in LW pig with high intramuscular fat content, which was consistent with current research. Compared with Duroc pig (lean type), LW pig has higher levels of intramuscular fat, saturated fatty acid (SFA) and monounsaturated fatty acid (MUFA) in its longissimus dorsi muscle [68]. Compared with Holstein cattle (low intramuscular fat level), Wagyu cattle (high intramuscular fat level) also has higher levels of MUFA in its intramuscular adipose tissue and longissimus dorsi muscle [69]. Since SFA and MUFA are key fatty acids determining meat quality, there should be significant difference in SAF and MUFA levels in adipose tissues between LW and LY pigs. The expression level of SCD is positively correlated with MUFA content [9], thus SCD should promote MUFA biosynthesis and deposition in the intramuscular adipose tissue of LW pig. In this study, lncRNAs including XLOC_061880, XLOC_025238 (PCC<-0.94) and XLOC_003916 (PCC> 0.94) ranked in the top three in terms of association with SCD expression. This suggested that these lncRNAs might be important for SCD expression. In sum, SCD might play a key role in intramuscular fat deposition and fatty acid biosynthesis in LW pig. Seven key lncRNAs including XLOC_025238, XLOC_027632 and XLOC_062192 might promote the intramuscular fat deposition and MUFA synthesis in LW pig by regulating genes, especially SCD gene, involved in PPAR signal pathway.

PLIN2 and FASN also play an important role in regulating intramuscular fat deposition [69, 11]. PLIN2, an adipose differentiation-related protein (ADRP), is a kind of perilipin which can promote the biosynthesis and storage of fatty acid by coating lipid droplets [69]. In this study, PLIN2 was down-regulated in the intramuscular adipose tissue of LW pig. This was inconsistent with previous research showing that PLIN2 was highly expressed in in the longissimus dorsi muscles of LW pig in comparison with that of LY pig [12]. FASN is a key enzyme of de novo fatty acid synthesis and can catalyze malonyl-CoA to form fatty acid [70]. The expression level of FASN was high in the longissimus dorsi muscle of Berkshire pig with high intramuscular fat content during the slaughter period [11]. In this study, FASN was highly expressed in the RNA extracted from intramuscular adipose tissue of LY pig and the intramuscular adipose tissue was separated in the part of longissimus muscle, which was inconsistent with the result of above-mentioned study.

By trans-regulation analysis and network analysis of lncRNAs and differentially expressed genes associated with lipid metabolism, it was found that TRIB3, NR1D2, BRCA1 and AKR1C4 were trans-regulated by XLOC_064871, XLOC_004398, XLOC_011001 and XLOC_015408, respectively. TRIB3, a kind of pseudokinase, can inhibit the synthesis of fatty acid by binding E3 ubiquitin ligase and protect against diet-induced obesity by stimulating fatty acid oxidation in the adipose tissue of mice during fasting [71]. In addition, TRIB3 can negatively regulate the transcriptional activation of PPARγ and inhibit adipogenesis in 3T3-L1 cells [72]. TRIB3 is expressed in porcine adipose tissue and skeletal muscle, and it is related to meat quality and production performance [73]. These indicate that TRIB3 might inhibit subcutaneous fat deposition in LY pig and be closely related to obesity-associated diseases. In this study, XLOC_064871 can trans-regulate TRIB3, indicating that XLOC_064871 might also play a vital role in adipocyte differentiation and fatty acid metabolism in pig. NR1D2 was highly expressed in skeletal muscle of Bantam chicken with high intramuscular fat content, and can induce the expression of PPARs, CEBPα and FABP4 genes, promoting adipogenic differentiation and affecting intramuscular fat deposition [74]. In this study, the expression of NR1D2 in intramuscular adipose tissue of LW pig was high, which was consistent with previous study. This indicated that NR1D2 might promote the adipogenic differentiation in pig and XLOC_004398 might play a role in intramuscular fat deposition in pig by interacting with NR1D2. GO annotation indicated that BRCA1 and AKR1C4 were mainly associated with fatty acid and lipid metabolism. Besides, previous study demonstrated that highly expressed BRCA1 in the subcutaneous fat of obese people might help to control fatty acid biosynthesis in adipocytes and adipose tissue from obese subjects [75], which was in according with our study that BRCA1 was highly expressed in LW pig. This implied that XLOC_011001 might target BRCA1, playing a vital role in fatty acid and lipid metabolism as well as obesity controlling.

In this study, we applied RNA-seq technology and bioinformatics methods to identify differentially expressed lncRNAs and mRNAs of intramuscular adipose tissues between Laiwu and Large White pigs for exploring the molecular mechanisms of intramuscular fat deposition. Results indicated that 55 lncRNAs and 513 known mRNAs were differentially expressed. LncRNAs can target mRNAs involved in PPAR and MAPK signal transduction pathways and play an important in intramuscular fat accumulation and adipogenic differentiation. XLOC_046142, XLOC_064871, XLOC_004398, XLOC_011001, XLOC_025238 which might target MAPKAPK2, TRIB3, NR1D2, BRCA1, and SCD, respectively, play key regulatory roles. This study provides useful information for understanding molecular mechanism of intramuscular fat deposition in the studied pig breeds.

X.Y.M. conceived and designed the study and wrote the paper. W.L.H performed the experiment, analyzed data and wrote the paper. A.L., X.X.Z. and L.L.X performed the experiment and interpreted the data. All the authors read and approved the final manuscript. This work was supported by a grant from the Agricultural Science and Technology Innovation Program (ASTIPIAS05), the National Basic Research and Development Program of China (973 Program) (No. 2015CB943100) and the Basic Research Fund for Central Public Research Institutes of CAAS (No. 2013ywf-zd-2).

The authors declare that they have no competing interests.