Spatiotemporal Regulation of Circular RNA Expression During Ningxiang Pig Skeletal Muscle, Subcutaneous Fat, and Liver Development

Background: In recent years, thousands of different circular RNA (circRNAs) have been discovered by comparative analysis of different pig breeds. However, very few studies have been conducted on the spatiotemporal expression pattern of circRNA during organ development which is crucial for functional and evolutionary analysis. Results: Here, we systematically analyzed circRNAs associated with fatty acid metabolism in the three main organs (muscle, fat and liver) for four growth time points (30d, 90d, 150d and 210d after birth) of Ningxiang pigs, a famous native pig breeds in China, known for its excellent meat quality. We identied 61,683 circRNAs and analyzed their molecular characteristics, potential functions, and interaction with miRNA. The circRNAs showed striking spatiotemporal specicity in the form of dynamic expression. At 90d, the number and complexity of differentially expressed circRNAs in the liver is the most obvious. We used STEM (Short Time Series Expression Miner) to perform time series analysis on differentially expressed circRNAs, and found that there are multiple model spectra that are signicantly enriched with time changes in all the three organs. All the signicantly changed circRNAs, miRNAs and mRNAs obtained from STEM analysis of each tissue were used to conrm that the circRNA and mRNA have a targeting relationship with miRNA. It is worth noting that we screened out 2384 (fat), 2672 (muscle) and 9187 (liver) circRNAs-miRNA-mRNA network pairs related to developmental time changes. Finally, we used uorescence in situ hybridization (cid:0) FISH (cid:0) to show the dynamic subcellular localization of circRNA during development. Conclusion: To a certain extent, these data indicate that circRNA is highly abundant in muscle, fat and liver, and is dynamically expressed in a spatiotemporal manner. These results suggest that circRNAs play an important role in fat metabolism of Ningxiang pigs. Our ndings also provide new ideas and perspectives for the study of Ningxiang pigs' lipid metabolism regulation pathway and meat quality traits. identied at four . We that a total of 9563 circRNAs were co-expressed in the three tissues, but the number of circRNAs uniquely expressed in the liver was far more than that of muscle and fat. At the same time, it was found that a total of 9202 circRNAs were co-expressed in four time periods. However, the number of circRNAs identied at 90d was the largest and the number of circRNAs that were uniquely expressed at 90d was far more than other time points. CircRNA also showed tissue-specic expression and time-specic expression in the skeletal muscle, fat, and liver during Ningxiang development.


Introduction
After years of breeding work, the lean meat rate of pigs has been greatly improved, but the improvement of living standards has also prompted people to pursue higher pork meat quality. Chinese indigenous breeds are known for its quality meat production, and there have been some reports of using Chinese native pigs as a biological model to study meat quality traits [1]. Ningxiang pigs are famous for their excellent meat quality and become one of the four famous local breeds in China. It has been reported that unsaturated fatty acids (UFA) in Ningxiang pork are up to 59.6% [2]. Moreover, He et al. [3] found that the serum of Ningxiang pigs contains more unsaturated lipids compared with lean pigs under the same environment and feeding conditions. Polyunsaturated fatty acids (PUFA) are indispensable essential fatty acids for the body, while long-chain (≥ C20) polyunsaturated fatty acids (LC-PUFA) are highly biologically active forms of PUFA, and their biological functions have been a hot topic in international researches [4,5]. This is a preciously ne trait of Ningxiang pig. It has been reported that the synthesis of unsaturated lipids is regulated by multiple genes, such as peroxisome proliferator activated receptor alpha (PPARα), fatty acid binding protein 4 (FABP 4), lipoprotein lipase (LPL) [6], CCAAT enhancer binding protein alpha (C/EBPα) [7], patatin like phospholipase domain containing 3 (PNPLA3) [8], malic enzyme 1 (ME1) [9] and 2,4-dienoyl-CoA reductase 1 (DECR1) [10], etc. Therefore, as a research model, Ningxiang pigs are of great signi cance for exploring genes related to meat quality traits and their expression patterns.
With the advancement of high-throughput sequencing technology, it is shown that biological processes are not only regulated by protein-encoding RNA(mRNA), but also by non-coding RNA (ncRNA), such as microRNA (miRNA) and circular RNA (circRNA). Circular RNAs (circRNAs) are a new group of noncoding RNA characterized by the existence of circular structures where the 3' ends are covalently bound to the 5' ends. Functional circRNAs have been shown to function as cytoplasmic microRNA sponges and nuclear transcriptional regulators, indicating that these RNAs participate in regulation of gene expression [11]. There have been many reports that there is a big difference between fatty and lean pig breeds in the in uence of circRNAs on pork quality. In the study of subcutaneous adipose tissue of Large White and Laiwu pigs, 275 differentially expressed circRNAs are found, in which the target genes of circRNA_26852 and circRNA_11897 are found to be involved in adipocyte differentiation and lipid metabolism [12]. The liver is an main place for metabolic processes such as lipid digestion, absorption, decomposition, synthesis and transportation. Adipose tissue is the main energy storage organ and energy source of animals [13]. The amount of intramuscular fat IMF is related to the number of adipocyte and the ability to accumulate fat. Intramuscular fatty acids, especially polyunsaturated fatty acids, are important precursors of pork avor substances, which indirectly or directly affect the indicators for evaluating pork quality, such as pH, water-holding capacity, tenderness, mouthfeel, meat color, etc [14]. However, the expression patterns and potential regulatory mechanisms of circRNAs related to meat quality or lipid anabolism in these organs across different developmental stages of pigs are rarely reported.
In this study, we rst compared circRNA expression levels among three major lipid metabolic organs/tissues of Ningxiang pigs, longissimus dorsi muscle (muscle), subcutaneous fat (fat) and liver across different developmental stages (30,90,150, and 240 days after birth). Their molecular characteristics, spatiotemporal expression patterns, potential functions and targeted miRNAs were systematically analyzed. To explore the dynamic expression pattern of circRNAs, we adopted a time-series analysis by STEM and constructed a network diagram of circRNAs-miRNAs-mRNAs interaction signi cantly related to time changes. Furthermore, screening circRNAs that affects fat metabolism throughout the developmental stage of pigs would have signi cant contribution to understand their biological and metabolic functions and relative importance.

