Abstract

Muscle growth and fat deposition are the two important biological processes in the development of pigs which are closely related to the pig production performance. Long intergenic noncoding RNAs (lincRNAs), with lack of coding potential and the length of at least 200nt, have been extensively studied to play important roles in many biological processes. However, the importance and molecular regulation mechanism of lincRNAs in the process of muscle growth and fat deposition in pigs are still to be further studied comprehensively. In our study, we used the data, including liver, abdominal fat, and longissimus dorsi muscle of 240 days’ age of two F2 full-sib female individuals from the white Duroc and Erhualian crossbreed, to identify 581 putative lincRNAs associated with pig muscle growth and fat deposition. The 581 putative lincRNAs shared many common features with other mammalian lincRNAs, such as fewer exons, lower expression levels, and shorter transcript lengths. Cross-tissue comparisons showed that many transcripts were tissue-specific and were involved in the important biological processes in their corresponding tissues. Gene ontology and pathway analysis revealed that many potential target genes (PTGs) of putative lincRNAs were involved in pig muscle growth and fat deposition-related processes, including muscle cell proliferation, lipid metabolism, and fatty acid degradation. In Quantitative Trait Locus (QTLs) analysis, some PTGs were screened from putative lincRNAs, MRPL12 is associated with muscle growth, GCGR and SLC25A10 were associated with fat deposition, and PPP3CA, DPYD, and FGGY were related not only to muscle growth but also to fat deposition. Therefore, it implied that these lincRNAs might participate in the biological processes related to muscle growth or fat deposition through homeostatic regulation of PTGs, but the detailed molecular regulatory mechanisms still needed to be further explored. This study lays the molecular foundation for the in-depth study of the role of lincRNAs in the pig muscle growth and fat deposition and further provides the new molecular markers for understanding the complex biological mechanisms of pig muscle growth and fat deposition.

1. Introduction

Muscle growth and fat deposition are complex quantitative traits and important economic traits in pig production. Muscle growth directly affects the meat quality and fat deposition of pigs that seriously affect pork quality and production efficiency. To a large extent, the two factors affect the consumer choices for pork [13]. In recent years, people have begun to pay more attention to pork quality rather than quantity. Duroc, a typical Western lean pig breed, has good characteristics of fast growth, high lean meat percentage, and high feed utilization rate [4]. Now it is widely used in commercial production. By contrary, Erhualian, a typical Chinese native dual-purpose pig breed, shows a good characteristic of high intramuscular fat content, which is delicate and juicy for processed foods [5]. Using these two pig breeds as parents, the F2 generations produced show different extreme phenotypes. Therefore, analyses of extreme performance individuals of the F2 generation provide a good opportunity to explore the molecular mechanism of muscle growth and fat deposition in these two breeds of pigs and further provide the new ideas for breeding pigs with better meat taste.

lincRNAs are a type of intergenic transcripts that are longer than 200bp in length and have almost no protein-coding ability [68]. Compared with protein-coding transcripts, lincRNAs have shorter transcript length, fewer exons number, stronger tissue specificity, and lower conservation [7, 9, 10]. So far, many lincRNAs have been discovered and identified for their corresponding functions. According to reports, there are 15,512 human lincRNAs and more than 10,000 mouse lincRNAs have been identified and some of them have been shown to play important roles in many biological processes [1113], such as gene regulation [14, 15], stem cell pluripotency [16, 17], X chromosome inactivation [18], and skeletal muscle development [7, 19]. However, compared with humans and mouse, many lincRNAs in pigs are still uncharacterized. The relationship between some lincRNAs and their potential target genes (PTGs) remains unclear too; meanwhile, there are relatively little understanding of the effects of lincRNAs on muscle growth and fat deposition.

