Comparative Transcriptomic Analysis of mRNAs, miRNAs and lncRNAs in the Longissimus dorsi Muscles between Fat-Type and Lean-Type Pigs

In pigs, meat quality and production are two important traits affecting the pig industry and human health. Compared to lean-type pigs, fat-type pigs contain higher intramuscular fat (IMF) contents, better taste and nutritional value. To uncover genetic factors controlling differences related to IMF in pig muscle, we performed RNA-seq analysis on the transcriptomes of the Longissimus dorsi (LD) muscle of Laiwu pigs (LW, fat-type pigs) and commercial Duroc × Landrace × Yorkshire pigs (DLY, lean-type pigs) at 150 d to compare the expression profiles of mRNA, miRNA and lncRNA. A total of 225 mRNAs, 12 miRNAs and 57 lncRNAs were found to be differentially expressed at the criteria of |log2(foldchange)| > 1 and q < 0.05. The mRNA expression of LDHB was significantly higher in the LD muscle of LW compared to DLY pigs with log2(foldchange) being 9.66. Using protein interaction prediction method, we identified more interactions of estrogen-related receptor alpha (ESRRA) associated with upregulated mRNAs, whereas versican (VCAN) and proenkephalin (PENK) were associated with downregulated mRNAs in LW pigs. Integrated analysis on differentially expressed (DE) mRNAs and miRNAs in the LD muscle between LW and DLY pigs revealed two network modules: between five upregulated mRNA genes (GALNT15, FKBP5, PPARGC1A, LOC110258214 and LOC110258215) and six downregulated miRNA genes (ssc-let-7a, ssc-miR190-3p, ssc-miR356-5p, ssc-miR573-5p, ssc-miR204-5p and ssc-miR-10383), and between three downregulated DE mRNA genes (IFRD1, LOC110258600 and LOC102158401) and six upregulated DE miRNA genes (ssc-miR1379-3p, ssc-miR1379-5p, ssc-miR397-5p, ssc-miR1358-5p, ssc-miR299-5p and ssc-miR1156-5p) in LW pigs. Based on the mRNA and ncRNA binding site targeting database, we constructed a regulatory network with miRNA as the center and mRNA and lncRNA as the target genes, including GALNT15/ssc-let-7a/LOC100523888, IFRD1/ssc-miR1379-5p/CD99, etc., forming a ceRNA network in the LD muscles that are differentially expressed between LW and DLY pigs. Collectively, these data may provide resources for further investigation of molecular mechanisms underlying differences in meat traits between lean- and fat-type pigs.


Introduction
Pork is a vital source of protein, energy and iron for humans, and occupies an important position in the world meat consumption [1]. Meat quality is attracting more and more attention in recent years. Intramuscular fat (IMF) content is a primary pork quality indicator, affecting meat quality traits including tenderness, juiciness, flavor and nutritional value [2]. IMF refers to the fat extracted from muscle tissue, which is deposited in muscle fiber, outer membrane and inner membrane and is considered to be a late developing fat storage during the process of fat deposition [3]. Meat quality is also influenced by the composition of different types of skeletal muscle fibers including type I (slow/oxidative),

Animals and Tissue Samples
A total of three castrated male Laiwu pigs with unrelated genetic background reared under identical feeding conditions and environment were randomly selected from the Conservation Center of Laiwu pigs (Laiwu, Shandong, China) at 150 d. Similarly, three castrated male DLY pigs of 150 d of age were randomly selected from Beihe Pig Breeding Co., Ltd. (Laiwu, Shandong, China). The animals were allowed access to food and water ad libitum and were housed under identical conditions. To minimize animal suffering, the pigs were weighed and electrically stunned before death. Immediately after slaughter, about 30 g sample of the Longissimus dorsi muscle (LD) at the third lumbar vertebra of each pig was collected into a 2 mL cryogenic vial (Corning, NY, USA) and frozen in liquid nitrogen for further study [15,23]. For RNA-seq, samples were transported in dry ice to BGI (Shenzhen, China) and stored at −80 • C for total RNA extraction. The same samples were also used for real-time quantitative PCR (RT-qPCR). All animal care and experimental procedures were in accordance with the guidelines for the care and use of laboratory animals prescribed by Shandong Agricultural University and the Ministry of Agriculture of China (No. SDAUA-2021-097).

