CircRNAs have been reported to be abundant in animal transcriptomes, and a large number of circRNAs have been identified in humans[22], mice[23], and nematodes[24].Pigs are important farm animals that provide meat for humans and represent an important animal modelfor medical research, and increasing numbers of studies have been conducted on pigs in recent years. Liang et al.[17] identified 149 circRNAs that may be related to muscle growth in the analysis of nine organs and three developmental stages of the Guizhou minipig and found that the host genes of these circRNAs were significantly involved in muscle development and function. These researchers alsoconstructed the first public S. scrofacircRNA database[17]. Morten et al. [25]profiled the expression of circRNA in five brain tissues at up to six time points during fetal porcine development, constituting the first report of circRNA in the brain development of a large animal. An unbiased analysis identified a highly complex regulation pattern of thousands of circRNAs with a distinct spatiotemporal expression profile, suggesting the important function of circRNAs in mammalian brain development. Accurate timing of gene expression is highly important for the developmental stage. Because circRNAs are expressed in a highly spatiotemporally specific manner, it is essential that studies of circRNAs in mammals assess various tissues, conditions, and developmental stages. Ningxiang pigs are known for their delicious meat, but the internal control mechanism governing this property have not been elucidated. To study the spatiotemporal expression pattern of circRNA on lipid synthesis, metabolism and transport in Ningxiang pig tissues, we performed total RNA sequencing across three organs (liver, skeletal muscle, and subcutaneous fat)at four developmental stages (30, 90, 150 and 240 d after birth). Compared with the reference genomes of other pig species, the results of using the Ningxiang pig genome assembled by our research group were more robust (data not published). The mapping rate of each sample exceeded 92%, and the percentage of Q30 bases was more than 92% (Table 1). We identified 61,683 circRNAs in all evaluated tissues. The identified circRNAs had similar distribution trends in chromosome, number, and length among the three tissues. As indicated in a previous study[17], more circRNAs were generated from chromosome 1, chromosome 6 and chromosome 13 than from the other chromosomes, and the largest number was derived from exons. In addition, the number of circRNAs identified at 90d was the largest, and the number of uniquely expressed circRNAs at 90 d was considerably greater than that at other time points. We also identified tissue-specific circular RNAs, which provided important information regarding their functions. For example, we identified 13871, 6274 and 20646 tissue-specific circRNAs in muscle, fat and liver, respectively(Fig. 1C). These findings suggest that specifically expressed circRNAs may play specific roles at specific times during the development of tissues or organs. This possibility is consistent with the findings of other studies, indicating that circRNAs have tissue-/stage-specific expression differences[25–27]. Subsequently, the identification of circRNAs by grouping statistics indicated that 3,362, 1,786 and 3,116 circRNAs with significant differential expression in Ningxiang pig liver, muscle and adipose tissue, respectively.
Since the number of muscle fibers has been determined after birth, skeletal muscle mass increases during postnatal development through a process of hypertrophy, and a similar process may be induced in adult skeletal muscle in response to contractile activity[28]. Each skeletal muscle in an animal contains a different type of fiber. Type I muscle fibers are rich in mitochondria and have a higher oxidative metabolism, while type IIB fibers have fewer mitochondria, and the glycolytic metabolism capacity is higher[29]. However, the role of circRNAs in skeletal muscle development is unknown. The identified circRNAs that were abundantly expressed in muscle include Chr05_81425549_81426478, Chr15_85761324_85761662, Scaffold180_2562774_2570664, Chr02_159842602_159843286, and Chr12_55773288_55775037(Fig. 2E and Additional Table S5). The contractile protein troponin I (TnI), a constituent protein of the troponin complex located on the thin filaments of striated muscle, is involved in calcium-mediated contraction and relaxation. Vertebrate TnI has evolved into three isoforms encoded by three homologous genes, which are expressed under muscle type-specific and developmental regulatory mechanisms[30, 31]. Among these genes, TNNI2 (fast-twitch skeletal muscle isoform, named TNI-fast) is a specific protein of the fast muscle fiber type. During muscle development, Chr02_159842602_159843286, produced by the host gene TNNI2, was abundantly expressed at 30d.It has been reported that the change inMYBPC1 abundance is related to the production of slow-twitch muscle fibers[32, 33]. Interestingly, we observed a peak expression of Chr05_81425549_81426478 transcribed by myosin binding protein C1 (MYBPC1) as a host gene in the developing 90 d pig muscle.Chr15_85761324_85761662, produced by the host gene titin (TTN), was demonstrated to be expressed at 30d and reached its peak at 210 d. It has been reported that TTN is related to intramuscular fat deposition[34], and recent studies have found that the circular RNA TTN acts as a miR-432 sponge to facilitate the proliferation and differentiation of cattle myoblasts via the IGF2/PI3K/AKT signaling pathway[35].In addition, the most abundant circRNAs expressed in muscles are also derived from members of the MYH myosin superfamily, such as MYH2 and MYH6. Recent studies have shown that these genes are associated with skeletal muscle development and influence the quality of pork[36]. These findings suggest that many circRNAs exhibit different expression patterns during muscle development, and changes in the abundance of these circRNAs may affect muscle fiber composition, thereby affecting the properties of meat. Notably, pathway analysis of the genes giving rise to circRNAs peaking at 30d indicated a significant predominance of genes associated with fatty acid biosynthesis, primary bile acid biosynthesis, unsaturated fatty acid biosynthesis, and peroxisome biosynthesis. The impact of circRNAs on the biosynthesis of unsaturated fatty acids will be of great concern and may be related to the higher proportion of unsaturated fatty acids in the diet of piglets. Malic enzyme 1 (ME1) plays an important role in the Krebs cycle for energy metabolism. The mRNA of ME1 was observed to be more abundant in obese Rongchang pigs than in lean Landrace pigs. Furthermore, mRNA abundance changes in ME1 have a marked significant positive correlation with adipocyte volume across the six adipose tissue types[37]. Chr01_176448874_176466186 is produced by the host gene ME1, which is abundantly expressed in adipose tissue (Fig. 3E and Additional TableS8). The expression of Chr01_176448874_176466186 was shown to increase from 90 d, peaking at 150 d and subsequently, with fluctuations, decrease until 210 d. This gene is surmised to be closely related to adipocyte development and adipodeposition. CircRNAs may also be involved in lipid transport; for example, Scaffold155_4811844_4828350 is abundantly expressed in the liver at 90d, and its host gene integer alpha-9 (ITGA9) is known to participate in the sphingolipid pathway(Fig. 4E and Additional TableS11).These lipids are required for active transport, a number of enzymatic processes, membrane formation, and cell signaling[38]. Fructose 2,6-bisphosphate (Fru-2,6-P2) is an important metabolite that controls glycolytic and gluconeogenic pathways in several cell types. The synthesisand degradation of this metabolite are catalyzed by the bifunctional enzyme 6-phosphofructo-2-kinase/fructose 2,6-bisphosphatase (PFKFB1)[39, 40]. Due to the abundant expression of Scaffold180_2562774_2570664 produced by the host gene PFKFB1 in muscle and liver, we surmise that Scaffold180_2562774_2570664 may affect the accumulation of lipids by participating in the regulation of the activity of related metabolic enzymes in gluconeogenesis during muscle and liver development, but this possibility requires follow-up experimental verification. In addition, we found that Chr01_143206410_143210729 was abundantly expressed in all three organs but presented different expression patterns. For example, the expression of Chr01_143206410_143210729 in developmental adipose tissue was shown to increase from 30 d, peak at 150 d and subsequently, with fluctuations, decrease towards 210 d. However, the expression of Chr01_143206410_143210729 in muscle and liver peaked at 90 d and decreased gradually thereafter. Although Chr01_143206410_143210729's host gene SAFB-like transcription modulator(SLTM ) has rarely been reported, its homologous family member SAFB1can interact with and repress the transcriptional activity of various other nuclear receptors, including PPARγ, FXRα, RORα1, PPARα, PPARβ, VDR, SF1, and LRH-1[41]. SAFB1 was also shown to regulate the expression of SREBP‐1c, a bHLH transcription factor that controls lipogenesis in the liver and is induced during over nutrition to facilitate the conversion of glucose to fatty acids and triglycerides for the storage of excess energy[42]. SLTM shares 34% overall identity with SAFB1. SLTM contains SAP and RRM domains, which share higher similarities of 60% and 70%, respectively, with SAFB1[43, 44]. The presence of these shared domains suggests at least some common functions between SAFB1 and STLM; however, it is reasonable to suspect that the circularization of these genes may play important roles in the additive, synergistic, or potentially antagonistic functions exhibited by family members.
Time-series analysis showed that certain differential circRNAs are dynamically expressed during skeletal muscle, liver and adipose tissue development. Although the molecular function of circRNAshas not been fully elucidated, the current research shows that in addition to the function of circRNAs in regulating host gene transcription, protein binding and translation, further research is needed to determine the role of miRNA sponges[45, 46].Competing endogenous RNAs (ceRNAs) are widely present in muscle, fat, and liver. For example, circular RNASNX29 sponges miR-744 to regulate the proliferation and differentiation of myoblasts by activating the Wnt5a/Ca signaling pathway[47]. CircSAMD4A regulates preadipocyte differentiation by acting as a miR-138-5p sponge and thus increasing EZH2 expression[48]. CircRNA_0046366 inhibits hepatocellular steatosis by abolishing the miR-34a-dependent inhibition of PPARα[49]. In our study, STEM analysis and ceRNA analysis were combined to screen out the interaction network pairs closely related to time point changes. For example, G protein signaling 9 (RGS9) is the targeted gene of ssc-miRNA-7137-3p, while the miRNA targeted 8 circRNAs(Chr15_57751128_57770746, Chr14_24730872_24760271, Chr03_19654503_19689919, Chr11_10784978_10824218, Chr06_124537421_124559486, Chr08_82174982_82248693, Chr07_23444137_23455147, Chr12_24249907_24258995) (Fig. 5D and Additional Table S13) in the muscle network diagram. As reported, RGS 9 knockout mice showed higher body weight and greater fat accumulation than wild-type mice[50, 51]. In addition, we also identified numerous well-known key markers in the muscle ceRNA networks, e.g., CDK16, CYC1, MYH3, HDC, and E2F1 for muscle cell growth and differentiation and FAAH, PLIN3, PNPLA2, and GAS7 for lipid synthesis and metabolism[52–54], suggesting that these ceRNA networks may play critical roles during muscle growth and development processes and providing ideas for further research studying muscle development and intramuscular fat deposition.
The differentiation of adipocytes is a highly complex biological process. While regulated by a series of transcription factor cascades, various cytokines can also participate in initiating and regulating the differentiation of adipocytes and maintaining cell characteristics through complex signal transduction. Adipose tissue is also an active endocrine organ that can secrete some hormones or cytokines (e.g., resistin, tumor necrosis factor-α and complement D) to participate in the immune response and treatment of obesity or cardiovascular diseases and other physiological and pathological processes[55, 56]. However, these functions may be inseparable from the regulation of circRNAs. For example, in our study, Chr15_78588493_78637426, produced by TLK1 as the host gene, was abundantly expressed in adipose tissue, starting at 30 d and decreasing after peaking at 150 d (Fig. 3E and Additional TableS8). Tousled-like kinase 1 (TLK1) is a serine/threonine protein kinase that is implicated in chromatin remodeling, DNA replication and mitosis[57]. Combined with the joint analysis results of SETM and ceRNA, it was found that among the predicted target miRNAs of circRNA-TLK1, only miRNA-204 was significantly expressed in the STEM analysis of miRNA. Among the target mRNAs of miRNA-204, only 4 mRNAs (LALMA4, COL6A1, ABCD2, ELMOD3) were significantly expressed in mRNA STEM analysis(Additional TableS13). It has been reported that lamininα4 (LAMA4) is located in the extracellular basement membrane that surrounds each individual adipocyte and affects the structure and function of adipose tissue in a depot-specific manner. Compared with wild-type (WT) mice, LAMA4−/− mice showed higher energy expenditure at room temperature and when exposed to a cold challenge. In addition, the mice had decreased adipose tissue mass and altered lipogenesis in a depot-specific manner[58, 59]. Interestingly, some scientists believe that obesity is an inflammatory condition that is associated with increased extracellular matrix (ECM) gene expression, whereas the collagen 6A1(COL6A1) gene is the primary gene in the ECM[60]. The ATP-binding cassette transporterABCD2 (D2) is a peroxisomal protein whose mRNA is highly expressed in fat and upregulated during adipogenesis. There are also articles reporting that D2 can not only promote the oxidation of long-chain monounsaturated fatty acids (LCMUFAs) but also inhibit the accumulation of dietary erucic acid (EA, 22:1ω9) in fat cells[61, 62]. These findings indicate that the same circRNA may play different roles in fat development and function.
Interestingly, combined with STEM analysis, we identified more ceRNA network pairs in the liver than in muscle and fat. Among these transcripts, many target mRNAs were determined to be related to fat metabolism and transport, such as ACSL3, ACSL5, APOA4, and PLIN4, and FADS, HACD3, and ACACA, for unsaturated fatty acid synthesis, indicating that these circRNAs are involved in fatty acid synthesis, metabolism and transport during the development of liver tissue. These results indicate that these ceRNA networks are active throughout tissue development. Although these in silico results should be investigated further in vivo, these results imply functional relatedness or a regulatory relationship between circRNAs and miRNAs during tissue development.