In our experiment, we performed RNA-seq data analysis on the transcripts of liver, abdominal fat, and longissimus dorsi muscle in the F2 full-sib pig individuals of Duroc and Erhualian pigs by using the data that have been published on the NCBI. A total of 581 putative lincRNAs were identified and the basic features of these lincRNAs were also characterized. We obtained differentially expressed lincRNAs (DELs) from different tissues of two pigs by differential expression analysis. By conducting gene ontology and pathway analysis of the PTGs of lincRNAs, we found that some DELs were significantly involved in the regulation of muscle growth and fat deposition-related biological processes. Through QTLs analysis of DELs, it was shown that most of the DELs could positively regulate their PTGs expression. This result provides directions and goals for further studies on molecular mechanism underling pig muscle growth and fat deposition; meanwhile, it also provides valuable data for the modern biotechnology to breed new pig breeds with the good characteristics of fast growth as well as high intramuscular fat deposition.

2. Materials and Methods

2.1. Ethics Statement and Data Acquisition

The experimental protocols in our study were approved by the Ethics Committee of Huazhong Agricultural University in 2013 (HZAUMU2013-0005), and all the data were downloaded from the NCBI SRA website with the accession number provided by Wang et al. (Table 1). These datasets samples contained 6 RNA-seq samples, including liver, abdominal fat, and longissimus dorsi muscle of 240 days’ age of two F2 full-sib female individuals from the white Duroc and Erhualian crossbreed. Each of the 6 RNA-seq samples has 3 technical replicates, and the 3 samples are randomly taken from different parts of the same sample for sequencing. Then the technically repeated data are combined for subsequent analysis [20]. The pig gene annotations files were downloaded from http://ftp.ensemblorg.ebi.ac.uk/pub/release-91/gtf/sus_scrofa/, and the human and mouse lincRNAs references data were downloaded from http://asia.ensembl.org/info/data/ftp/index.html.

2.2. RNA-Seq Reads Mapping and Transcriptome Assembly

The quality of sequencing reads was evaluated by Fastqc command. The raw reads were filtered and trimed by Trimmomatic (version 0.36) with default parameters [21]. Then, clean reads were mapped to the pig reference genome (Sus scrofa 11.1) using HISAT2 (version 2.0.2) [22] with default parameters. The mapped reads were sorted and removed duplicates by SAMtools (version 0.1.19) [23]. Usually, StringTie (version 1.2.2) [22] was used to assemble the mapped reads with default parameters, and we set the “−G” option of StringTie for novel transcript assembly. Then, we used StringTie to merge the 6 assembled transcript files (GTF format) of tree groups into a nonredundant transcriptome.

2.3. Pipeline for lincRNA Identification

We processed nonredundant transcriptome through the lincRNAs detection pipeline to identify the lincRNAs, and the main steps are as follows (Figure 1): (1) filter the transcripts with ‘U’ category categorized (represent intergenic transcripts) by gffcompare program. (2) Reserve transcripts with exon >1 and length >200 bp. (3) The coding potential of the transcripts in both strands can be calculated by a coding potential calculator (CPC) tool [24], and the transcripts with cpc>0 were removed. (4) In order to evaluate whether the remaining transcripts contain any known protein-coding domain, the HMMER-3 was used to identify the transcripts translated in all six possible frames with homologs that were concluded in any of the known protein family domain in the Pfam database, and transcripts that matched to the Pfam hit (E-value < 1 × 10–5) were excluded. (5) The transcripts with the similarity to the known protein in against the NCBINR and UniRef90 database [25] would be filtered out by the BLASTX program with an E-value < 1 × 10–5. (6) Reserve the transcripts expressed in at least one sample [22, 26].

2.4. Differential Expression lincRNAs Analysis

We used “htseq-count” program to count the number of reads across six samples from three tissues, and then, in order to identify the effects of different lincRNAs from different tissues on muscle growth and fat deposition, we used the “DESeq2” package to screen out the differential lincRNAs between any two tissues [27]. In the screening results, transcripts with (FoldChange) greater than 1 or less than -1, and the corrected padj value less than 0.05, are considered to be the differential lincRNAs in these two groups [28].

2.5. Comparisons between lincRNAs and Protein-Coding Transcripts

