Identification of Candidate Genes and Regulatory Competitive Endogenous RNA (ceRNA) Networks Underlying Intramuscular Fat Content in Yorkshire Pigs with Extreme Fat Deposition Phenotypes

Intramuscular fat (IMF) content is vital for pork quality, serving an important role in economic performance in pig industry. Non-coding RNAs, with mRNAs, are involved in IMF deposition; however, their functions and regulatory mechanisms in porcine IMF remain elusive. This study assessed the whole transcriptome expression profiles of the Longissimus dorsi muscle of pigs with high (H) and low (L) IMF content to identify genes implicated in porcine IMF adipogenesis and their regulatory functions. Hundreds of differentially expressed RNAs were found to be involved in fatty acid metabolic processes, lipid metabolism, and fat cell differentiation. Furthermore, combing co-differential expression analyses, we constructed competing endogenous RNAs (ceRNA) regulatory networks, showing crosstalk among 30 lncRNAs and 61 mRNAs through 20 miRNAs, five circRNAs and 11 mRNAs through four miRNAs, and potential IMF deposition-related ceRNA subnetworks. Functional lncRNAs and circRNAs (such as MSTRG.12440.1, ENSSSCT00000066779, novel_circ_011355, novel_circ_011355) were found to act as ceRNAs of important lipid metabolism-related mRNAs (LEP, IP6K1, FFAR4, CEBPA, etc.) by sponging functional miRNAs (such as ssc-miR-196a, ssc-miR-200b, ssc-miR10391, miR486-y). These findings provide potential regulators and molecular regulatory networks that can be utilized for research on IMF traits in pigs, which would aid in marker-assisted selection to improve pork quality.


Introduction
Pork is the most widely consumed meat worldwide, accounting for 40% of human red meat consumption, and its quality directly affects human health [1]. Intramuscular fat (IMF) content is an important characteristic of pork quality and is positively associated with its palatability, marbling score, tenderness, juiciness, and overall meat quality [2]. It is generally accepted that a higher IMF content tends to be an indicator of high-grade pork [3]. Therefore, breeding pigs with more IMF can produce more palatable pork. With relevance to meat quality in pig and human health, IMF content is now an important economic trait in pig breeding programs [4,5].
In pigs, IMF has high heritability, with estimated values varying from 0.21 to 0.86 and an approximate average of 0.5 [6,7]. This high heritability makes porcine IMF a suitable and important indicator for selection breeding programs focused on obtaining high-quality pork. However, the direct selection of IMF content by traditional breeding methods during pig breeding is extremely difficult to implement, partly because of phenotyping complexity and partly because it is influenced by different factors, such as species, breed, age, sex, and nutritional levels [8]. Molecular breeding methods could be a valid method for swine IMF improvement; therefore, identifying potential molecular markers for IMF is now an important task in genetic research and pig molecular breeding programs [4][5][6][7]. Moreover, pigs are an ideal model for human obesity-related research because they share many similarities with humans at the physiological and genomic levels [9]. Therefore, studies on the molecular mechanisms of IMF deposition are important for improving the economic efficiency of pigs and human health.
As a polygenic trait, IMF is a complex metabolic process determined by the hyperplasia and hypertrophy of adipocytes and involves many biological processes and pathways, which are regulated by various molecules, such as hormones, transcription factors, proteincoding genes, and non-coding RNAs [10][11][12]. However, when compared with other traits, dissecting the genetic basis of IMF content is still limited, owing to its complexity. Previous findings have shown that porcine intramuscular adipogenesis is regulated not only by protein-coding RNAs (mRNAs) but also by non-coding RNAs (ncRNAs). For example, miR-146a-5p targets SMAD4 and TRAF6 to inhibit porcine intramuscular preadipocytes adipogenesis through the TGF-β and AKT/mTORC1 signaling pathways [13]. Further, LMFlnc1 (lncRNA) promotes adipogenesis in porcine IMF cells by upregulating CAV-1 expression via miR-199a-5p [14]. Moreover, Sus_circPPARA promotes differentiation and hinders proliferation in porcine intramuscular preadipocytes [11]. Despite the increasing number of transcriptomic studies aimed at dissecting IMF deposition, cascade events related to this trait remain mostly unknown. Furthermore, non-coding RNAs related to IMF in pigs are largely unknown, and the mechanisms by which they regulate IMF deposition are unclear and require further study. Therefore, to fully elucidate the complex gene networks and molecular cascade pathways that regulate IMF deposition in pigs, it is crucial to expand basic knowledge by further characterizing coding and non-coding transcriptomes.
Transcriptome studies comparing individuals with extreme phenotypes of a trait are useful in identifying gene pathways and networks with divergent expression among groups of livestock species. IMF content can vary from 2 to 10% in different swine breeds [15,16]. For example, it is typically higher in indigenous Chinese pig breeds than in Western pig breeds and commercial pigs [17]. Pig breeds in China have high contents of subcutaneous and intramuscular fat and good meat quality [18]. However, commercial pig breeds, such as Duroc, Landrace, and Yorkshire, are selected for better lean muscle growth, leading to a decrease in the deposition of fat, including IMF, resulting in reduced meat palatability [12]. Studies have also shown that IMF varies considerably among pure commercial breeds [7,11,16]. For example, Yorkshire pigs are a typical lean-type Western breed and have a low level of fat deposition as a whole. However, phenotypic variation in IMF content is still observed in this population [16]. Previous transcriptome studies using pigs of different breeds have provided relevant results related to the genetic expression patterns and networks underlying IMF production [17,19]; however, few studies have been conducted on individuals of the same breed with distinct IMF contents to identify consistent candidate genes [15,16]. Improving IMF production from a genetic perspective in lean pigs is one of the major goals of pig breeding programs to improve meat quality, and it would be particularly interesting to analyze changes in the transcriptome and regulatory factors of finishing lean pigs with divergent IMF content.
In this study, we collected tissues of the longissimus dorsi muscle of castrated finishing Yorkshire males of similar ages and body weights, but with divergent back fat thickness and IMF. Whole-transcriptome sequencing of this tissue was performed to explore the gene profiles and identify candidate lncRNAs, circRNAs, miRNAs, and mRNAs associated with IMF deposition by evaluating differences in expression and identifying the pathways in which the differentially expressed genes are involved. Furthermore, potential lncRNA/circRNA-miRNA-mRNA co-regulatory networks were constructed to elucidate the complex genetic architecture related to swine IMF. This research could provide a com-prehensive and bioinformatic resource to study the regulatory mechanisms of pig IMF deposition mediated by non-coding RNAs. These results might also contribute to improvements in pork quality at the genetic and molecular levels and provide insight into human obesity and related diseases.