IMF Content Measurement
The IMF content in the LD muscle tissues of pigs was measured by Soxhlet petroleum ether extraction method [24]. Briefly, after removing the fascia, blood vessels and connective tissue, skeletal muscle sample was ground using a meat grinder. Ten grams of the sample was dried to a constant weight in a drying oven (Yiheng Scientific Instrument Co., Ltd., Shanghai, China) at 103 • C. Dried meat sample with filter paper package was treated with an appropriate amount of petroleum ether into the Soxhlet bottle and extracted continuously for 6-8 h. Then, the sample was placed in a desiccator overnight and placed in a vacuum drying oven at 80 • C for 12 h and finally weighed. The IMF was calculated by dividing (the sample of dried meat wrapped in the filter paper (g) minus the weight of the dried meat of the filter paper after extraction (g)) by the sample of dried meat (g) multiplied by 100.

Total RNA Isolation, Library Preparation, and Sequencing
Total RNA was extracted from the LD muscle tissues using Trizol (Invitrogen, Carlsbad, CA, USA) according to manufacturer's instructions. The quantity and purity of the total RNA were evaluated using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Total RNA was divided into two samples after preparation, and each was used for mRNA and lncRNA cDNA library and miRNA library construction, respectively. For mRNA and lncRNA cDNA library, the first step involves the removal of ribosomal RNA (rRNA) using target-specific oligos and RNase H reagents. The cleaved RNA fragments were copied into first strand cDNA using reverse transcriptase (MGIEasy RNA, Shenzhen, China) and random primers, followed by second strand cDNA synthesis using DNA Polymerase I (MGIEasy RNA, Shenzhen, China) and RNase H (MGIEasy RNA, Shenzhen, China). The products were enriched with PCR to create the final cDNA library. For miRNA cDNA library, total RNA was purified to obtain 18-30 nt small RNA, which was ligated to a 5 -adaptor and a 3 -adaptor. The adaptor-ligated small RNAs were subsequently transcribed into cDNA by SuperScript II Reverse Transcriptase (Invitrogen, Carlsbad, CA, USA) and then PCR amplified to enrich the cDNA fragments. The library quality was checked by analyzing the distribution of the fragments size using the Agilent 2100 bioanalyzer, and quantified by RT-qPCR. The qualified libraries were sequenced pair end on the BGISEQ-500 platform (BGI, Shenzhen, China).
The statistically significant differentially expressed (DE) mRNAs, miRNAs and lncRNAs were obtained by a q value threshold of <0.05 and |log 2 (fold change)| >1 using the DEseq2 software (v1.36, Simon Anders, Heidelberg, Germany) [30]. Analyses of the differentially expressed mRNAs by GO were performed using Gene Ontology database (http://www.geneontology. org/, accessed on 1 May 2021) [31,32], and KEGG pathway analysis was performed by KEGG database [33]. Pathway enrichment statistical analyses were performed using R package phyper (https://stat.ethz.ch/R-manual/R-devel/library/stats/html/Hypergeometric.html, accessed on 1 May 2021) to calculate p value [34]. Multiple testing correction of p value was applied using R package qvalue (https://bioconductor.org/packages/release/bioc/html/qvalue.html, accessed on 1 May 2021) with corrected q values < 0.05 by Bonferroni [35]. The plots of KEGG enrichment were prepared using R package ggplot2 (https://www.rdocumentation. org/packages/ggplot2/versions/3.3.6, accessed on 1 May 2021) [36]. The heatmap was drawn by pheatmap R package (https://cran.r-project.org/web/packages/pheatmap/index.html, accessed on 1 May 2021) according to the gene expression level in different samples. A proteinprotein interaction network was constructed by the Search Tool for the Retrieval of Interacting Genes (STRING) online software against the Sus scrofa database. Additionally, Qiagen's Ingenuity Pathway Analysis (IPA) was used to identify the predicted upstream regulators of DEGs.

Immunohistochemistry
The LD muscle tissue sections from LW and DLY pigs were processed according to standard protocols following the procedures as described by Zhang et al. [38]. Mouse monoclonal antibody to S46 (myosin heavy chain, slow developmental; 1:200, AB_528376, DSHB, Iowa city, IA, USA) and mouse monoclonal antibody to F59 (myosin heavy chain, all fast isoforms; 1:200, AB_528373, DSHB, Iowa city, IA, USA) were used to incubate the above prepared sections. After overnight storage at 4 • C, the sections were incubated with goat anti-mouse IgG (PV-9002, ZSGB-Bio Co., Ltd., Beijing, China) for 1 h in a 37 • C thermostat box. The sections were washed with PBS (3 × 5 min) after these incubations. All samples were incubated in diaminobenzidine tetrachloride (DAB kit, PN3122, G-CLONE Co., Ltd., Beijing, China) for 1~3 min and counterstained with hematoxylin and color developed in tap water. The sections were then dehydrated, sealed in clear resin, mounted, and observed microscopically for the distribution of positive cells using a bright field of view.

Overexpression and Knockdown Assay
The entire coding region of murine GALNT15 gene was amplified, and polymerase Gflex (TaKaRa, Dalian, China) was used to ensure high fidelity (Table S1). The polymerase chain reaction (PCR) fragments were generated by double enzyme (Hind III and Xho I) digestion and ligated with pcDNA3.1(+) expression vectors (Invitrogen, Carlsbad, CA, USA) by T4 DNA ligase (Thermo Fisher Scientific, Waltham, MA, USA), which were transformed into DH5α (TaKaRa, Dalian, China) competent cells. After being confirmed by bidirectional sequencing and purified using an EndoFree Plasmid Midi Kit (Aidlab, Beijing, China), these plasmids were used for transfecting 3T3-L1 cells. The obtained plasmid was named pcDNA3.1(+)-GALNT15. Empty pcDNA3.1(+) vector was used as the control.
For knockdown (KD) assay, siRNA was designed according to murine GALNT15 mRNA (Shanghai GenePharma Co., Ltd., Shanghai, China). The negative control of the siRNA had the same composition with the siRNA sequence but had no homology with GALNT15 mRNA (Table S2). Three pairs of siRNAs were designed. The most effective siRNA was used to analyze the knockdown effect of siRNA on GALNT15 gene.

RT-qPCR
The total RNA used for above transcriptome analysis was also used for RT-qPCR. The cDNA was synthesized using a Primescript RT Mix Kit with gDNA Clean (Accurate Biotechnology, Changsha, China), and the resultant cDNA was stored at −20 • C for mRNA and lncRNA expression analysis. Synthesis of the cDNA of miRNA was performed using a miRNA 1st strand cDNA synthesis kit (Accurate Biotechnology, Changsha, China) according to the manufacturer's protocol. Specific primers of DE mRNAs and DE ncRNAs were designed based on the gene sequences of pigs reported in GenBank using DNAMAN 8.0 and synthesized by Shanghai Bioengineering Co., Ltd., Shanghai, China.
The qPCR for mRNAs and lncRNAs was conducted by mixing 10 µL of 2×SYBR Green Pro Taq HS Premix (Accurate Biotechnology, Changsha, China), 0.4 µL of each primer (forward and reverse, Table S3), 100 ng of cDNA template and RNase free water up to 20 µL and run on a Light Cycler 96 real-time PCR system (Roche, Basel, Switzerland) with the following program: 95 • C for 30 s, followed by 40 cycles of denaturation at 95 • C for 5 s and annealing and extension at gene-specific temperature (Table S3) for 30 s. The qPCR for miRNA was conducted by mixing 10 µL of 2X SYBR Green Pro Taq HS Premix II (Accurate Biotechnology, Hunan, China), 0.4 µL of miRNA specific primer (Table S3), 0.4 µL miRNA qPCR 3 primer (Accurate Biotechnology, Changsha, China), 100 ng of cDNA template and RNase free water up to 20 µL and run on a Light Cycler 96 real-time PCR system (Roche, Basel, Switzerland) with the following program: 95 • C for 30 s, followed by 40 cycles of denaturation at 95 • C for 5 s and annealing and extension at annealing temperature (Table S3) for 30 s. The melting curves were obtained, and quantitative analysis of the data was performed using the 2 −∆∆CT relative quantification method [39]. The GAPDH and U6 were used as internal reference gene for the relative quantification of mRNA and lncRNA, and miRNA, respectively.

Statistical Analysis
Statistical analysis was performed using Mann-Whitney U test to compare the live weight and IMF content of LW and DLY pigs (n = 3). Statistical analysis was performed using Student's t test to compare the mRNA expression level of the genes, each experiment was repeated four times, and at least three independent experiments were performed. The quantitative data were presented as the mean ± SEM and the p value below 0.05 is considered as significant. Correlation analysis and plot between RT-qPCR and RNA-seq was performed using GraphPad Prism 8.0. p value below 0.05 is considered as significant.

Overview of the Transcriptome Sequencing Data of Porcine LD Muscles
At 150 d, the live weight of DLY pigs was significantly higher than that of LW pigs (p < 0.05, Figure 1A), while the IMF content in the LD muscle of DLY pigs was significantly lower than that of LW pigs (p < 0.05, Figure 1B). Oil red O staining showed that the Longissimus dorsi muscle of LW pigs contained higher fat deposition than that of DLY pigs ( Figure 1C) and immunohistochemical analysis showed that there were a large number of F59 (fast/type2/glycolytic) myofibers and less F46 (slow/type1/oxidative) myofibers in the muscle tissue of DLY pigs compared to LW pigs ( Figure 1D). Three LD muscle tissue samples collected from each of these LW and DLY pigs were used for transcriptome sequencing to identify genes underlying IMF variation and meat quality traits. For each sample, an average of 114.66 million paired-end clean reads and 11.46 Gb of clean data were produced for mRNA and lncRNA cDNA libraries on DNBseq platform. All six samples had at least 90.68% reads equal to or exceeding Q30 and the clean reads ratio of each sample was greater than 92.93% after removing adapters, low quality reads and data containing N, and the total mapping ratio of each sample was greater than 92.38% (Table 1). For miRNA sequencing, an average of 23.89 million clean tag counts of miRNA cDNA libraries were obtained for each LD sample, with more than 98.5% equal to or exceeding Q20, and the total mapping ratio of each sample was greater than 89.07% (Table 2). Furthermore, a total of 30,681 expressed genes were mapped to pig genome, including 20,727 mRNAs, 6702 lncRNAs and 3252 miRNAs. These sequence data were uploaded to NCBI with the accession number PRJNA815878.

Differentially Expressed mRNAs in the LD Muscle between LW and DLY Pigs
Compared with DLY pigs, a total of 225 differentially expressed (DE) mRNAs (105 upregulated and 120 downregulated) were revealed in the LD muscle of LW pigs ( Figure 2A). According to the q value threshold of <0.05 and log 2 (foldchange) > 3, the upregulated and downregulated DE mRNAs in the LD muscle of LW pigs were shown in Tables 3 and 4, respectively. The greatest fold change in mRNA expression was LDHB (upregulated) ( Table 3) and LOC110258854 (downregulated) ( Table 4), respectively. According to the expression level, all of the DE mRNAs in porcine LD muscle were clustered in two clades of LW and DLY as visualized by the heatmap plot ( Figure 2B). The upregulated and downregulated DE mRNAs between LW and DLY pigs were further analyzed by GO and KEGG pathway approaches. For the upregulated DE mRNAs expressed in the LD muscle of LW pigs, GO terms were mainly enriched in chromatin DNA binding, nucleoside binding and aromatase activity, organelle membrane, glycerol biosynthetic process, inosine catabolic process and fatty acid homeostasis ( Figure 2C), and KEGG pathways were related to processes required for metabolic pathways, glucagon signaling and adipocytokine signaling pathways ( Figure 2D). For the downregulated DE mRNAs expressed in the LD muscle of LW pigs, GO terms were mainly enriched in actin-binding and DNA integration ( Figure 2E), and KEGG pathways were related to processes required for metabolic pathways and MAPK signaling pathway ( Figure 2F). Ten of the randomly selected DE mRNAs were validated by RT-qPCR, showing a significant correlation of 0.857 (p < 0.01) ( Figure 2G). Metabolic pathway plays a critical role in skeletal muscle contraction, development and IMF deposition [40]. The DE mRNAs enriched into metabolic pathways include PTGES2, PNP, GALNT15, KYAT1, NNT, LDHB, ST8SIA5, GOT1, PLA2G4E, PLPP1, TKTL2, DHCR24, MTHFD1, AK5, LOC106505238, etc. ( Figure 2H).    A network of protein-protein interactions (PPI) was constructed to reveal the interactions among DE mRNA genes. The highest number of interactions was observed for estrogen related receptor alpha (ESRRA) in upregulated DE mRNAs in the LW pigs. For the downregulated DE mRNAs in LW pigs, the highest number of interactions was observed for versican (VCAN) and proenkephalin (PENK) ( Figure 3A). As is predicted by the ingenuity pathway analysis (IPA), ESRRA and VCAN were regulated by peroxisome proliferative activated receptor, gamma, coactivator 1 alpha (PPARGC1A), whereas FKBP prolyl isomerase 5 (FKBP5) was regulated by vascular endothelial growth factor (VEGF) family and growth hormone ( Figure 3B).

Differentially Expressed miRNAs in the LD Muscle between LW and DLY Pigs
Six upregulated and six downregulated DE miRNAs were identified in the LD muscles of LW pigs compared to DLY pigs ( Figure 4A) and were clustered into two clades of LW and DLY pigs, respectively ( Figure 4B). The greatest fold change in miRNA expression was ssc-miR1379-5p (upregulated) and ssc-miR204-5p (downregulated), respectively, and all of the DE miRNAs in porcine LD muscle are shown in Tables 5 and 6. Predictive analysis using the ncRNA targeting database revealed 341 target genes of DE miRNAs, among which 208 genes were targeted by upregulated miRNAs and 212 genes by downregulated miRNAs, and 79 genes were overlapped. GO and KEGG enrichment analysis on these genes indicated that they were mainly enriched in histone methyltransferase complex ( Figure 4C), metabolic pathways ( Figure 4E) and ubiquitin mediated proteolysis ( Figure 4D), and cAMP, RAP1 and Ras signaling pathways ( Figure 4F). All of the DE miRNAs were validated by RT-qPCR, showing a significant correlation of 0.859 (p < 0.01) ( Figure 4G).

Differentially Expressed lncRNAs in the LD muscle between LW and DLY Pigs
A total of 57 DE lncRNAs, 32 upregulated and 25 downregulated, were identified in the LD muscles of LW pigs compared to DLY pigs ( Figure 5A) and were clustered into two clades of LW and DLY pigs, respectively ( Figure 5B). The greatest fold change in lncRNA expression was LOC110257307 (upregulated) and LOC110259141 (downregulated), respectively; and all of the DE lncRNAs in porcine LD muscle are shown in Tables 7 and 8. Furthermore, 239 genes targeted in cis by DE lncRNAs were predicted, including 139 genes targeted by upregulated DE lncRNAs and 111 genes targeted by downregulated lncRNAs (11 genes were overlapped). A total of 163 genes targeted in trans by DE lncRNAs were predicted, including 40 genes targeted by upregulated DE lncRNAs and 123 genes targeted by downregulated DE lncRNAs.   GO and KEGG analysis revealed that defense responses to virus ( Figure 5C) and metabolic pathways ( Figure 5D) were mainly enriched by cis targeted genes of DE lncR-NAs upregulated in LW pigs. DBird complex ( Figure 5E) and necroptosis and MAPK signaling pathways ( Figure 5F) were mainly enriched by cis targeted genes of DE lncRNAs downregulated in LW pigs. SH3/SH2 adaptor activity ( Figure 5G) and adipocytokine signaling pathway and JAK-STAT signaling pathway ( Figure 5H) were mainly enriched by trans targeted genes of DE lncRNAs upregulated in LW pigs. Structural molecule activity ( Figure 5I) and purine metabolism and endocytosis ( Figure 5J) were mainly enriched by trans targeted genes of DE lncRNAs downregulated in LW pigs. Ten of the randomly selected DE lncRNAs were validated by RT-qPCR, showing a significant correlation of 0.807 (p < 0.01) ( Figure 5K).

Integrated Analysis on DE mRNAs, miRNAs, and lncRNAs in the LD Muscle between LW and DLY Pigs
The mRNA-miRNA regulatory networks were constructed according to gene expression patterns in the LD muscle of LW pigs: "mRNA up-miRNA down" and "mRNA down-miRNA up". Five upregulated mRNA genes, including GALNT15, FKBP5, PPARGC1A, LOC110258214 and LOC110258215, and six downregulated miRNA genes, including ssc-let-7a, ssc-miR190-3p, ssc-miR356-5p, ssc-miR573-5p, ssc-miR204-5p, and ssc-miR-10383, formed one regulatory network ( Figure 6A,B). Three downregulated DE mRNA genes, including IFRD1, LOC110258600 and LOC102158401, and six upregulated DE miRNAs, including ssc-miR1379-3p, ssc-miR1379-5p, ssc-miR397-5p, ssc-miR1358-5p, ssc-miR299-5p and ssc-miR1156-5p, formed another network ( Figure 6A,C). The expression changes of these mRNA ( Figure 6D) and miRNA genes ( Figure 4G) in the LD muscle between LW and DLY pigs were validated by RT-qPCR. Based on the mRNA and ncRNA binding site targeting database, we constructed a regulatory network with miRNA as the center and mRNA and lncRNA as the target genes (Figure 7), forming a ceRNA network in the LD muscles that are differentially expressed between LW and DLY pigs.

The Role of GALNT15 in Lipid Deposition
By prediction, GALNT15, PPARGC1A and FKBP5, are targets of ssc-let-7a ( Figure 6B), among which the role of PPARGC1A and FKBP5 in lipid deposition has beenlinebreak confirmed [41,42]. We further tested the role of GALNT15 in lipid deposition. 3T3-L1 cells transfected with overexpression vector pcDNA3.1(+)-GALNT15 showed increased expression of GALNT15 ( Figure 8A), more cells containing oil droplets and significantly upregulated mRNA expression of adipogenic marker genes CEBPα and FASN ( Figure 8B). Knockdown (KD) of GALNT15 dramatically downregulated the level of oil droplets in 3T3-L1 cells and the mRNA expression of adipogenic marker genes (PPARγ, CEBPα, FASN and SCD) ( Figure 8C), while overexpression of GALNT15 in these cells recovered the oil droplet level. (Figure 8D). These results indicated that GALNT15 plays a positive role in the process of lipid deposition.

Discussion
Analysis of transcriptome profiling of mRNA and noncoding RNA is more and more widely used as a strategy to investigate the mechanism of IMF deposition and muscle development [43]. IMF is the total lipid associated with all cells present in a meat tissue, including extramyocellular lipids and intramyocellular lipids [44]. We found that, at 150d, the IMF of the Longissimus dorsi muscle between LW and DLY pigs was different, and LW pigs had more IMF content and slow myofibers and less fast myofibers in the LD muscle. Studies suggest that the difference in IMF is caused by differences in transcript abundance of genes [45], and there is a significant positive correlation between the content of IMF and the expression of lipid metabolism genes [46]. Due to that IMF deposition greatly differs between fat-and lean-type pig breeds, in this study, to identify other genes especially noncoding RNA genes, we further systematically compared the transcriptomes of mRNA, miRNA and lncRNA in the LD muscle between LW and DLY pigs of 150 d, a time point when intramuscular fat is rapidly deposited in LW pigs [23].
Lean-and fat-type pigs exhibit various differences except IMF, which is caused by genetic background as well as other factors including nutrition and farm management. Omics comparisons between lean-and fat-type pigs were reported in several pig breeds [10][11][12][13][14].
To reveal genes underlying differences in IMF and other meat quality traits between LW and DLY pigs, the individuals of both breeds were sampled from the farms within the same region, and the pigs were reared with the same diet (crude protein 16.5% and digestive energy 13.63 MJ/Kg) and slaughtered at the same date (150 d). To control variations within groups, three individuals with similar live weight and IMF were used for transcriptome analysis. By using these samples with great difference in IMF between LW and DLY pigs, we identified differentially expressed mRNAs, miRNAs and lncRNAs genes that are re-lated to IMF deposition and skeletal muscle development and constructed a regulatory network of these three types of genes, which provides resource data for further elucidating mechanisms underlying meat quality variations between lean-and fat-types of pigs.
Comparison of the mRNA transcriptome changes in the LD muscle between LW and DLY pigs identified 225 DE genes, with LDHB showing the greatest upregulation in LW pigs. LDHB encodes the B subunit of lactate dehydrogenase enzyme, which catalyzes the interconversion of pyruvate and lactate with concomitant interconversion of NADH and NAD+ in a post-glycolysis process [47]. Consistently, the mRNA expression of LDHB is significantly lower in Pietrain pigs than Duroc and Duroc-Pietrain crossbred pigs [48] and LDHB expression is positively correlated with IMF in pigs crossed by (Pietrain × Duroc) boars and (Landrace × Yorkshire) sows [49], suggesting an important role of LDHB in porcine fat deposition traits. KEGG enrichment analysis showed that the upregulated mR-NAs were mainly enriched in metabolic pathways (PTGES2, PNP, GALNT15, KYAT1, NNT, LDHB, ST8SIA5, GOT1, LOC100525112, LOC100737183, LOC100739101, LOC110255237) and adipocytokine signaling pathways (PPARGC1A). GO enrichment analysis showed that the upregulated mRNAs were mainly enriched in glycerol biosynthetic process (GOT1) and fatty acid homeostasis (GOT1). Adipocytokine signaling is a crucial pathway for IMF deposition and lipid metabolism [50]. GOT1 (glutamic-oxaloacetic transaminase 1)-related pathways are glucose metabolism and metabolism. One recent study showed that GOT1 regulates adipocyte differentiation by altering NADPH content [51]. A recent study showed that PTGES2 (Prostaglandin E synthase 2) is essential for effective skeletal muscle stem cell function, augmenting regeneration and strength [52]. NNT (mitochondrial nicotinamide nucleotide hydrogenase) is the major enzyme that generates NADH in mitochondria, and catalyzes the transfer of a hydride between NADH and NADP+ [53]. KEGG enrichment analysis showed that the downregulated mRNA genes were mainly enriched in metabolic pathways (PLA2G4E, PLPP1, TKTL2, DHCR24, MTHFD1, AK5, LOC106505238) and the MAPK signaling pathway (PLA2G4E, FLNC, HSPA2, PFN2, HSP70.2, HSPA6, IGF2). It has been reported that the MAPK pathway is the main regulator of skeletal muscle development [54]. PLA2G4E is one of the important members of the PLA2 family, and it regulates skeletal muscle and metabolic diseases [55]. AK5 (adenylate kinase 5) and AMP signaling are necessary for energy communication between mitochondria, myofibrils and nuclei, as well as metabolic programs that promote cardiac differentiation in stem cells [56]. GO enrichment analysis showed that the downregulated mRNA genes were mainly enriched in actin binding (XIRP1, SLC6A2, FLNC, ENAH, MYH15, DBN1, PFN2, CNN1, ANKRD1, DNASE1), consistent with a report that actin participates in maintaining muscle function and ensuring muscle contraction [57].
To better understand the potential relationship between DE genes, we subsequently performed PPI network analysis and revealed that the most frequently involved gene in the interaction was ESRRA that was upregulated in LW pigs. ESRRA acts as a site-specific transcription factor and interacts with members of the PGC-1 family of transcription cofactors to regulate the expression of most genes involved in cellular energy production as well as mitochondrial biogenesis [58]. By co-activating with ESRRA, PPARGC1A promotes LDHB transcription and regulates skeletal muscle metabolism [59]. PPARGC1A also enhances lipid oxidation to provide energy for sustained muscle contraction and regulates glucose metabolism [60]. In this study, the mRNA expression of PPARGC1A in the LD muscle was higher in LW pigs compared to DLY pigs, suggesting a regulatory role of PPARGC1A in skeletal muscle development. In the gene network constructed by downregulated genes in LW pigs, VCAN and PENK were most frequently involved in the interaction. VCAN was involved in cell adhesion, proliferation, migration and angiogenesis [61,62]. PENK encodes small endogenous opioid peptides and plays a critical role in cell proliferation and differentiation [63]. As an endoplasmic reticulum localized protein, SHISA2 regulates the fusion of muscle satellite cell-derived primary myoblasts [64]. Furthermore, VEGF factors and growth hormone were predicted as upstream regulators of FKBP5, which is closely related to the growth and development of skeletal muscle cells or tissues. In 3T3-L1 cells, overexpression of FKBP5 promotes lipid deposition [42]. In pigs, the expression of FKBP5 is responsive to glucocorticoid receptor NR3C1 [65,66]; it is involved in activating T lymphocyte [67] and significantly contributes to the breeding value for residual feed intake [68]. Their roles in porcine skeletal muscle development require further study.
The regulation of miRNA in a variety of biological processes has attracted more attention. KEGG enrichment analysis showed that the target genes of downregulated miRNA were enriched in metabolic pathways (IMPDH1, GNE, HGSNAT, SAT1, MTMR8, NT5C3A, PANK2, PHYKPL, GALNT11) and the PPAR signaling pathway (ANGPTL4). The PPAR pathway is closely related to glucose homeostasis and lipid metabolism [69]. ANGPTL4 (angiopoietin-like protein family 4) has been shown to regulate lipoprotein metabolism through the inhibition of lipoprotein lipase (LPL) [70]. GO enrichment analysis showed that the target genes of downregulated miRNA included adipose tissue development (SLC25A25). SLC25A25 is reported to regulate lipid metabolism by PGC1-α [71]. KEGG enrichment analysis showed that the target genes of upregulated miRNA were enriched in the cAMP pathway (GNAS, GLI3, AKT1, LOC106505329) and Ras signaling pathways (AKT1, RALA, BCL2L1). The cAMP signaling pathway plays an important role in regulating muscle development and skeletal muscle differentiation [72]. Activation of the PI3K/protein kinase B (Akt) pathway by Ras induces muscle growth but does not alter fiber-type distribution [73]. Studies have shown that AKT1 and AKT2 are indispensable for the regulation of preadipocyte and adipocyte number [74]. GO enrichment analysis showed that the target genes of upregulated miRNA were enriched in histone methyltransferase complex (KDM6A, ZFP64, NCOA6, KMT2C, ZNF335, LOC100626655). Histone methyltransferase complexes are widely involved in the regulation of muscle development [75,76]. Mutations in KDM6A and KMT2D can reduce myocyte differentiation in vitro and damage muscle fiber regeneration in vivo [77].
Integrated analysis of DE mRNA and DE miRNA revealed a targeted relationship among five upregulated DE mRNAs and six downregulated DE miRNAs, of which ssc-let-7a is the most abundant and stable miRNA in porcine muscle development [78,79] and it was predicted to target FKBP5, PPARGC1A and GALNT15 genes. Among these genes, an increase in GLANT15 expression was reported in adipose-derived stromal cells through a differentiation period of 21 d [80]. In this study, we tested the role of GALNT15 in 3T3-L1 cells and found that it played a positive role in lipid deposition, which is consistent with the higher expression of GALNT15 in LW pigs. Among these miRNAs, loss of miR-204 increases insulin secretion and regulates lipid metabolism [81], while downregulated miR-204 promotes skeletal muscle regeneration [82]. However, further studies are needed to verify the role of ssc-miR204-5p in porcine IMF deposition and skeletal muscle development. Three downregulated DE mRNAs and six upregulated DE miRNAs constitute another target relationship, of which IFRD1 was the common target of both ssc-miR1379-5p and ssc-miR1379-3p. It has been reported that IFRD1 can stimulate skeletal muscle regeneration and, as a regulator of MyoD and NF-κB, participate in myoblast differentiation [83], suggesting a possible role in skeletal muscle growth.
LncRNA could participate in a variety of biological processes through different mechanisms [84][85][86]. In this study, we analyzed the interaction between lncRNA and mRNA in cisand trans-interactions. In GO and KEGG analysis, the enrichment of cis-target genes of DE lncRNA mainly includes metabolic pathways (NT5M, OPLAH, ME3, GLUL, ACO2, PEMT, LOC100525869, LOC110259864, PLA2G4E, PLA2G4D, AMY2, LOC100624333), the PI3K-Akt signaling pathway (ITGB6) and the MAPK signaling pathway (PLA2G4E, PLA2G4D, TNFRSF1A). The metabolic pathways mainly include glycerol and phospholipid metabolism. PI3K-Akt and MAPK pathways have been shown to be involved in mediating muscle fiber types and glucose metabolism [87,88]. Some trans-target genes of lncRNA are enriched in the AMPK signaling pathway (IRS1, LEPR) and adipocytokine signaling pathway (IRS1, LEPR), which have been shown to be involved in the regulation of lipid deposition and muscle development [89], suggesting a possible role of these DE lncRNA in regulating porcine meat quality traits. IRS1 has been shown to be necessary for myoblast differentiation and glucose metabolism [90]. LEPR (leptin receptor) mutations have been shown to be associated with lipid metabolism [91]. In this study, we constructed an mRNA-miRNA-lncRNA regulatory network and found potential pathways that may regulate intramuscular fat deposition and muscle development in pigs. Similarly, a comparison of the expression profiling of mRNAs, lncRNAs and circRNAs in the LD muscle of Beijing Black and Yorkshire pigs at 210 d identified DE mRNAs that are mainly enriched in the ECM-receptor interaction, focal adhesion, AMPK signaling pathway, PI3K-Akt signaling pathway, adipocytokine signaling pathway, fatty acid metabolism, and PPAR signaling pathway [92].
Although some genes and pathways identified by this study, such as FKBP5, PI3K-Akt signaling and adipocytokine signaling pathways, were also reported in other pig breeds, still a lot of other genes of mRNA, miRNA and lncRNA expressed in the LD muscle were only identified by comparison between LW and DLY pigs. This is likely caused by specific differences in LW pigs that have super capability of fat deposition and in the age of pigs used for comparison. Whether these genes play similar roles in IMF deposition and other meat quality traits needs further investigations with more individuals of different pig breeds.

Conclusions
In this study, we performed a comparative transcriptome analysis of mRNA, miRNA and lncRNA in the Longissimus dorsi muscle between lean-and fat-type pigs, and found that some upregulated DE mRNAs including LDHB, GALNT15, FKBP5, PPARGC1A, ES-RRA, and their interacting DE miRNAs and lncRNAs were associated with intramuscular fat deposition, and some downregulated DE mRNAs, including IFRD1, VCAN, PENK, LOC110258600, LOC102158401, and their interacting DE miRNAs and lncRNAs were associated with skeletal muscle development. It is necessary to elucidate the role of ESRRA and PPARGC1A in regulating LDHB transcription in porcine skeletal muscle metabolism and intramuscular fat deposition using more samples and pig breeds. The miRNA ssc-let-7a may regulate the expression of GALNT15, PPARGC1A and FKBP5 to affect fat deposition in pigs, which requires further investigations using 3T3-L1 and porcine preadipocytes. Furthermore, it is essential to uncover the regulatory mechanism of these genes and to identify variations in important cis DNA elements to elucidate the differences in meat quality traits between lean-and fat-types of pigs. The results of this study are helpful for the identification of genes underlying IMF variations in pigs.

Data Availability Statement:
The transcriptome data can be accessed from NCBI (https://www. ncbi.nlm.nih.gov/) with the accession number PRJNA815878.

Conflicts of Interest:
The authors declare no conflict of interest.