We selected the transcripts annotated as “protein-coding” in gene annotation file, and the obtained lincRNAs were screened with “known” and “novel” by “blastn” command. Then we compared the lincRNAs (“known” and “novel”) with these protein-coding transcripts in the following aspects: transcript length, exon number and length, expression level, and FPKM.

2.6. Screen the Adjacent Genes of lincRNAs

Many of the lincRNAs are unknown, and the annotations are not comprehensive. To determine the relevant functions of linRNAs, first of all, functional prediction is required, and then further functional verification is performed. At present, the main prediction methods are cis-type prediction (adjacent gene method) and trans-type prediction (correlation prediction) [29, 30]. It has been reported that lincRNAs may have a greater regulatory effect on adjacent genes. In mammals, lincRNAs can regulate some biological processes such as development and transcriptional regulation through its adjacent target genes (<10kb) (e.g., GATA2, GZF1, and NEUROG2) and then participate in relevant biological processes [31]. Therefore, in our study, we used adjacent gene method to predict the function of lincRNAs. BEDTools (version 2.17.0) can be used to screen the adjacent gene (10k or 100k) of each lincRNAs locus. By analyzing the functions of their adjacent genes (such as pathway analysis and QTLs analysis), we could predict the function of corresponding lincRNAs [32]. In order to make the results better, in the process of late functional prediction, the selected adjacent genes needed to be expressed in at least one sample.

2.7. Gene Ontology and Pathway Analysis

In order to query each protein-coding gene and understand their functions, the DAVID database was used to perform gene function enrichment analysis via Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. And GO terms or KEGG pathways with corrected P-value less than 0.05 were considered to be enriched clusters. Due to the limitation of genes annotation in Sus scrofa, the BIOMART module in Ensembl was used to convert all genes into human homologues genes suitable for DAVID analysis [6].

2.8. QTLs Analysis of Differentially Expressed lincRNAs

To further predict the functions of all 581 putative lincRNAs, we selected 272 DELs from these putative lincRNAs to perform the QTLs analysis. The pig quantitative trait loci (QTLs) database was downloaded from the Animal QTLdb (Pig QTLdb) and the download path is https://www.animalgenome.org/cgi-bin/QTLdb/SS/download?file=bedSS_11.1. When doing QTL analysis, the main command was the “intersectBed” [33]. The p-value and Pearson correlation coefficient between each pair of lincRNAs and its PTGs obtained in RNA-seq were calculated with cor.test() function in the R software. When p-value < 0.05 [34, 35], it is considered to be statistically significant, and when the Pearson correlation coefficient is closer to 1, it indicates that the correlation between lincRNAs and its PTGs is higher [36].

2.9. Correlation Validation between lincRNAs and PTGs

We verified the quantitative relationship between lincRNAs and PTGs and selected 11 RNA samples from 3 tissues, including liver, abdominal fat, and longissimus dorsi muscle. For quantitative verification, total RNA was extracted using Trizol reagent (Invitrogen, Life Technologies, CA, USA) and performed according to the manufacturer's instructions. After checking the quality of RNA samples, cDNA synthesis for lincRNAs and PTGs detection was performed using the RevertAid First Strand cDNA Synthesis Kit (Thermo, Wuhan, Cat# k1622). According to the manufacturer's instructions, qPCR for lincRNAs and PTGs detection in Roche LightCyler 480 system (Roche, Mannheinm, Germany) was performed using SYBR Green (CWBIO, Beijing, China, CW0957). The qPCR data were analyzed using the 2-∆∆CT method and R scripts were used to perform related linear regression analysis.

3. Results

3.1. Transcripts Assembly and lincRNAs Identification

In order to identify and analyze lincRNAs associated with pig muscle growth and fat deposition, RNA-seq data involving three tissues of white Duroc and Erhualian F2 full-sib pig individuals were obtained from a previously published study and the main process of identifying lincRNAs was shown in Figure 1(a).