Characterization of the Longissimus Dorsi Transcription
Two cDNA libraries, an rRNA-depleted library and an miRNA library, were constructed. After the redundant and low-quality reads were removed, the rRNA-depleted library contained 82,402,852, 74,367,466, 124,581,402, 90,505,406, 12,780,608, and 137,892,896 clean reads with greater than 93.05% Q30 scores from the H1-3 and L1-3 samples, respectively, which were used to identify the mRNAs, lncRNAs, and circRNAs (Table S2). Among them, 94.94%, 95.94%, 95.69%, 95.86%, 93.86%, and 95.76% of reads from H1-3 and L1-3, respectively, were mapped to the pig reference genome (S. scrofa 11.1; Table S2). In addition, for the small RNASeq libraries, there were 32,998,980 raw reads (H: 9,685,832-12,734,725 for each library) and 34,161,122 raw reads (L: 9,295,100-13,349,696 for each library) in the H and L groups, respectively. After trimming and filtering, 32,168,194 (H) and 33,271,882 (L) clean reads were obtained. In total, 2,127,279 (H) and 2,143,330 (L) known miRNA reads and 23,953 (H) and 23,846 (L) novel miRNA reads were obtained after a series of analyses (Table S3). In total 21,353 mRNAs, 1088 miRNAs, 9980 lncRNAs, and 13,521 circRNAs were obtained from the H and L groups. Further, 341 mRNAs, 36 miRNAs, 91 lncRNAs, and 178 circRNAs were identified as differentially expressed (DE) RNA molecules between the two comparison groups (Table 1 and Figure 1A-D). In this study, we collected tissues of the longissimus dorsi muscle of castrated finishing Yorkshire males of similar ages and body weights, but with divergent back fat thickness and IMF. Whole-transcriptome sequencing of this tissue was performed to explore the gene profiles and identify candidate lncRNAs, circRNAs, miRNAs, and mRNAs associated with IMF deposition by evaluating differences in expression and identifying the pathways in which the differentially expressed genes are involved. Furthermore, potential lncRNA/circRNA-miRNA-mRNA co-regulatory networks were constructed to elucidate the complex genetic architecture related to swine IMF. This research could provide a comprehensive and bioinformatic resource to study the regulatory mechanisms of pig IMF deposition mediated by non-coding RNAs. These results might also contribute to improvements in pork quality at the genetic and molecular levels and provide insight into human obesity and related diseases.