Ethics Statement
This protocol involving pig handling and treatments was carried out in accordance with the recommendations of Laboratory Animals-Guideline of Welfare and Ethics of China. The study was approved by the "Institutional Animal Care and Use Committee of Hunan Agricultural University".

Experimental Animals and Sample Collection
The 12 male purebred Ningxiang pigs used in this study (all pigs are half-sibs). The experimental pigs were all provided by Ningxiang pig breeding farm of Dalong animal husbandry technology Co. Ltd., in Hunan province, China. They are housed under standard environmental conditions (including unregulated room temperature and natural light) and fed three times on a standard diet with free access to water every day. Next, 3 healthy individuals with no disease and similar body weight were randomly selected in the 4 age groups from 30 days to 210 days in accordance with standard procedures. The pigs were fasted and humanely slaughtered the following day. The longissimus dorsi muscle, subcutaneous fat, and liver samples were collected within 30 minutes after slaughter, immediately frozen in liquid nitrogen and stored at -80 °C until analysis.

Library Preparation, and Illumina Hiseqxten Sequencing
RNA-seq transcriptome strand library was prepared following TruSeqTM stranded total RNA Kit from Illumina (San Diego, CA) using 5μg of total RNA. Shortly, ribosomal RNA (rRNA) depletion instead of poly(A) puri cation was performed by Ribo-Zero Magnetic kit and then fragmented by fragmentation buffer rstly.
Secondly rst-stranded cDNA was synthesized with random hexamer primers. Then we removed the RNA template and synthesized a replacement strand, incorporating dUTP in place of dTTP to generate ds cDNA. AMPure XP beads were used to separate the ds cDNA from the second strand reaction mix. A single'A' nucleotide was added to the 3' ends of these blunt fragments to prevent them from ligating to one another during the adapter ligation reaction. Lastly, multiple indexing adapters were ligated to the ends of the ds cDNA. Libraries were size selected for cDNA target fragments of 200-300 bp on 2% Low Range Ultra Agarose followed by PCR ampli ed using Phusion DNA polymerase (NEB) for 15 PCR cycles.After quanti ed by TBS380, paired-end RNA-sequencing library was sequenced with the Illumina HiSeqxten(2 × 150bp read length). In addition, 3μg of total RNA was ligated with sequencing adapters with TruseqTM Small RNA sample prep Kit (Illumina, San Diego, CA, USA). Subsequently, cDNA was synthesized by reverse transcription and ampli ed with 12 PCR cycles to produce libraries. After quanti ed by TBS380, deep sequencing was performed by Shanghai Majorbio Bio-Pharm Biotechnology Co., Ltd. (Shanghai, China).

Identi cation of circRNAs
The CIRI 2 (CircRNA Identi er 2) tool was used to identify circRNA. It scans SAM les twice and collects su cient information to identify and characterize circRNAs. Brie y, during the rst scanning of SAM alignment, CIRI 2 detects junction reads with PCC signals that re ect a circRNA candidate. Preliminary ltering is implemented using paired-end mapping (PEM) and GT-AG splicing signals for the junctions. After clustering junction reads and recording each circRNA candidate, CIRI 2 scans the SAM alignment again to detect additional junction reads and meanwhile performs further ltering to eliminate false positive candidates resulting from incorrectly mapped reads of homologous genes or repetitive sequences. Finally, identi ed circRNAs are output with annotation information. The expression level of each circRNA was calculated according to the Spliced Reads per Billion Mapping (SRPBM) method.