After removing the unqualified reads in raw reads with the “trimmomatic” command, approximately 217.2 of 225.5 million clean reads were obtained and mapped to the pig reference genome (Sus scrofa 11.1) by HISAT2 (Table 1). Then, we assembled the transcriptome for each sample (liver, abdominal fat, and longissimus dorsi muscle) by StringTie and synthesized all the transcripts into nonredundant transcripts by using StringTie-Merge program. After merging the nonredundant transcripts, about 15.38% (14,392 of 93,571) of the transcripts were defined as intergenic transcripts. The 581 putative lincRNAs encoded by 500 gene loci of 14,392 intergenic transcripts were obtained using the selection conditions shown in Figure 1(a), and the annotation files of all lincRNAs were shown in File S1. As all the samples were from female pigs, these 500 gene loci were not distributed on the Y chromosome and the detailed distribution map was shown in Figure 1(c). However, 93 of these 581 lincRNAs were novel lincRNAs (Figure 1(b)).

3.2. Characterization of Identified lincRNAs

There are many differences between protein-coding transcripts and lincRNAs, such as the length of transcripts, exon number, and exon length. To verify these differences, we performed the comparison analysis. In this process, the putative lincRNAs (including novel lincRNAs and known lincRNAs) were compared with the protein-coding transcripts, and the characteristics of novel lincRNAs were analyzed during the comparison process. We obtained 45,788 protein-encoding transcripts, including 22,342 annotated genes in the Ensembl pig sequence database by reconstructing transcripts. In addition, the pig lincRNAs annotation file contains 12,103 known lincRNAs transcripts that corresponded to 7,381 lincRNAs genes. After identification, we found that the average transcripts lengths of protein-coding, novel lincRNAs, and known lincRNAs were 3285bp, 669bp, and 953bp, respectively. And the average transcripts length of novel lincRNAs was shorter than the other two transcripts (Figure 2(a)). In addition, the average exon lengths of the protein-coding transcripts, novel lincRNAs, and known lincRNAs were 283, 286, and 376, respectively. And it shows that the exon length of novel lincRNAs is shorter than that of known lincRNAs but longer than that of protein-coding transcripts (Figure 2(b)). Meanwhile, the average exon number of the protein-coding transcripts, novel lincRNAs, and known lincRNAs are 11.6, 2.3, and 2.5, suggesting that the average number of exons is approximately the same as that of known lincRNAs but much smaller than that of protein-coding transcripts (Figure 2(c)). The comparison results (shorter transcript length, smaller exon number, and longer exon length) of novel lincRNAs with protein-coding genes are basically consistent with previous reports [37, 38].

3.3. Expression Analysis of lincRNAs

Previous reports showed that lincRNAs had lower expression level when compared with protein-coding transcripts [39]. To test this conclusion, we compared the average expression levels of the 581 putative lincRNAs (488 known lincRNAs and 93 novel lincRNAs) with the protein-coding transcripts to explore the lincRNAs expression profile. We found that putative lincRNAs (novel lincRNAs vs. known lincRNAs, 0.76FPKM vs. 1.34FPKM) showed indeed significantly lower expression level than that of protein-coding transcripts (3.95FPKM) (Figure 3(a)). Then, in order to identify the effects of differentially expressed lincRNAs from different tissues on muscle growth and fat deposition, we compared the three tissues (abdominal fat vs. longissimus dorsi (af vs. ldm), abdominal fat vs liver (af vs. liver), and longissimus dorsi muscle vs. liver (ldm vs. liver)) through the “DEseq2” package in the R software to get the differentially expressed genes (DEGs) between tissues. Finally, we obtained the numbers of the DEGs between different tissues which were 128 (af vs. ldm), 185 (af vs. liver), and 181(ldm vs. liver), respectively (Figures 3(b)3(d)), and the specific expression of these DEGs between different tissues was shown in Tables S1S3. There are more differentially expressed lincRNAs between different tissues of the same individual, and the lincRNAs expressed in the same tissue of different individuals have little difference (Figures 3(b)3(d)). Meanwhile, 271 differentially expressed protein-coding genes were also identified (Figure 3(e)).