Transcriptional Profiling of mRNAs
In the H vs. L group set analyses, 341 DE mRNAs (DEMs) were detected, with expression levels of 134 mRNAs downregulated and those of 207 mRNAs upregulated in the H group (Table 1 and Supplementary File S1), including a series of DE genes directly related to the regulation of IMF deposition, such as ADIPOQ (log2FC = −1.325), CEBPA (log2FC = −1.173), ALOX12B (log2FC = −3.624), LEP (log2FC = −3.148), ACACA (log2FC = −1.226), DGAT2 (log2FC = −1.608), ACLY (log2FC = −6.554), IP6K1 (log2FC = −7.907), PLIN1 (log2FC = −1.759), CYB5R1 (log2FC = −1.508), ACBD7 (log2FC = 2.336), SPP1 (log2FC = 2.864), ADAMTS8 (log2FC = 2.921). GO analysis revealed that these DEMs were significantly enriched (p < 0.05) in terms related to fatty acid biosynthesis and lipid metabolism, such as the lipid metabolic process (GO:0006629), brown fat cell differentiation (GO:0050873), white fat cell differentiation (GO:0050872), regulation of cholesterol transport (GO:0032374), fatty acid metabolic process (GO:0006631) ( Figure 2A, Table 2, and Supplementary File S1). KEGG analysis revealed that these DEMs were significantly enriched (p < 0.05) in steroid biosynthesis, steroid hormone biosynthesis, glycerolipid metabolism, adipocytokine signaling pathway, and steroid hormone biosynthesis ( Figure 3A and Supplementary File S1). We then performed a comprehensive bioinformatic analysis of PPI networks among DEMs to further extract relevant information from the identified transcriptome data. As shown in Figure 4, the PPI network of DEMs comprised 32 nodes, 62 edges, and four significantly enriched pathways/terms related to lipid metabolism (the regulation of the triglyceride biosynthetic process, fatty acid metabolic process, lipid metabolic process, and adipocytokine signaling pathway). From this integrated analysis of DEM, GO, and KEGG pathway results, we focused on DEMs that interacted with three or more other genes and were associated with one or more GO/KEGG terms. LEP, ACACA, ACLY, PLIN1, HSD17B7, and CPT1A were identified as hub genes in the network. Based on these results, we assumed that all of these could be promising candidate genes that affect porcine longissimus dorsi fatty acid and lipid metabolism and ultimately the accumulation of IMF.   We then performed a comprehensive bioinformatic analysis of PPI networks among DEMs to further extract relevant information from the identified transcriptome data. As shown in Figure 4, the PPI network of DEMs comprised 32 nodes, 62 edges, and four significantly enriched pathways/terms related to lipid metabolism (the regulation of the triglyceride biosynthetic process, fatty acid metabolic process, lipid metabolic process, and adipocytokine signaling pathway). From this integrated analysis of DEM, GO, and KEGG pathway results, we focused on DEMs that interacted with three or more other genes and were associated with one or more GO/KEGG terms. LEP, ACACA, ACLY, PLIN1, HSD17B7, and CPT1A were identified as hub genes in the network. Based on these results, we assumed that all of these could be promising candidate genes that affect porcine longissimus dorsi fatty acid and lipid metabolism and ultimately the accumulation of IMF.   Circle nodes, genes/proteins; rectang nodes, KEGG pathway or GO terms. Pathways or GO terms are colored with gradient color fro yellow to blue, with yellow for smaller p-value and blue for bigger p-value. According to trend an ysis, genes/proteins are colored in red (representing up-regulation) and green (representing dow regulation). Interactions are shown as solid lines between genes/proteins, and edges of KEGG pa ways/Go terms are presented as dashed lines.