Differential Expression Analysis and Functional Enrichment
To identify DEGs (differential expression genes) between two different samples, the expression level of each transcript was calculated according to the transcripts per Kilobase of exon model per Million mapped reads (TPM) method. RSEM (http://deweylab.biostat.wisc.edu/rsem/) was used to quantify gene abundances. R statistical package software EdgeR (Empirical analysis of Digital Gene Expression in R,http://www.bioconductor.org/packages/2.12/bioc/html/edgeR.html) was utilized for differential expression analysis. In addition, functional-enrichment analysis including GO and KEGG were performed to identify which DEGs were signi cantly enriched in GO terms and metabolic pathways at P-value ≤ 0.05 compared with the whole-transcriptome background. GO functional enrichment and KEGG pathway analysis were carried out by Goatools (https://github.com/tanghaibao/Goatools) and KOBAS (http://kobas.cbi.pku.edu.cn/home.do).

Time-series Analysis
Time-series analysis was performed by STEM (Short time-series Expression Miner) [15]. Time-series analysis studies the dynamic behavior of gene expression and measures a series of processes that have strong correlations with time points. There is the color of the background statistically signi cant pro le, with white background that is not statistically signi cant pro le. Model pro les of the same color belong to the same cluster of pro les. Signi cantly enriched model pro les are indicated by different colors (Bonferroni-adjusted P-value ≤ 0.05). The corrected P value is sorted from small to large.

Analysis of ceRNAs Regulatory Network
To reveal the role and interactions among ncRNAs and mRNAs, we constructed a ncRNAs-mRNAs regulatory network. MiRanda was used to predict the circRNA-miRNA-mRNA pairs (score cutoff ≥ 160 and energy cutoff ≤ -20). The co-expression network of circRNA-miRNA-mRNA was constructed using the Cytoscape software (v3.2.1) to explore the function of key circRNAs [16].

RT-PCR
To validate the identi ed circRNAs in Ningxiang pig, total RNA was extracted using Animal Total RNA Kit (Tiangen, China) and treated with Ribonuclease R(RNase R). cDNA was synthesized by reverse transcription using RevertAid First Strand cDNA Synthesis Kit (Thermo Scienti c, USA). The divergent PCR primers were designed using the "out-facing" strategy to exclude linear mRNA from ampli cation (synthesized by Sangon Biotech Co., Ltd., Shanghai, China). PCR products were validated by Sanger sequencing.

Fluorescence in situ hybridisation (FISH)
Muscle samples xed by 4% paraformaldehyde were cut sagitally into thick sections and dewaxing para n sections to water, the sections were boiled in the repair solution for 10 min. Proteinase K (20 μg/mL) was added dropwise and digested at 37℃ for 30 min and wash with PBS for 5 min at 3 times. Then the prehybridisation solution was added and incubated at 37℃ for 1 h. The Chr01_143206410_143210729 probe (5'-GCTTCCTGTTTTTTACTTGGGCTGTTAG-3') (General Biol China containing a hybridization solution was added at a concentration, and hybridisation was performed at 37℃ overnight. The nucleus stained by DAPI (Servicebio, Wuhan, China) was blue under ultraviolet excitation, and the positive expression was red uorescence. The images were observed using a upright uorescence microscope (Nikon, Tokyo, Japan).

Identi cation and Characteristics of circRNAs
To perform comprehensive pro ling of circular RNA in pigs, we carried out total RNA sequencing in three major organs/tissues of lipid anabolic metabolism (fat, liver, and skeletal muscle) and across four developmental stages (30,90,150 and 210 days after birth) of pigs. In this study, after low-quality and adaptorpolluted reads were rstly removed from the raw data, the Q30 value of the clean reads of the three organs/tissues at 30, 90, 150, and 210 days was greater than 95% (Table 1). Among them, mapped with the Ningxiang pig genome (accession number PPJNA531381) assembled by our research group, the mapping rate of each sample exceeded 92%. In addition, HISAT is used for sequence segmentation alignment to accurately identify circRNA. We identi ed 61,683 unique circRNAs in all assessed biological tissues (Additional Table S1). Analysis of chromosome distribution indicates that the identi ed circRNAs are transcribed widely and unevenly on the chromosome. Compared with other chromosomes, most of the circRNAs identi ed were distributed on chromosome 1, chromosome 6 and chromosome 13, which is consistent with the characteristics of circRNA reported by others [17]. In addition, we also detected circRNA in the sex chromosomes ( Figure 1A). As shown in previous studies, most circRNAs contain multiple exons, and some circRNAs also retain introns. In this study, most of the circRNAs identi ed (approximately 68.6% to 69.7%) were exons, followed by intergenic (approximately 16.1% to 17.8%) circRNAs, and only a small portion (approximately 12.8% to 15.3%) were located in introns. The length analysis showed that most circRNAs in three tissues are longer than 3000 bp ( Figure 1B). We further constructed Venn maps of the identi ed circRNAs in three tissues and at four time points Figure 1D-1E . We found that a total of 9563 circRNAs were co-expressed in the three tissues, but the number of circRNAs uniquely expressed in the liver was far more than that of muscle and fat. At the same time, it was found that a total of 9202 circRNAs were co-expressed in four time periods. However, the number of circRNAs identi ed at 90d was the largest and the number of circRNAs that were uniquely expressed at 90d was far more than other time points. CircRNA also showed tissue-speci c expression and time-speci c expression in the skeletal muscle, fat, and liver during Ningxiang pig development. The values represent the reads and proportion that were compared to those in the Ningxiang pig reference genome (PRJNA531381) using the Hisat2 program

The Expression Pattern of circRNAs in the Muscle Growth and Development of Ningxiang Pigs
Intramuscular fat (IMF) content and fatty acid composition are important meat quality characteristics [18]. However, the molecular mechanisms regulating intramuscular fat accumulation and fatty acid composition are still not understood clearly. In order to understand the regulation of circRNA in skeletal muscle, we analyzed the expression pro le of circRNAs in skeletal muscle at 30, 90, 150 and 210 days after birth. This analysis can assess the dynamic changes in circRNAs expression from lactation to fattening, and identify circRNAs that may be related to intramuscular fat (IMF) deposition and fatty acid composition. First, we designed ve comparison groups (1: 30d vs. 90d; 2: 30d vs. 150d; 3: 30d vs. 210d; 4: 90d vs. 150d; and 5: 150d vs 210d) to discover the differential expression of circRNAs (DECs) during development. Differentially expressed circRNAs were identi ed by the DEseq2 with P-value < 0.05 and had fold-change > 2 between any two groups. Total of 1786 circRNAs were de ned as DECs among these comparison groups ( Figure 2A).
To evaluate the potential function of differentially expressed circRNAs (DECs) between two closed time points. Compared with the piglets (30d after birth), among these DECs, 156 circRNAs were upregulated at 90d, 204 circRNAs were upregulated at 150d, and 188 circRNAs were upregulated at 210d ( Figure 2B, Additional Table S2). Under the assumption that the functions of circRNAs are related to the functions of their host genes, we performed KEGG pathway analysis to predict the functions of DECs during growth and development of skeletal muscle. KEGG pathways analysis showed that compared with lactation period (30d), the biological functions of these three closed time points were different. For example, skeletal muscle signi cantly upregulated pathways at 90d are mainly involved in Propanoate metabolism, Valine, leucine and isoleucine degradation, and Starch and sucrose metabolism, etc; The signi cantly upregulated pathways at 150d are mainly involved in Ubiquitin mediated proteolysis, Cellular senescence, and Starch and sucrose metabolism, etc; The signi cantly upregulated pathways at 210d are mainly related to genetic information processing such as RNA degradation, Ubiquitin mediated proteolysis, and Protein processing in endoplasmic reticulum, etc; However, the down-regulated pathways of these three groups of closed time points are all related to lipid metabolism and transport such as Fatty acid biosynthesis, Primary bile acid biosynthesis, Biosynthesis of unsaturated fatty acids and Peroxisome, etc(p<0.05, Additional Figure 1 and Additional Table S3). Overall studies of the host genes of DECs at different time points indicate that the circularization of these genes may be important for skeletal muscle function. Among them, we also found that 55 circRNAs were identi ed as common DE genes throughout the muscle development process ( Figure 2C and Additional Table S4). These common DE genes play an important role in development. KEGG pathway analysis found that these common DE genes are mainly enriched in Lysosome, Cellular senescence, VEGF signaling pathway, and ABC transporters, etc ( Figure 2D).
In order to study the role of circRNAs in the transformation of skeletal muscle function, in addition to 30d vs 90d, we added 90d vs 150d and 150d vs 210d.
The results represent circRNAs that were up-or downregulated only between two consecutive time points during muscle growth. Between 90d and 150d, signi cantly upregulated pathways are mainly involved in Longevity regulating pathway, EGFR tyrosine kinase inhibitor resistance, and Lysine degradation, etc., and downregulated pathways are mainly involved in Valine, leucine and isoleucine degradation, Ubiquitin mediated proteolysis, Regulation of actin cytoskeleton, etc.; Compared with 150d, the up-regulated pathways at 210d are mainly related to Rheumatoid arthritis, Antigen processing and presentation, Cytokine-cytokine receptor interaction, and downregulated pathways are mainly related to metabolic pathways such as Starch and sucrose metabolism, Synthesis and degradation of ketone bodies, Valine, leucine and isoleucine degradation etc. (p<0.05, Additional Figure 1 and Additional Table S3).These upand downregulated circRNAs at different time periods may mean regulating the start/stop of growth and/or physiological processes at speci c developmental stages.
Next, we focused on the changing rules of top20 abundant circRNAs in skeletal muscle. It is worth noting that several of the most abundant circRNAs in skeletal muscle originate from protein coding genes with pivotal roles in skeletal muscle growth, intramuscular fat deposition, and intramuscular sugar metabolism (e.g. MYH2, MYH6, TTN, PFKFB1, TNNI2, MYBPC1, PRKAR1A) ( Figure 2E and Additional Table S5).

The Expression Pattern of circRNAs in the Adipose Tissue Growth and Development of Ningxiang Pigs
Adipose tissue is mainly composed of fatty acids (FA) and triglycerides, which plays an important role in the quality of meat [19]. However, the type and proportion of fatty acids (FA) are also affected by factors such as breed and genetic conditions [20,21].We clustered 3,116 differential circRNAs identi ed in the ve comparison groups to obtain an overview of DECs in adipose tissue ( Figure 3A).
Among these closed time points, we detected 274 circRNAs upregulated at 90d, 392 circRNAs upregulated at 150d, and 396 circRNAs upregulated at 210d, compared with 30d ( Figure 3B and Additional Table S6). KEGG analysis showed that compared with lactation period (30d), the up-regulated pathways of these three closed time points were signi cantly enriched in AMPK signaling pathway, these host genes are mainly generated from PIK3R1, IGF1R, PPP2R3A, PPARG, RPS6KB1, ACACA, and PPP2R3A, which encode proteins associated with adipocyte development and lipid anabolism. However, the down-regulated KEGG pathways of these three groups of closed time points are all related to Complement and coagulation cascades, Fatty acid biosynthesis, and Steroid hormone biosynthesis, etc. In addition, between 90d and 150d, signi cantly upregulated pathways are mainly involved in Synthesis and degradation of ketone bodies, Basal transcription factors Acute myeloid leukemia, and downregulated pathways are mainly involved in Complement and coagulation cascades, Steroid hormone biosynthesis, Primary bile acid biosynthesis; Compared with 150d, the upregulated pathways at 210d are mainly related to metabolic pathways such as Citrate cycle (TCA cycle), Lysine degradation and Glutathione metabolism, and downregulated pathways are mainly related to Protein processing in endoplasmic reticulum, Starch and sucrose metabolism, Adherens junction (p<0.05, Additional Figure 1 and Additional Table S3).
Compared to Piglets of 30d, 256 circRNAs were identi ed as common DE genes during the entire development process of adipose tissue indicated that the KEGG pathways were mainly enriched in Complement and coagulation cascades, Lysine degradation, PPAR signaling pathway, TGF-beta signaling pathway, AMPK signaling pathway, and PI3K-Akt signaling pathway, etc ( Figure 3C-3D and Additional Table S7). Among the most abundant circRNAs expressed in adipose tissue, several circRNAs derived from protein-coding genes that affect the key roles of adipose tissue development and fatty acid metabolism are highly expressed in adipose tissue (e.g. NADP-dependent malic enzyme (ME1), protein kinase A type 1-α regulatory subunit (PRKAR1A), and 6-phosphofructo-2-kinase/fructose-2, 6-bisphosphatase 1 (PFKFB1) ( Figure 3E and Additional Table S8).

Spatiotemporal Dynamic Expression Pattern of circRNAs in Liver of Ningxiang Pigs
The liver is the main organ for triglyceride metabolism, and it plays a role as a regulator of energy metabolism together with adipose tissue. It regulates the expression of various lipid metabolism-related genes and takes up extracellular free fatty acids to promote the transmembrane transport of fatty acids. To understand the regulation of circRNAs in liver, we performed expression pro ling of circRNAs in liver at postnatal days 30, 90, 150 and 210. This pro ling allowed evaluation of dynamic changes in circRNA expression from lactation to fattening and identi cation of circRNAs associated with fatty acid anabolic and transport.
We clustered the 3,362 differential circRNAs identi ed in the ve comparison groups ( Figure 4A). Between two closed time points, we detected 598 circRNAs upregulated at 90d, 580 circRNAs upregulated at 150d, and 572 circRNAs upregulated at 210d, compared with 30d ( Figure 4C and Additional Table S9). KEGG analysis showed that between two closed time points in liver tissue, the signi cantly upregulated pathways are mainly enriched in lipid metabolism and amino acid metabolism such as Tryptophan metabolism, Primary bile acid biosynthesis, Valine, leucine and isoleucine degradation Retinol metabolism, etc; However, the signi cantly downregulated pathways are mainly related to signal transduction, such as Notch signaling pathway, ECM-receptor interaction, and MAPK signaling pathway, etc (p<0.05, Additional Figure 1 and Additional Table S3).
In addition, between two consecutive time points, the number of differentially expressed circRNAs in the two groups of 90d vs. 150d and 150d vs. 210d is far less than that of 30d vs. 90d ( Figure 4B). The KEGG pathway analysis also shows that the circRNA expressed between 30d and 90d participates in more biological pathways or more complex biological processes (p<0.05, Additional Figure 1 and Additional Table S3). These results indicate that during liver development, circRNAs may play an important role in the regulation of amino acid metabolism, fatty acid metabolism and carbohydrate metabolism.
Compared to 30d, 559 circRNAs were identi ed as common DE genes during the entire liver development process, far exceeding the number of muscle and adipose tissues ( Figure 4D, Additional Table S10). KEGG pathway analysis also found that these common DE genes are mainly enriched in Complement and coagulation cascades, ECM-receptor interaction, Tryptophan metabolism, Notch signaling pathway, and Steroid hormone biosynthesis, etc. Next, we observed that several of the most abundant circRNAs in liver originate from protein coding genes with pivotal roles in lipid synthesis and metabolism. (e.g. Scaffold180_2562774_2570664, Scaffold155_4811844_4828350, and Chr01_143206410_143210729) ( Figure 4E and Additional Table S11).
3.5 Constructing of the circRNA-miRNA-mRNA Co-expression Networks through Timeseries Analysis According to their dynamic expression patterns across the four development stages, we found that all the identi ed circRNAs were classi ed into 4, 4 and 5 cluster pro les in the muscle, fat, and liver, respectively. Each tissue included at least 5 signi cantly enriched model pro les (Figure 5A-5C). Finally, we observed colored modules and only considered the largest module. Functional enrichment analysis showed that the largest module in three tissues are enriched in a variety of biological processes. Some of them were related to signal transduction, Cell growth and death, Amino acid metabolism, lipid metabolism, etc(Additional Table S12).
Previous studies have shown that circRNAs act as miRNAs sponges and indirectly regulate gene transcription. In order to explore whether there are ceRNAs with functional correlation or regulatory relationship between circRNAs that are strongly correlated with time points, we also added STEM analysis of miRNA and mRNA (Additional gure 2). In the circRNA STEM analysis, we selected the circRNAs in the colored modules to predict the miRNA and took the intersection with the miRNAs that are signi cant in the miRNA STEM analysis; Then target analysis was performed with mRNAs that are signi cant in mRNA STEM analysis. We screened out 2384 (fat), 2672 (muscle) and 9187 (liver) circRNAs-miRNAs-mRNAs network pairs related to developmental time changes (Additional Table S13). We found that there are far more ceRNA network pairs in the liver than muscle and fat. More importantly, we have discovered many ceRNA network pairs related to cell cycle, cell differentiation, fatty acid biosynthesis metabolism markers such as CDK16, MYH3, ACSL3, ACSL5, APOA4, FADS, etc. (Figure 5D-5F). It indicates that these ceRNA networks may play an important role in the regulation of adipogenesis and metabolism during the development of Ningxiang pigs.

Veri cation of Identi ed circRNAs
In order to validate the circRNAs identi ed from the circRNA-Seq analysis, reverse transcription RT-PCR and Sanger sequencing experiments were performed for 11 randomly selected circRNAs. A pair of divergent primerswere designed for each circRNA, and both cDNA and gDNA (negative control) were used as template for PCR ampli cation ( Figure 6 and Additional Table S14). Finally, among the 11 circRNAs, 10 were successfully con rmed (90.90 ), suggesting the reliability of our circRNAs identi cation result. In addition, as an alternative visualization method, we used uorescence in situ hybridization (FISH) for higher resolution subcellular localization of Chr01_143206410_143210729. At four time points of development, the Chr01_143206410_143210729 signal may be related to both nucleus and cytoplasm. Although signal strength cannot be used as a quantitative measure of circRNA abundance, we detected more signal strength in 90d and 150d, which is consistent with the Chr01_143206410_143210729 sequencing data in muscle( Figure 6C and Figure 4E).

Discussion
CircRNAs have been reported to be abundant in animal transcriptomes, and a large number of circRNAs have been identi ed in humans [22], mice [23], and nematodes [24]. The pig is an important farm animal that provides meat for humans and an important medical animal model, and research on pigs has become more and more abundant in recent years. Liang et al. [17] identi ed 149 circRNAs that may be related to muscle growth in the analysis of the nine organs and three developmental stages of the Guizhou minipig, and found that their host genes were signi cantly involved in muscle development and function. They also constructed the rst public S. scrofacircRNA database [17]. Morten et al. [25] pro le the expression of circRNA in ve brain tissues at up to six time-points during fetal porcine development, constituting the rst report of circRNA in the brain development of a large animal. An unbiased analysis reveals a highly complex regulation pattern of thousands of circRNA, with a distinct spatiotemporal expression pro le, suggesting the important function of circRNA in mammalian brain development. Accurate timing of gene expression is very important for the developmental stage. Because circRNAs are expressed in a highly spatiotemporal speci c 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 is still unclear. In order to study the spatiotemporal expression pattern of circRNA on lipid synthesis, metabolism and transport in Ningxiang pig tissues, we carried out total RNA sequencing across three organs (liver, skeletal muscle, and subcutaneous fat)at four developmental stages (30, 90, 150 and 240 days after birth). Compared with the reference genomes of other pig species, the results of using 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 were more than 92% (Table 1). We identi ed 61,683 circRNAs in all evaluated tissues. The identi ed circRNAs have 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 in comparison with the other chromosomes, and the largest number was derived from exons. In addition, the number of circRNAs identi ed at 90d was the largest, and the number of uniquely expressed circRNAs at 90d was far more than other time points. We also identi ed tissue-speci c circular RNA, which provided important information for their functions.
For example, we identi ed 13871, 6274 and 20646 tissue-speci c circRNAs in muscle, fat and liver, respectively( Figure 1C). These ndings suggest that speci cally expressed circRNAs may play speci c functions at speci c times during the development of tissues or organs. This is consistent with other studies, indicating that circRNA has tissue-/stagespeci c expression differences differences [25][26][27]. Subsequently, the identi cation of circRNAs by grouping statistics revealed 3,362, 1,786 and 3,116 circRNAs with signi cant differential expression in Ningxiang pig liver, muscle and adipose tissue, respectively.
Since the number of muscle bers 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 ber. Type I muscle bers are rich in mitochondria and have a higher oxidative metabolism, while type IIB bers have fewer mitochondria, and the glycolytic metabolism capacity is higher [29]. However, the role of circRNAs in skeletal muscle development is unknown. Of the identi ed muscle abundant expressed circRNAs for example, Chr05_81425549_81426478, Chr15_85761324_85761662, Scaffold180_2562774_2570664, Chr02_159842602_159843286, and Chr12_55773288_55775037 were highly expressed in muscle( Figure 2E and Additional Table S5). The contractile protein troponin I (TnI), a constituent protein of the troponin complex located on the thin laments 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-speci c and developmental regulations [30,31]. Among them, TNNI2 (Fast-Twitch Skeletal Muscle Isoform, named TNI-fast) is a speci c protein of fast muscle ber type. During muscle development, it was observed that Chr02_159842602_159843286, produced by host gene TNNI2, was abundantly expressed at 30d. It has been reported that the change of MYBPC1 abundance is related to the production of slow muscle bers [32,33]. Interestingly, we observed a peak expression of Chr05_81425549_81426478 transcribed by Myosin binding protein C1 MYBPC1 as host gene in the developing 90d pig muscle. Chr15_85761324_85761662, produced by host gene Titin TTN , was shown to be expressed at 30d and reached its peak at 210d. It has been reported that TTN is related to intramuscular fat deposition [34], and recent studies have found that circular RNA TTN acts as a miR-432 sponge to facilitate 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 MYH myosin superfamily, such as MYH2 and MYH6. Recent studies have shown that they are associated with skeletal muscle development and in uence the quality of pork [36]. These ndings suggest that many circRNAs exhibit different expression patterns during muscle development, and changes in the abundance of these circRNAs may affect muscle ber composition, thereby affecting meat character. Notably, pathway analysis of the genes giving rise to the circRNAs peaking at 30d reveals a signi cant predominance of genes associated with fatty acid biosynthesis, primary bile acid biosynthesis, unsaturated fatty acid biosynthesis and peroxisome. An impact of circRNAs on the biosynthesis of unsaturated fatty acids will be of great concern, which 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 more abundant in obese Rongchang pigs than in lean Landrace pigs. Furthermore, mRNA abundance changes of ME1 have the remarked signi cant positive correlation with adipocyte volume across the six adipose tissue types [37]. Chr01_176448874_176466186 is produced by host gene ME1, which is abundantly expressed in adipose tissue ( Figure 3E and Additional Table S8). The expression of Chr01_176448874_176466186 is shown to increase from 90d, peak at 150d and thereafter, with uctuations of 210d. It is speculated to be closely related to adipocyte development and adipodeposition. CircRNA 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( Figure 4E and Additional Table S11). 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. Its synthesis and degradation are catalyzed by the bifunctional enzyme 6-phosphofructo-2kinase/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 speculate 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 requires follow-up experimental veri cation. In addition, we found that Chr01_143206410_143210729 is abundantly expressed in all three organs but presented different expression patterns. For example, expression of Chr01_143206410_143210729 in developmental adipose is shown to increase from 30d, peak at 150d and thereafter, with uctuations, decrease towards the time of 210d. However, the expression of Chr01_143206410_143210729 in muscle and liver peaked at 90d and decreased gradually after that. Although Chr01_143206410_143210729's host gene SAFB like transcription modulator (SLTM ) has been rarely reported, its homologous family member SAFB1 can interact with and repress 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 which is induced during over nutrition to facilitate the conversion of glucose to fatty acids and triglycerides for the storage of the 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 an important role on additive, synergistic, or potentially antagonist functions among 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 circRNA is not completely clear, the current research shows that in addition to the function of circRNA to regulate host gene transcription, protein binding and translation, more research is to play the role of miRNA sponge [45,46]. Competing endogenous RNAs (ceRNAs) are widely present in muscle, fat, and liver. For example, circularRNA SNX29 sponges miR-744 to regulate 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) ( Figure 5D and Additional Table S13), in the muscle network diagram. As reported, RGS 9 knockout mice showed higher body weight and more fat accumulation, compared with wild-type mice [50,51]. In addition, we also found many well-known key markers in the muscle ceRNA networks, e.g., CDK16, CYC1, MYH3, HDC, E2F1, etc. for muscle cells growth and differentiation, FAAH, PLIN3, PNPLA2, GAS7, etc. for lipid synthesis and metabolism [52][53][54], suggesting that these ceRNA networks may play critical roles during muscle growth and development process, and provide ideas for studying muscle development and intramuscular fat deposition.
The differentiation of adipocytes is an extremely complex biological process. While being 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, which can secrete some hormones or cytokines (resistin, tumor necrosis factor-α, complement D, etc.) 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 30d and falling after peaking at 150d( Figure 3E and Additional Table S8). 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 signi cant in the STEM analysis of miRNA. Among the target mRNAs of miRNA-204, only 4 mRNAs (LALMA4, COL6A1, ABCD2, ELMOD3) are signi cant in mRNA STEM analysis(Additional Table S13). It has been reported that Laminin α4 (LAMA4) is located in the extracellular basement membrane that surrounds each individual adipocyte and affecting the structure and function of adipose tissue in a depot-speci c 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-speci c manner [58,59]. Interestingly, Some scientists believe that Obesity is an in ammatory condition t, which is associated with increased extracellular matrix (ECM) gene expression, whereas collagen 6A1(COL6A1) gene is the primary gene ECM [60]. The ATP binding cassette transporter, ABCD2 (D2), is a peroxisomal protein whose mRNA is highly expressed in the fat and upregulated during adipogenesis. There are also articles reporting that D2 can not only promote the oxidation of long-chain monounsaturated fatty acids (LCMUFA), but also inhibit the accumulation of dietary erucic acid (EA, 22:1ω9) in fat cells [61,62]. These analysis results indicate that the same cirRNA may play different roles in fat development and function.
Interestingly, combined with STEM analysis, we found more ceRNA network pairs in the liver, compared to muscle and fat. Among them, many target mRNAs are related to fat metabolism and transport, such as ACSL3, ACSL5, APOA4, PLIN4, etc., and FADS, HACD3, ACACA, etc., for unsaturated fatty acids 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 have always existed during tissues development. Although these in silico results should be investigated in vivo, these results imply functional relatedness or a regulatory relationship between circRNAs and miRNAs during tissue development.

Conclusions
In summary, this study provides a clear and accurate dynamic pro le of circRNA in the skeletal muscle, subcutaneous fat and liver based on our Ningxiang pig genome. Moreover, our study shows the role of circRNAs in lipid synthesis and metabolism to a certain extent, and provides high-quality resources and novel clue for studying the mechanism of the molecular regulation of pork quality traits.

Availability of data and materials
The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.

Ethics approval
All animal experiments in this study were approved by the "Institutional Animal Care and Use Committee of Hunan Agricultural University"(Changsha, China).

Consent for publication
Not applicable.

Con ict of interests
The authors have no con icts of interest to declare.           Detailed information can be found in supplement Table S20. Black arrows represent the PCR ampli cation orientation of the divergent primer. (C) Subcellular localization of Chr01_143206410_143210729 in the skeletal muscle at four development time points. Subcellular localization was visualized using panomics probes for high-resolution ISH and DAPI for nuclear localization.