3.4. Nearest Neighbor Analysis of lincRNAs

Previous studies indicate that lincRNAs may act in a cis-manner to regulate the expression levels of their neighboring genes [40]. It is valuable to predict the corresponding function of lincRNAs by identifying the function of a protein-coding gene transcribed near lincRNAs (<100 kb). To predict the functional characteristics of putative lincRNAs, first of all, we obtained the adjacent protein-coding genes (<100 kb) of all 272 differentially expressed genes (DEGs) (lincRNAs) (Table S4). And then, we queried each protein-coding gene in the DAVID database and performed related pathway analysis [11] (Table S5). The results of the DAVID analysis showed that 216 of the 272 DEGs were significantly involved in 33 biological processes and 15 pathways. DELs had 147 nearby protein-coding genes involved in biological processes and pathways involved in muscle growth and fat deposition, such as skeletal muscle tissue development, actin cytoskeleton reorganization, lipid metabolic process, beta-alanine metabolism, fatty acid degradation, and glucagon and signaling pathway (Figure 4(a)). Interestingly, there are some links between these major pathways involved in muscle growth and fat deposition and several bridge genes among these pathways such as CPTIA, GCGR, PPP3CA, and PPARA. Some genes in the ADH and AOC family genes are related to the biological processes of muscle growth or fat deposition (Figure 4(b)).

3.5. QTLs Analysis of lincRNAs and Functional Prediction

A total of 3,794 quantitative trait loci (QTLs) were obtained by DELs screening, involving regeneration, production, meat, carcass quality, health, and exterior (Table S6). About 18.1% (686 / 3,794) of QTLs were related with fat deposition, and 5.81% (219/ 3,794) of QTLs related to the muscle growth (Figure 5(a)). The chromosome distribution of linRNAs which were associated with muscle growth and fat deposition QTLs was shown in Figure 5(b). Even if we narrowed down all of the QTLs fragments obtained (<10kb), we still obtained some QTLs associated with muscle growth and fat deposition. After that, we continued to screen out the lincRNAs corresponding to QTLs related to muscle growth and fat deposition within 10k of the QTLs fragments and obtained a total of 7 lincRNAs related to muscle growth and fat deposition. Then, the functions of these 7 lincRNAs were predicted by analyzing the functions of corresponding adjacent genes (<100k), and the correlationships between these lincRNAs and their adjacent genes were shown in Figure 5(c). To obtain more accurate functions of lincRNAs, the range of adjacent genes of these 7 lincRNAs was then narrowed down to 10k, and the functions of these adjacent genes were statistically analyzed. These adjacent genes mainly included FGGY, GCGR, SLC25A10, PPP3CA, MRPL12, and DPYD. MRPL12 are associated with muscle growth and GCGR and SLC25A10 with fat deposition, while PPP3CA, DPYD, and FGGY are related not only to muscle growth but also to fat deposition (Figure 5(d)). Further analysis of these six target genes revealed that PPP3CA and GCGR existed in the glucagon signaling pathway and DPYD were involved in the beta-alanine metabolism pathway. This result further underscored the potential mechanisms by which differentially expressed lincRNAs are associated with muscle growth and fat deposition [36, 41].

3.6. Correlation Validation between lincRNAs and PTGs