Expression Patterns of lncRNAs
In total, 9980 lncRNAs (including 9367 annotated lncRNAs and 613 novel lncRNA were identified on all chromosomes. Five types of lncRNAs were detected, as follows: least one splice junction shared with a reference transcript (j), a transfrag falling entire within a reference intron (i), generic exonic overlap with a reference transcript (o), u known intergenic transcript (u), and exonic overlap with a reference on the oppos strand (x). Most (49.92%) lncRNAs were type j, and the minority (9.79%) were type i (Fi ure 5A). Most (76.92%) of lncRNAs were more than 1200 bp in length, whereas the leng of mRNAs was evenly distributed, ranging from 0 to > 3600 bp ( Figure 5B). Compared protein-coding genes, lncRNAs mostly contained two or three exons ( Figure 5C). In add tion, the expression levels of lncRNAs were lower than those of protein-coding genes (Fi ure 5D).
In addition, 91 DE LncRNAs (DELs) were detected in the H vs. L group compariso with expression levels of 44 lncRNAs upregulated and those of 47 downregulated in t H group (Table 1 and Supplementary File S2). To determine the functions of the identifi DELs in IMF deposition, three independent algorithms (antisense, mRNA sequence com plementarity; cis, genomic location; and trans, expression correlation) were used to pr dict the target genes of all DELs. GO and KEGG enrichment analyses were then perform for the target genes. GO analysis revealed that many GO terms related to triglyceride m abolic processes and lipid metabolism were significantly enriched (p < 0.05), including t cholesterol metabolic process (GO:0008203), positive regulation of triglyceride metabo process (GO:0090208), and triglyceride catabolic process (GO:0019433) ( Figure 2B, Tab 3, and Supplementary File S2). KEGG analysis revealed that these DELs were significant Circle nodes, genes/proteins; rectangle nodes, KEGG pathway or GO terms. Pathways or GO terms are colored with gradient color from yellow to blue, with yellow for smaller p-value and blue for bigger p-value. According to trend analysis, genes/proteins are colored in red (representing up-regulation) and green (representing down-regulation). Interactions are shown as solid lines between genes/proteins, and edges of KEGG pathways/Go terms are presented as dashed lines.

Expression Patterns of lncRNAs
In total, 9980 lncRNAs (including 9367 annotated lncRNAs and 613 novel lncRNAs) were identified on all chromosomes. Five types of lncRNAs were detected, as follows: at least one splice junction shared with a reference transcript (j), a transfrag falling entirely within a reference intron (i), generic exonic overlap with a reference transcript (o), unknown intergenic transcript (u), and exonic overlap with a reference on the opposite strand (x). Most (49.92%) lncRNAs were type j, and the minority (9.79%) were type i ( Figure 5A). Most (76.92%) of lncRNAs were more than 1200 bp in length, whereas the length of mRNAs was evenly distributed, ranging from 0 to > 3600 bp ( Figure 5B). Compared to proteincoding genes, lncRNAs mostly contained two or three exons ( Figure 5C). In addition, the expression levels of lncRNAs were lower than those of protein-coding genes ( Figure 5D). digestion and absorption, steroid hormone biosynthesis, sphingolipid metabolism, and cholesterol metabolism ( Figure 3B and Supplementary File S3).    (Table 1 and Supplementary File S2). To determine the functions of the identified DELs in IMF deposition, three independent algorithms (antisense, mRNA sequence complementarity; cis, genomic location; and trans, expression correlation) were used to predict the target genes of all DELs. GO and KEGG enrichment analyses were then performed for the target genes. GO analysis revealed that many GO terms related to triglyceride metabolic processes and lipid metabolism were significantly enriched (p < 0.05), including the cholesterol metabolic process (GO:0008203), positive regulation of triglyceride metabolic process (GO:0090208), and triglyceride catabolic process (GO:0019433) ( Figure 2B, Table 3, and Supplementary File S2). KEGG analysis revealed that these DELs were significantly enriched (p < 0.05) in pathways involved in fat and cholesterol metabolism, such as fat digestion and absorption, steroid hormone biosynthesis, sphingolipid metabolism, and cholesterol metabolism ( Figure 3B and Supplementary File S3).

Detection of miRNA Expression
In total, 36 DE miRNAs (DEMiRs) were defined between the H and L groups, with expression levels of 10 miRNAs downregulated and those of 26 upregulated in the H In total, 178 DE CircRNAs (DECs) were identified between the H and L groups, with expression levels of 92 circRNAs upregulated and those of 86 downregulated in the H group (Table 1 and Supplementary File S4). To understand the functional distribution of DECs, functional annotation of DEC parent genes was performed by GO and KEGG pathway analyses. GO analysis revealed that these target genes were significantly enriched (p < 0.05) in terms related to lipid metabolism, such as lipid kinase activity (GO:0001727), CoA carboxylase activity (GO:0016421), regulation of lipid kinase activity (GO:0043550), phosphatidylinositol biosynthetic process (GO:0006661), and phosphatidylinositol-translocating ATPase activity (GO:0004012) ( Figure 2C, Table 4, and Supplementary File S4). KEGG analysis revealed that these DEC parent genes were significantly enriched (p < 0.05) in pathways involved in lipid metabolism, such as purine metabolism and galactose metabolism ( Figure 3C and Supplementary File S5).

qRT-PCR Validation
To validate potential interactions in the networks, the expression levels of ssc-miR-196a and ssc-miR7138-5p, related ceRNAs, were measured by qRT-PCR. The expression of ssc-miR-196a was upregulated in the H group. In contrast, that of its target genes, including novel_circ_011355, ENSSSCT00000066779, ZMYND19, PARM1, and ADAMTS8 was downregulated (Figure 8). The expression level of ssc-miR7138-5p was downregulated in the H group. In contrast, that of its target genes, including ENSSSCT00000066779, ENSSSCT00000076340, ENSSSCT00000068476, ETV4, NPR3, and STARD3, was upregulated ( Figure 8). In addition, validation of the RNA-seq results was carried out using qRT-PCR for another randomly selected 15 genes, including two DEMs (CHRNA3 and

qRT-PCR Validation
To validate potential interactions in the networks, the expression levels of ssc-miR-196a and ssc-miR7138-5p, related ceRNAs, were measured by qRT-PCR. The expression of ssc-miR-196a was upregulated in the H group. In contrast, that of its target genes, including novel_circ_011355, ENSSSCT00000066779, ZMYND19, PARM1, and ADAMTS8 was downregulated (Figure 8). The expression level of ssc-miR7138-5p was downregu-lated in the H group. In contrast, that of its target genes, including ENSSSCT00000066779, ENSSSCT00000076340, ENSSSCT00000068476, ETV4, NPR3, and STARD3, was upregulated ( Figure 8). In addition, validation of the RNA-seq results was carried out using qRT-PCR for another randomly selected 15 genes, including two DEMs (CHRNA3 and TMEM38B), five DECs (novel_circ_002804, novel_circ_008940, novel_circ_001557, novel_circ_011588, novel_circ_003997), four DEMiRs (ssc-miR-200b, miR-1983-z, miR-486-y, miR-10-x), and four DELs (ENSSSCT00000080712, ENSSSCT00000090965, MSTRG.12825.1, ENSSSCT00000070023). The expression patterns of these transcripts were highly consistent with those obtained by RNA-seq (Table S4), indicating the high reproducibility and reliability of the gene expression profiles obtained in this study. were highly consistent with those obtained by RNA-seq (Table S4), indicating the high reproducibility and reliability of the gene expression profiles obtained in this study. Figure 8. The expression patterns of ssc-miR-196a, ssc-miR-7138-5p, related target genes were verified by qRT-PCR. The data represented the Mean ± SD from 5 biological replicates, and each measurement was repeated 3 times. ** indicates p < 0.01, * indicates p < 0.05.

Discussion
During the last decade, meat producers have started to focus more on pork quality. The IMF content, which represents the amount of fat, including phospholipids, triglycerides, and cholesterol within muscles, is an important factor that is positively associated with overall meat quality [20,21]. Pork with a higher IMF tends to be of better quality, resulting in higher overall acceptability [22]. Given the importance of IMF in the economics of pig meat production, clarifying the molecular mechanisms underlying IMF deposition in pigs is of great significance. In the present study, longissimus dorsi whole-transcriptomes of two groups of Yorkshire purebred finishing castrated boars divergent in IMF content were comprehensively analyzed to compare differences in the expression profiles of mRNAs and non-coding RNAs. This was performed to identify candidate genes and pathways related to the divergent IMF deposition and to interpret the complex molecular cascade events related to the variability observed in this trait. In addition, potential non-coding RNA-miRNA-mRNA co-expression ceRNA networks were constructed and preliminarily validated, providing new insights into the molecular mechanisms of porcine IMF deposition.
Comparing gene expression between individuals with divergent traits in the same population can reduce noise owing to different genetic backgrounds [12]. In the present study, all pigs were selected based on the same breed, sex, and weight to allow candidate genes to be more dependable, thereby excluding unreliable factors identified owing to different breeds or different feed conditions. Yorkshire pigs are a typical lean-type western breed that has been intensively selected over the past few decades to increase lean meat production and reduce fat deposition; however, they still exhibit considerable phenotypic variation in fatty traits in the population [16]. Too much intrapopulation variation is a disadvantage for commercial production [4]. Moreover, Yorkshire breeds often serve

Discussion
During the last decade, meat producers have started to focus more on pork quality. The IMF content, which represents the amount of fat, including phospholipids, triglycerides, and cholesterol within muscles, is an important factor that is positively associated with overall meat quality [20,21]. Pork with a higher IMF tends to be of better quality, resulting in higher overall acceptability [22]. Given the importance of IMF in the economics of pig meat production, clarifying the molecular mechanisms underlying IMF deposition in pigs is of great significance. In the present study, longissimus dorsi whole-transcriptomes of two groups of Yorkshire purebred finishing castrated boars divergent in IMF content were comprehensively analyzed to compare differences in the expression profiles of mRNAs and non-coding RNAs. This was performed to identify candidate genes and pathways related to the divergent IMF deposition and to interpret the complex molecular cascade events related to the variability observed in this trait. In addition, potential non-coding RNA-miRNA-mRNA co-expression ceRNA networks were constructed and preliminarily validated, providing new insights into the molecular mechanisms of porcine IMF deposition.
Comparing gene expression between individuals with divergent traits in the same population can reduce noise owing to different genetic backgrounds [12]. In the present study, all pigs were selected based on the same breed, sex, and weight to allow candidate genes to be more dependable, thereby excluding unreliable factors identified owing to different breeds or different feed conditions. Yorkshire pigs are a typical lean-type western breed that has been intensively selected over the past few decades to increase lean meat production and reduce fat deposition; however, they still exhibit considerable phenotypic variation in fatty traits in the population [16]. Too much intrapopulation variation is a disadvantage for commercial production [4]. Moreover, Yorkshire breeds often serve as the first dam line in crossbred (Duroc × (Landrace × Yorkshire)) pig production [23], the first sire line in crossbred (Duroc × (Yorkshire × Landrace)) pig production [24], and the sire line in the binary hybridization production of Chinese local pigs [25], which means that they contribute approximately one-third or one-half of the genetic material passed on to the offspring. Therefore, it will be beneficial and efficient to genetically improve the IMF of Yorkshire pigs. Because of the particular characteristics of the Yorkshire pig, it is particularly interesting to analyze changes in the transcriptome and regulatory networks of finishing Yorkshire pigs with divergent IMF content.
In total, 341 DEMs were identified from the sequencing data, of which the expression levels of 207 DEMs were upregulated in the H group, including a series of lipogenic genes, such as ACACA, ACLY, PTGR1, DGAT2, and CYB5R1. These differentially expressed genes are involved in various aspects of de novo fatty acid synthesis, including direct and indirect regulation, and play catalytic roles in fatty acid biosynthesis [18,[26][27][28]. Moreover, the expression differences in ACACA, ACLY, PTGR1, DGAT2, and CYB5R1 between the H and L groups were consistent with the increasing IMF deposition trend (Supplementary File S1), indicating that these genes might act as positive regulators of the IMF deposition process. Furthermore, many DEMs (LEP, PIP5K1B, SREBF2, CPT1A, ALOX12B, IP6K1, NR1H2, ACBD7, DGKI, SPP1, CYP1A1, LPIN1, HSD17B7, FFAR4, ADIPOQ, ADAMTS8, and HSD11B1) were significantly enriched in fatty acid biosynthesis and lipid metabolism-related GO terms and KEGG pathways. Among these genes, LEP encodes a protein, leptin, primarily secreted by white adipocytes but also expressed in other tissues including skeletal muscle, was also reported to be associated with porcine longissimus dorsi IMF content [29]. Individuals with a high body fat composition have higher levels of leptin [30]. In our trial, LEP expression tended to be upregulated in the H group, which agrees with the typical fat composition of this group when compared to that of the L group. The ADIPOQ gene encodes a protein hormone, adiponectin, which is secreted by adipocytes and is involved in the regulation and inhibition of lipogenesis and the stimulation of fatty acid oxidation [31]. Consistent with our results, the expression of ADIPOQ was higher in Lantang, a high-IMF pig breed, than in Landrace, a low-IMF pig breed [32]. Meanwhile, HSD17B7 encodes an enzyme involved in cholesterol biosynthesis [33], whereas PLIN1 belongs to the periplasmin family of proteins and is associated with the formation of intracellular lipid droplets [34]. Further, CPT1A encodes carnitine palmitoyltransferase-1, an enzyme responsible for transporting long-chain fatty acids for Î 2 -oxidation. In our trial, CPT1A expression was upregulated in pigs with a high IMF content, which is consistent with a previous report showing that its expression level is positively correlated with IMF content [35]. Interestingly, in our PPI network results, we identified six hub genes, including LEP, ACACA, ACLY, PLIN1, HSD17B7, and CPT1A, most of which are involved in fatty acid and lipid metabolism. These results suggest that genes responsible for fatty acid and lipid metabolism in the longissimus dorsi tissue significantly differ between high-and low-IMF Yorkshire pigs.
In this study, 91 DELs were identified from sequencing data, including 60 known lncR-NAs and 31 novel lncRNAs. To further investigate the difference in IMF deposition between H and L groups, we predicted the target genes of these DELs (Supplementary File S8). Some of these differentially expressed lncRNA-target genes have been reported to be important for IMF deposition. LPL (lipoprotein lipase), the target gene of ENSSSCT00000050073, is involved in fatty acid catabolism, and its expression level is positively associated with swine IMF content [36]. The PPARG gene for MSTRG.14355.1, encoding the core regulator of the PPAR signaling pathway, regulates lipid metabolism and glucose homeostasis, promotes adipocyte differentiation and fat deposition, and harbors polymorphisms that could significantly affect IMF deposition in pigs [37]. HMGCR, the target gene of ENSSSCT00000072909, and MSTRG.4737.2, a cholesterol-synthesis-limiting enzyme, harbor one polymorphism that shows a positive relationship with IMF in pigs [38].  : SREBF1), are also closely related to fat metabolism [29,39]. KEGG and GO analyses revealed that many fat and cholesterol metabolism-related pathways and terms were significantly enriched, and some genes related to lipid metabolism were enriched multiple times, including APOA1, APOA4, APOB, MTTP, ABCG5, and ABCG8 [40][41][42][43][44]. Furthermore, via trans prediction, we found that ENSSSCT00000062623, ENSSSCT00000069577, ENSSSCT00000069784, ENSSSCT00000075507, ENSSSCT00000077869, MSTRG.129.2, and MSTRG.5755.1 all formed targeting relationships with APOA1, APOA4, APOB, MTTP, ABCG5, and ABCG8 (Supplementary File S8), and the expression levels of these seven DELs were all higher in the high IMF group, suggesting that these DELs might facilitate lipid transport and metabolism. Based on these results, we suspected that one of the main roles of these lncRNAs is to regulate IMF deposition in pigs, and further highlighting their detailed mechanisms in this process would provide a strong basis for future investigations.
CircRNAs have been found to play important roles in adipogenesis and lipid metabolism [45,46]. In this study, we identified 178 DECs. The functions of circRNAs are assumed to be related to those of their host genes. Some of these DEC host genes have been reported to be important for fat deposition (Supplementary File S4). LRP6, the host gene of novel_circ_000878, is a well-established factor that governs lipid generation and secretion, which regulates body fat and glucose homeostasis by modulating nutrient-sensing pathways and mitochondrial energy expenditure [47]. PLIN1, the host gene of novel_circ_011588, plays an important role in regulating lipolysis and lipid storage in adipocytes and is reported to be a candidate gene affecting porcine IMF content [48]. Furthermore, KEGG and GO analyses revealed that some lipid metabolic process-related pathways and terms were significantly enriched, and some host genes were enriched multiple times, including ALOX15 (novel_circ_006541), PDGFA (novel_circ_005676), DAB2IP (novel_circ_008787), and PCYT1A (novel_ circ_009957). ALOX15, a lipoxygenase isoform, is involved in the metabolism of linoleic and arachidonic acids [49]. DAB2IP (novel_circ_008787), a Ras GAP, regulates lipid droplet homeostasis by acting as a GAP for RAB40C [50]. Further, PCYT1A (novel_ circ_009957) is the major isoform of the key enzyme CTP (choline phosphate cytidylyltransferase), which is essential for phosphatidylcholine synthesis during lipid metabolism [51]. PDGFA (novel_circ_005676) plays a vital role in the proliferation and maintenance of adipocyte progenitors in dermal adipose tissue through the PI3K-Akt pathway [52] and was found to be a potential candidate marker that regulates pig growth and fat deposition [53]. Meanwhile, in the GO analyses, glycerolipid biosynthetic process (GO:0045017) and arachidonate 12-lipoxygenase activity (GO:0004052) were significantly enriched for both DEMs and DECs; further, some GO terms were also enriched for both DEMs and DECs but only significantly enriched for one of them, such as the lipid metabolic process (GO:0006629) and regulation of lipid localization (GO:1905952) (Supplementary File S4). These results indicate that the differences in lipid metabolism and fat deposition between the high-and low-IMF pigs are regulated not only by mRNAs but also by circRNAs.
Studies have shown that mRNAs, lncRNAs, and circRNAs regulate the expression of each other (when sharing the same miRNA-binding sites) by functioning as competing endogenous RNAs (miRNA sponges) during adipocyte differentiation. For example, the lncRNA ADNCR inhibits adipocyte differentiation by functioning as a ceRNA for miR-204, thereby augmenting expression of the miR-204-target gene SIRT1 [67]. Further, circFLT1 and lncCCPG1 sponge miR-93 to regulate the proliferation and differentiation of adipocytes by promoting lncSLC30A9 expression [68]. In this study, combined with the co-differentially expressed DEMs, DELs, DECs, and DEMiRs, we constructed ceRNA regulatory networks, which showed that 30 lncRNAs and 61 mRNAs exhibited crosstalk with each other through 20 miRNAs, and that five circRNAs and 11 mRNAs showed crosstalk through four miRNAs. Consistent with previous studies [19,69], this also indicates that IMF deposition in pigs results from a balance in gene expression. Furthermore, from the ceRNA networks, we observed a series of ceRNA subnetworks that might play key roles in the regulation of IMF deposition. For example, MSTRG.12440.1 and its target FFAR4 exhibited crosstalk through miR-141-y and miR-200-y response elements, ENSSSCT00000085972 and its target LEP showed crosstalk through miR-10391, ENSSSCT00000050073 and its target CEBPA exhibited crosstalk through miR486-y, and novel_circ_011355 and its targets PARM1, ZMYND19, and ADAMTS8 showed crosstalk through miR-196a. FFAR4, LEP, PARM1, CEBPA, ZMYND19, and ADAMTS8 have been reported to be involved in fat metabolism [70][71][72]. In addition, several nodes were found to be shared by both lncRNA-miRNA-mRNA co-regulatory networks and circRNA-miRNA-mRNA co-regulatory networks, namely ssc-miR-196a, ZMYND19, PARM1, and ATAMTS8. From these data, it could be inferred that the identified non-coding RNAs participate in the intramuscular adipogenesis process by acting as ceR-NAs and that the genes involved in ceRNA regulatory networks might play an important role in IMF deposition through molecular synergism and the upregulation of important pathways, which should be studied in the future.

The Experimental Animals and Sample Collection
A Yorkshire finishing pig resource population was housed at Anhui Lvjian Breeding Pig Co., Ltd. (Quanjiao, China) under consistent and standard environmental conditions. In total, 75 healthy castrated males of similar ages and weights (approximately 170-daysold and 125 kg, respectively) were selected. The live backfat thickness was measured 5 cm from the left dorsal midline between the last third and fourth ribs using real-time Bmode ultrasonography. All 75 live backfat thickness phenotypic values were fitted normal distribution. The two tails of the distribution, including 5 samples for each, were established by two investigation groups. According to these values, five pigs with extremely high (16.66 ± 0.78 mm) and five pigs with extremely low (9.14 ± 0.74 mm) live backfat thickness were slaughtered in the same batch. After slaughter, the average three-point backfat thicknesses of carcasses at the shoulder end of the dorsal midline, the thoracolumbar junction, and the lumbar spine junction were measured using vernier calipers as the carcass backfat thickness. Samples of longissimus dorsi muscle at the 3rd/4th last rib were collected [73], with a portion of the samples stored at −20 • C for IMF determination and the rest frozen in liquid nitrogen until RNA isolation. The IMF was measured as a percentage using the Soxhlet extraction method [74]. The average IMF values were 7.18% (SD = ±0.013) and 1.64% (SD = ±0.006) in the high and low live backfat thickness groups, respectively. The carcass backfat thickness values were 38.20 mm (SD = ±1.29) and 18.06 mm (SD = ±1.42) for these groups, respectively. The significance of the differences in these three fatness traits between the two groups was assessed using a t-test in SPSS 20.0, with all p values < 0.001. Subsequently, three individuals were randomly selected from the extremely high and low live backfat thickness, IMF, and carcass backfat thickness groups to comprise the extremely high IMF (H) and low IMF (L) groups for RNA-seq analysis.

mRNA/lncRNA/circRNA Sequencing and Data Analysis
Total RNA was isolated from each longissimus dorsi sample (six individuals) using the TRIzol reagent kit (Invitrogen, Carlsbad, CA, USA). After the integrity, purity, and quality of the isolated RNA were tested, samples with an RNA integrity number greater than seven were used for further analysis. Ribosomal RNA (rRNA) was removed from the DNA-free RNA using the Ribo-ZeroTM kit (Epicenter, Madison, WI, USA). DNA and rRNA-free RNA were used to create a library using the NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, E7530L). The library was sequenced using an Illumina Novaseq6000 instrument (Gene Denovo Biotechnology Co., Guangzhou, China). Raw data from the six pig transcriptomes were uploaded to the NCBI SRA under the accession number PRJNA821451.
Clean data were obtained by filtering the adapters, unknown bases, and low-quality reads. HISAT2 (2.10) [75] was used to map the clean reads to the Sus scrofa 11.1 reference genome (http://ftp.ensembl.org/pub/release-103/fasta/sus_scrofa/ (accessed on 2 May 2021)). Gene abundance was quantified with RSEM [76]. The fragments per kilobase of transcript per million (FPKM) value was then used to represent the expression levels of the genes.

Identification of circRNAs
After the clean reads were aligned to the porcine reference genome, junctions of the unmapped reads were identified using a back-splice algorithm. Findcirc software (v1.0) [81] was used to predict circRNAs. The expression levels of circRNAs were reflected by the number of mapped back-splicing junction reads per million mapped reads.

Small RNA Library Construction and Sequencing
Total RNA from each longissimus dorsi sample (six individuals) was isolated. RNA quality and quantity were determined using a Bioanalyzer 2100 (Agilent Technologies, Palo Alto, CA, USA). The adapters were then added, and the 36-48 nt RNAs were enriched by polyacrylamide gel electrophoresis. The final products were amplified using reverse transcription (RT)-PCR to construct a cDNA library. They were then sequenced using an Illumina Novaseq6000 by Gene Denovo Biotechnology Co. (Guangzhou, China). Raw miRNA sequencing data were uploaded to the NCBI Biotechnology Information database (PRJNA824228).

Alignment and Identification of Small RNA
To obtain clean tags, raw reads were further filtered, as with conventional processing [82]. All clean tags were aligned with small RNAs in the GenBank database (Release 209.0) and the Rfam database (Release 11.0) to identify and remove rRNA, scRNA, snoRNA, snRNA, and tRNA. All clean tags were also aligned to the reference genome to remove tags mapped to exons, introns, or repeat sequences. Next, the filtered tags were searched against the miRBase database (Release 22) to identify the known porcine miRNAs. The unannotated tags were predicted to be and identified as novel miRNAs using mirdeep2 software, according to the tag positions in the genome and their hairpin structures.

Prediction of Target Genes of miRNAs
Miranda (v3.3a) and TargetScan (v7.0) were used to predict miRNA targets. The intersection of the results was selected as the predicted miRNA-target genes.

Differentially Expressed RNA Identification and Pathway Analysis
DESeq2 software [83] was used to detect differentially expressed mRNAs (DEMs), circRNAs (DECs), miRNAs (DEMiRs), and lncRNAs (DELs), with values of p < 0.05 and |log2fold-change (FC)| ≥ 1. Gene Ontology (GO) term and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed for DELs (using the DEL target genes, which were predicted in a cis, trans, and antisense manner), DEMiRs (using the DEMiR target genes), DECs (using the DEC parent genes), and DEMs using the DAVID tool, and p < 0.05 was considered significant.

Protein-Protein Interaction (PPI) Analysis of Differentially Expressed mRNAs
PPI analysis of DEMs was performed using the STRING database (http://string-db. org/cgi/input.pl (accessed on 30 July 2022)), which is a well-known tool used to predict PPIs. The following network model was generated based on information gained from up to four levels of functional analysis: fold-change of genes or proteins, PPIs, KEGG pathway enrichment, and GO enrichment. DEMs were mapped to STRING with a confidence score > 0.4, which is a reasonable score used to limit the number of interactions to those with higher confidence that are much more likely to be true positives. PPI networks were generated using Cytoscape software (v3.3.0) [84].

Construction of the lncRNA/circRNA-miRNA-mRNA Network
To reveal the functions of and interactions among ncRNAs and mRNAs, we constructed an ncRNA-mRNA regulatory network. The ceRNA network was constructed based on ceRNA theory. The prediction of target genes for differentially expressed miRNAs was the first step. The mRNA-miRNA, lncRNA-miRNA, or circRNA-miRNA expression correlations were evaluated using the Spearman rank correlation coefficient (SCC). Pairs with an SCC value < −0.7 were selected as negatively co-expressed lncRNA-miRNA, mRNA-miRNA, or circRNA-miRNA pairs. mRNA, lncRNA, and circRNA were miRNAtarget genes, and all RNAs were differentially expressed. The lncRNA-mRNA and circRNA-mRNA expression correlations were evaluated using the Pearson correlation coefficient (PCC). Pairs with a PCC value > 0.9 were selected as co-expressed lncRNA-mRNA pairs or circRNA-mRNA pairs. Both the mRNA and lncRNA in such pairs were targeted and negatively co-expressed with a common miRNA or both mRNA and circRNA in such pairs were targeted and negatively co-expressed with a common miRNA. The ceRNA regulatory network was constructed by assembling all co-expressed competing triplets, which were identified previously herein and visualized using Cytoscape software (v3.3.0) [84].
The primers used are listed in Supplemental Table S1. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH; for mRNA, lncRNA, and circRNA) and U6 (for miRNA) were used as endogenous controls. The synthesized cDNA was used as a template for RT-qPCR using a CFX96 TouchTM Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). All reactions were performed in triplicate for each sample. The 2 −∆∆CT method was used to quantify changes in relative gene expression. Significant differences were analyzed by a t-test using SPSS 22.0 with the following criteria: p < 0.05 (*), p < 0.01 (**).

Conclusions
The present study provides a comprehensive understanding of the differences in the whole-transcriptome profiles of longissimus dorsi tissues between high-and low-IMF Yorkshire finishing pigs. In summary, 341 DEMs, 91 DELs, 178 DECs, and 36 DEMiRs were characterized and found to be widely involved in terms and pathways related to fatty acid metabolic processes, lipid metabolism, and fat cell differentiation, indicating that non-coding RNAs are abundant in the longissimus dorsi muscle and are involved in porcine adipocyte adipogenesis and lipid metabolism. Furthermore, we constructed ceRNA regulatory networks and found a series of ceRNA subnetworks that might play a key role in the regulation of IMF deposition. Specifically, the results showed that functional lncRNAs and circRNAs (such as MSTRG.12440.1, ENSSSCT00000066779, ENSSSCT00000076340, ENSSSCT00000050073, MSTRG.12825.1, ENSSSCT00000085972, ENSSSCT00000080916, and novel_circ_011355, novel_circ_011355) act as ceRNAs of important fat deposition-related mRNAs (LEP, ZMYND19, CEBPA, ADAMTS8, ZMYND19, IDH3, CYB5R1, DGKI, CRTC3, IP6K1, FFAR4, and STARD3, etc.) by sponging functional miRNAs (such as ssc-miR-196a, ssc-miR-200b, ssc-miR10391, miR486-y, ssc-miR7138-5p, and novel-m0064-3p, miR-2683-z). These results provide a strong foundation for future investigations, as well as potential regulators and molecular regulatory networks for in vitro and in vivo experiments to study the detailed mechanisms underlying IMF deposition in pigs. This work also has an important significance for the pig industry, especially for breeding pigs based on the IMF trait.