Based on the expression analysis of “htseq-count,” we randomly selected 4 pairs of lincRNAs and PTGs to verify the relationship between putative lincRNAs and the corresponding PTGs by performing the qPCR experiments in corresponding tissues that included longissimus dorsi muscle, liver, and abdominal fat. All of the selected lincRNAs and their PTGs were differentially expressed, and their primer sequences for PCR experiments are shown in Table S7. First, we calculated r0 and p0 between each pair of lincRNAs and its PTGs by the cor.test() function in the R software. r0 and p0, respectively, represented the Pearson correlation coefficient and p-value of each pair of lincRNAs and its PTGs obtained in RNA-seq. Then, through the qPCR experiment, we obtained the relevant r and p, which represented the Pearson correlation coefficient and p-value of each pair of lincRNAs and its PTG calculated from the experimental data, and the biological significance of r and p is similar to that of r0 and p0, respectively. Finally, through comparison, we found that the results of 4 pairs of quantitative verifications were consistent with the trend of RNA-seq sequencing (Figures 7(a)7(d)). 2 of the 4 lincRNAs (MSTRG.7054 vs. GCGR, MSTRG.29576 vs. DPYD) were randomly selected from the lincRNAs that corresponded to QTLs related to muscle growth or fat deposition within 10k of the QTL fragment. The correlations were 0.811 and 0.825, respectively, and the P-values were all less than 0.01.

4. Discussion

Numerous genomic studies have shown that the majority of the long noncoding RNAs in mammalian genomes are lincRNAs [17, 42]. Due to the similarity between pigs and humans in physiology process, organ development, disease research, and so on, they are widely used as important animal models, but the lincRNAs identified in pigs are far less perfect than that in humans and mice [13]. Many types of lincRNAs in pigs are still unidentified and their function remains unknown [37]. In particular, the mechanisms of the biological processes of lincRNAs underlying the regulations of the muscle growth and fat deposition still require further research. In our study, we comprehensively identified and investigated lincRNAs associated with muscle growth and fat deposition in pigs based on RNA-seq data published on NCBI.

In this study, we presented the systematical transcriptome profiling of the raw reads obtained from six samples from three tissues, including liver, abdominal fat, and longissimus dorsi muscle, using high throughput RNA-seq technology [43]. After the “fastqc” command operation, the number of clean reads obtained from the six samples was 32831240-40055406, and about 77.4% of the reads uniquely located in the reference genome, which not only further improved the gene annotation in the pig genome but also made a certain contribution to increase the number of genomes. It has been reported in the literature that most lincRNAs exhibit high tissue specificity [44, 45]. Therefore, it can be preliminarily concluded that some of these 581 putative lincRNAs are not only associated with pig muscle growth and fat deposition, but also specifically expressed. Due to people's increasing attention to pork quality, the two biological processes of muscle growth and fat deposition have gradually become the focus of later studies, and the related lincRNAs will also become the focus of research.

We characterized the identified 581 putative lincRNAs and found that they were consistent with previously reported results. Compared to protein-coding transcripts, there was shorter transcript length, smaller exon number, and longer exon length. Based on the FPKM expression levels of putative lincRNAs and protein-coding transcripts, we analyzed the specificity of liver, abdominal fat, and longissimus dorsi muscle, and then it was also demonstrated that the lincRNAs in the three tissues were more tissue-specific than the protein-coding transcripts.

To further predict whether putative lincRNAs are involved in muscle growth and fat deposition-related biological processes, we selected all DELs adjacent genes (<100 kb) for gene ontology and pathway analysis. During the analysis, we are indeed enriched with pathways and biological processes associated with muscle growth and fat deposition (Figure 4(a)). In terms of fat deposition, we have enriched the two KEGG pathways of fatty acid degradation and glucagon signaling pathway, as well as the biological processes of lipid metabolic process. Both in fatty acid degradation and glucagon signaling pathway, the CPT1A gene has been screened. In some related literature reports, the carnitine palmitoyltransferase 1 (CPT1) family gene members are necessary to transport the long-chain fatty acids into the mitochondria for oxidation. The CPT1A gene can induce mouse obesity through diet, and its expression positively correlated with BMI (R = 0.46) [4649]. Related literature reports that PPARA that was screened together in the lipid metabolism process and the glucagon signaling pathway was a known lipid metabolism regulator [50]. In addition, PPP3CA and GCGR are also enriched in the glucagon signaling pathway. It is reported that PPP3CA is associated with intramuscular fat [51], and activation of GCGR results in decreased plasma triglyceride and cholesterol levels [52]. In terms of muscle growth, it is worth mentioning that not only do both PPP3CA and GCGR exist in the glucagon signaling pathway, but PPP3CA also exists in the skeletal muscle tissue development, and GCGR exists in the actin cytoskeleton reorganization. In addition, we have enriched two remarkable biological processes, including regulation of skeletal muscle tissue development and phospholipase C-activating G-protein coupled receptor signaling pathway and a KEGG pathway of beta-alanine metabolism. Similarly, the HOXD13 gene was simultaneously screened in both biological processes. By studying the enriched skeletal muscle tissue development, it was found that some of the target genes enriched in the pathway were significantly associated with muscle growth traits, such as TAZ, and MYL (MYL3 and MYL6) family genes. It has been reported that mutations in the TAZ gene will cause muscle weakness in mice [53]. MYL3 influence the differentiation and growth performance in quail muscle [54] and participate in the contraction of the chicken's muscles [55]. In pigs, muscle growth and muscle fiber differences are also regulated by MYL3 [56], and MYL6 was associated with fibroblast development in mice [57] (Figure 4(b)). In the phospholipase C-activating G-protein coupled receptor signaling pathway, the ESR and PTH1R were screened. It has been reported that knockout, knockdown, and overexpression of ESR will lead to overall muscle hypertrophy [58], and PTH1R accelerates myocyte differentiation [59]. In the beta-alanine metabolism, the selected DPYD gene is related not only to muscle growth but also to fat deposition. It is reported that DPYD has an effect on the longissimus muscle of Angus cattle [60], and it is related to the increase of beef marbled fat [61]. All target genes enriched in the corresponding pathways listed above are adjacent genes of DELs, and they are obviously involved in the two biological processes of muscle growth or fat deposition. So, it is speculated that some putative lincRNAs may participate in the biological processes of muscle growth and fat deposition by regulating their adjacent protein-coding genes. However, the mechanism by which individual lincRNAs regulate their adjacent genes is worthy of experimental verification and further study.

After performing the analysis of the adjacent gene pathway of DELs, and in order to further verify the function of the putative lincRNAs, we made further functional prediction of DELs through QTLs analysis. When the range of QTLs was reduced to 10kb, we screened out 7 linRNAs corresponding to QTLs that related to muscle growth and fat deposition. Then, we screened out all the 6 adjacent genes (10kb) of 7 lincRNAs. The obtained adjacent genes were analyzed and most of them were related to muscle growth or fat deposition, which was basically consistent with our expected results. Among all 6 the adjacent genes obtained, MRPL12 is related to muscle growth. Relevant literature reports that reduced MRPL12 levels will result in overall mitochondrial translation defects in fibroblasts [62]. The genes GCGR and SLC25A10 are associated with fat deposition. It is reported that GCGR gene plays an important role in lipid metabolism [63]. At the same time, other literature reports that glucagon can bind to its receptor GCGR to regulate glucose levels and fatty acid oxidation in patients with type 2 diabetes by regulating the level of cAMPs and pka-independent pathways (Figure 6(a)) [64]. Related literature reports that the biological process by which citrate transports in exchange for malate across the mitochondrial membrane is the start of fatty acid synthesis in adipose tissue or the liver. Moreover, experiments have shown that SLC25A10 plays an important role in providing malate for citric acid transport required for fatty acid synthesis, and SiRNA was used for knockdown of SLC25A10 which is effective in reducing lipid accumulation in adipose tissue (Figure 6(b)) [65, 66]. The genes related to muscle growth and fat deposition are PPP3CA, DPYD, and FGGY. Some literature reports that calcineurin (PPP3CA) can promote the dephosphorylation of the nuclear factor of activated T cells (NFAT), allowing it to smoothly enter the nucleus from the cytoplasm, and the target gene of NFAT can participate in skeletal muscle differentiation and Hypertrophic Gene Program (Figure 6(c)) [67, 68]. The DPYD gene is involved in β-alanine metabolism and participated in the modification of lysine, and lysine plays an important role in lipid peroxidation and muscle and muscle bond development (Figure 6(d)). There are reports in the literature that DPYD has an effect not only on the longissimus muscle of cattle but also on the fat content of beef marble [60, 61]. FGGY gene is associated with amyotrophic lateral sclerosis [69], and FGGY mRNA is stronger in mouse WATs than in brown adipose tissue and is enhanced in gonadal fat by diet-induced obesity [70].

Both muscle and fat cells are derived from mesenchymal stem cells and are closely related during the development of the organism [71, 72]. At present, it has been found that some functional genes have a regulatory effect on both muscle and fat growth, such as PPP3CA, DPYD, and FGGY in this study. Therefore, it is possible that some related lincRNAs also regulate muscle and fat growth at the same time, which provides a good reference for the later studies of lincRNAs in muscle and fat.

Combined with pathway analysis and QTLs results, it was found that there were three of the six genes screened during QTL analysis enriched in pathways related to muscle growth or fat deposition. Both PPP3CA and GCGR are enriched in the glucagon signaling pathway, and DPYD is enriched in the beta-alanine metabolism pathway. It is also worth noting that the 7 lincRNAs and the 6 target genes in the process of QTLs analysis are almost one-to-one correspondence (Figure 5(d)). For example, only FGGY and PPP3CA are in the 10k range of MSTRG.35876.1 and MSTRG.40844.1, respectively. The different splice isoforms of MSTRG.7054 (MSTRG.7054.1, MSTRG.7054.3, and MSTRG.7054.4) all fall near the GCGR; meanwhile, the different splice isoforms of MSTRG.29576 (MSTRG.29576.3 and MSTRG.29576.2) all fall near the DPYD. Thus, we can more strongly speculate that the putative lincRNAs may participate in two biological processes of muscle growth or fat deposition by regulating its PTGs.

In this experiment, we identified and characterized all the putative lincRNAs in the liver, abdominal fat, and longissimus dorsi of pigs. After a series of analyses, especially pathway analysis and QTLs analysis, and according to quantitative verification, it was found that putative lincRNAs had a good correlation with the PTGs. Thus, we further inferred that many putative lincRNAs might be involved in muscle growth and fat deposition-related processes. Moreover, our study provides new insights into the discovery and annotation of lincRNAs related to pig muscle growth and fat deposition, especially in the analysis of differential expression of lincRNAs. Several genes are related to muscle growth and fat deposition, such as FGGY, GCGR, and PPP3CA, which are ideal candidates for later experimental verification.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

No potential conflicts of interest were disclosed.

Authors’ Contributions

Changchun Li conceived and designed the experiments and explained the data. Guoting Chen analyzed main content of the data with the help of Cheng Zou, Chengchi Fang, Gaoli Shi, and Lin Chen. Xiaofang Cheng performed the experiment with the help of Mengxun Li and Jingxuan Li. Guoting Chen wrote the paper with the help of Changchun Li and Chengchi Fang.

Acknowledgments

The work was supported by National Natural Science Foundation of China (NSFC, 31872322), the Fundamental Research Funds for the Central Universities (2662017PY030), National High Technology Research and Development Plan (863, 2011AA100302), National Natural Science Foundation of China (NSFC, 31601859), and China Postdoctoral Science Foundation Funded Project (2016M592344).

Supplementary Materials

Supplementary Materials: File S1: All lincRNAs comment files. Table S1: The specific expression of the DEGs between two tissues (af vs. ldm). Table S2: The specific expression of the DEGs between two tissues (af vs. liver). Table S3: The specific expression of the DEGs between two tissues (ldm vs. liver). Table S4: The basic information of the 272 DEGs. Table S5: Results of functional enrichment in gene ontology and pathway analysis. Table S6: Results of relevant QTLs screened by DELs. Table S7: PCR primers of 4 pairs of lincRNAs and PTGs. (Supplementary Materials)