Spatiotemporal regulation of circular RNA expression during the developmentof skeletal muscle, subcutaneous fat, and liver in Ningxiang pigs

Background In recent years, thousands of different circular RNAs (circRNAs) have been identied through comparative analysis of different pig breeds. However, very few studies have investigated the spatiotemporal expression patterns of circRNA during organ development, which is crucial for functional and evolutionary analysis. In this study, we systematically analyzed circRNAs associated with fatty acid metabolism in the three main organs (muscle, fat and liver) a tfour growth time points (30 d, 90 d, 150 d and 210 d after birth) for Ningxiang pigs, awell-known native pig breed in China known for its excellent meat quality. We identied 61,683 circRNAs and analyzed their molecular characteristics, potential functions, and interactions with miRNAs. The circRNAs exhibited notable spatiotemporal specicity in the form of dynamic expression. At 90 d, the number and complexity of differentially expressed circRNAs in the liver was the highest.We used STEM (Short Time-series Expression Miner) to perform time series analysis on differentially expressed circRNAs and observed that there were multiple model spectra that were signicantly enriched with time changes in all three organs. All the signicantly changed circRNAs, miRNAs and mRNAs obtained from STEM analysis of each tissue were utilized to conrm that the circRNA and mRNA have a targeting relationship with miRNA. Notably, we screened 2384 (fat), 2672 (muscle) and 9187 (liver) circRNA-miRNA-mRNA network pairs related to developmental time changes. Finally, we used uorescence in situ hybridization (FISH) to demonstrate the dynamic subcellular localization of circRNA during development. These data suggest that circRNA is highly abundant in muscle, fat and liver and is dynamically expressed in a spatiotemporal manner. These results suggest that circRNAs play important roles in fat metabolism in Ningxiang pigs. Our ndings also provide new ideas and perspectives for the study of Ningxiang pigs' lipid metabolism regulation pathway and meat quality traits. and at


Background
In recent years, thousands of different circular RNAs (circRNAs) have been identi ed through comparative analysis of different pig breeds. However, very few studies have investigated the spatiotemporal expression patterns of circRNA during organ development, which is crucial for functional and evolutionary analysis.

Results
In this study, we systematically analyzed circRNAs associated with fatty acid metabolism in the three main organs (muscle, fat and liver) a tfour growth time points (30 d, 90 d, 150 d and 210 d after birth) for Ningxiang pigs, awell-known native pig breed in China known for its excellent meat quality. We identi ed 61,683 circRNAs and analyzed their molecular characteristics, potential functions, and interactions with miRNAs. The circRNAs exhibited notable spatiotemporal speci city in the form of dynamic expression. At 90 d, the number and complexity of differentially expressed circRNAs in the liver was the highest.We used STEM (Short Time-series Expression Miner) to perform time series analysis on differentially expressed circRNAs and observed that there were multiple model spectra that were signi cantly enriched with time changes in all three organs. All the signi cantly changed circRNAs, miRNAs and mRNAs obtained from STEM analysis of each tissue were utilized to con rm that the circRNA and mRNA have a targeting relationship with miRNA. Notably, we screened 2384 (fat), 2672 (muscle) and 9187 (liver) circRNA-miRNA-mRNA network pairs related to developmental time changes. Finally, we used uorescence in situ hybridization (FISH) to demonstrate the dynamic subcellular localization of circRNA during development.

Conclusion
These data suggest that circRNA is highly abundant in muscle, fat and liver and is dynamically expressed in a spatiotemporal manner. These results suggest that circRNAs play important roles in fat metabolism in Ningxiang pigs. Our ndings also provide new ideas and perspectives for the study of Ningxiang pigs' lipid metabolism regulation pathway and meat quality traits.

Background
After years of breeding work, the lean meat rate of pigs has been considerably improved, but the improvement of living standards has also prompted consumers to seek pork with higher meat quality. Chinese indigenous breeds are known for their 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 well-known for their excellent meat quality and have become one of the four best-known local breeds in China. It has been reported that unsaturated fatty acids (UFAs) in Ningxiang pork account for up to 59.6% of all fatty acids [2].Moreover, He etal. [3] found that the serum of Ningxiang pigs contains more unsaturated lipids than that of lean pigs under the same environment and feeding conditions. Polyunsaturated fatty acids (PUFAs) are indispensable essential fatty acids for the body, while long-chain (≥ C20) polyunsaturated fatty acids (LC-PUFAs) are highly biologically active forms of PUFAs, and their biological functions have become a heavily researched topic in international research [4,5].This characteristic is a valuable 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].Therefore, as a research model, Ningxiang pigs are of considerable signi cance for studies exploring genes related to meat quality traits and their expression patterns.
With the development of high-throughput sequencing technology, it has been shown that biological processes are regulated not only by protein-encoding RNAs (mRNAs) but also by noncoding RNAs (ncRNAs), such as microRNAs (miRNAs) and circular RNAs (circRNAs). CircRNAs are a new group of noncoding RNAs characterized by the existence of circular structures in which the 3' ends are covalently bound to the 5' ends. Functional circRNAs have been demonstrated to function as cytoplasmic microRNA sponges and nuclear transcriptional regulators, indicating that these RNAs participate in the regulation of gene expression [11].There have been many reports indicating a large difference between fatty and lean pig breeds regarding the in uence of circRNAs on pork quality.In a study of subcutaneous adipose tissue of Large White and Laiwu pigs, 275 differentially expressed circRNAs were identi ed, in which the target genes of circRNA_26852 and circRNA_11897 were found to be involved in adipocyte differentiation and lipid metabolism [12].The liver is a primary location for metabolic processes, such as lipid digestion, absorption, decomposition, synthesis and transportation. Adipose tissue is the primary energy storage organ and energy source of animals [13]. The amount of intramuscular fat (IMF) is related to the number of adipocytes and the ability to accumulate fat. Intramuscular fatty acids, especially polyunsaturated fatty acids, are important precursors of pork avor substances that indirectly or directly affect indicators of meat quality, such as pH, water-holding capacity, tenderness, mouth feel, and meat color [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 have rarely been reported.
In this study, we rst compared circRNA expression levels among three major lipid metabolic organs/tissues of Ningxiang pigs, namely, longissimus dorsi muscle (muscle), subcutaneous fat (fat) and liver, atvarious developmental stages (30,90,150, and 240 d after birth). The molecular characteristics, spatiotemporal expression patterns, potential functions and targeted miRNAs of the organ samples 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 circRNA-miRNA-mRNA interactions that signi cantly changed over time. Furthermore, screening circRNAs that affect fat metabolism throughout the developmental stages of pigs could signi cantly help to elucidate their biological and metabolic functions and relative importance.

Ethics Statement
The protocol for pig handling and treatment was in accordance with the recommendations of the 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
Twelve male purebred Ningxiang pigs were used in this study (all pigs were half-sibs). The experimental pigs were all provided by the Ningxiang pig breeding farm of Dalong Animal Husbandry Technology Co. Ltd. in Hunan Province, China. The 12pigswere examined at each of the four developmental stages in this study, including 30 d (30 d after birth), 90d(90 d after birth),150 d (150 d of age) and 210 d (210 d of age). Three biological repeats were employed per time point. The pigs were housed under standard environmental conditions (including unregulated room temperature and natural light) and fed three times per day on a standard diet with free access to water. Next, 3 healthy individuals with no disease and similar body weights were randomly selected in the 4 age groups from 30 d to 210 d in accordance with standard procedures. The pigs were fasted the following day. Next, the pigs were sacri ced by exsanguination under iso urane anesthesia at the surgical plane (4.5% tidal volume by mask). The longissimus dorsi muscle, subcutaneous fat, and liver samples were collected within 30 min after the pigs were sacri ced, immediately frozen in liquid nitrogen and stored at -80 °C until analysis.

Library Preparation and Illumina HiSeq X TenSequencing
The RNA-seq transcriptome strand library was prepared following the TruSeq TM Stranded Total RNA kit from Illumina (San Diego, CA) using 5 µg of total RNA. Brie y, ribosomal RNA (rRNA) depletionin stead of poly(A) puri cation was performed using an IlluminaRibo-Zero Magnetic kit, rst employing fragmentation by a fragmentation buffer. Second, rst-stranded cDNA was synthesized with random hexamer primers. Next, we removed the RNA template and synthesized a replacement strand, incorporating dUTP in place of dTTP to generate ds cDNA. AMPure XP beads were utilized to separate the ds cDNA from the secondstrand reaction mix. Asingle 'A' nucleotide was added to the 3' ends of these blunt fragments to block them from ligating to one another during the adapter ligation reaction. Finally, 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 cation using Phusion DNA polymerase (NEB) for 15 PCR cycles. After quanti cation by TBS380, the paired-end RNA-sequencing library was sequenced with Illumina HiSeq X Ten(2 × 150-bp read length). In addition, 3 µg of total RNA was ligated with sequencing adapters with the TruSeq™ 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 cation by TBS380, deep sequencing was performed by Shanghai Majorbio Bio-Pharm Biotechnology Co., Ltd. (Shanghai, China).

Identi cation of circRNAs
The CIRI2 (CircRNA Identi er 2) tool was used to identify circRNAs. This tool scans SAM les twice and collects su cient information to identify and characterize circRNAs. Brie y, during the initial 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 2scans the SAM alignment again to detect additional junction reads and performs further ltering to eliminate false positive candidates that may result 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. Signi cantly differentially expressed(DE) circRNAs were extracted with |log2FC| > 1 and P-values ≤ 0.05 by Deseq2.

Differential Expression Analysis and Functional Enrichment
To identify DEGs (differentially expressed 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 employed 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 analyses, was performed to identify which DEGs were signi cantly enriched in GO terms and metabolic pathways at Pvalues ≤ 0.05 compared with the whole-transcriptome background. GO functional enrichment and KEGG pathway analysis were performedusing 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. This analysis produces a color for backgrounds with signi cant pro les, while producing a white background for pro les that are not signi cant. 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-values ≤ 0.05). The corrected P-values are sorted from small to large.

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

RT-PCR
To validate the identi ed circRNAs in Ningxiang pigs, total RNA was extracted using an Animal Total RNA Kit (Tiangen, China) and treated with ribonuclease R (RNase R). cDNA was synthesized by reverse transcription using a 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 situhybridization (FISH)
Muscle samples xed with 4% paraformaldehyde were cut sagittally into thick sections and dewaxed in para n sections in water. Next, the sections were boiled in repair solution for 10 min. Proteinase K (20 µg/mL) was added dropwise, allowed to react with the samples at 37 °C for 30 min, and washed away by rinsing with PBS for 5 min 3 times. Then, the prehybridization solution was added and incubated at 37 °C for 1 h. The Chr01_143206410_143210729 probe (5'-GCTTCCTGTTTTTTACTTGGGCTGTTAG-3') (General Biol, China) containing a hybridization solution was added at a concentration, and hybridization was performed at 37 °C overnight. The nuclei stained with DAPI (Servicebio, Wuhan, China) appeared blue under ultraviolet excitation, and positive expression was indicated by red uorescence. The images were observed using an upright uorescence microscope (Nikon, Tokyo, Japan).

Statistical analysis
Differential gene expression analysis based on the negative binomial distribution and used the DESeq2, which a method for differential analysis of count data, using shrinkage estimation for dispersions and fold changes to improve stability and interpretability of estimates. This enables a more quantitative analysis focused on the strength rather than the mere presence of differential expression. Expression level difference analysis is a statistical inference process to determine the occurrence of differential expression of all detected genes. Due to the large number of genes involved, statistical tests need to be performed multiple times (the number of genes to be tested). In order to control the false discover rate, the pvalue obtained from statistical tests is need corrected, which perform multiple correct tests. Then the p-adjust value smaller than 0.05 was signi cant. With EDSeq2, the Wald test is commonly used for hypothesis testing when comparing two groups. A Wald test statistic is computed along with a probability that a test statistic at least as extreme as the observed value were selected at random. This probability is called the p-value of the test. Use software Goatools to perform GO enrichment analysis on the genes in the gene set, so as to obtain the main GO functions of the genes in the gene set. The method of use is Fisher's exact test. When the corrected P value (p-adjust) < 0.05, the GO function is considered to be signi cantly enriched. Use the R script to perform KEGG PATHWAY enrichment analysis on the genes in the gene set. The calculation principle is the same as the GO function enrichment analysis. When the corrected P value (p-adjust) < 0.05, the KEGG PATHWAY function is considered to be signi cantly enriched. RT-PCR data were analyzed by using unpaired 2-tailed Student's t tests; p < 0.05 (*), p < 0.01 (**). Data were expressed as mean ± S.E.M (standard error of mean).

Identi cation and Characteristicsof circRNAs
To perform comprehensive pro ling of circular RNA in pigs, we performed 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 dafter birth) in pigs. In this study, after low-quality and adapterpolluted reads were rst removed from the raw data, the Q30 value of the clean reads of the three organs/tissues at 30, 90, 150, and 210 d was observed to be greater than 95% (Table 1). Among the samples, as 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 was used for sequence segmentation alignment to accurately identify circRNAs. We identi ed 61,683 unique circRNAs in all assessed biological tissues (Additional Table S1).Analysis of chromosome distribution indicated that the identi ed circRNAs were transcribed widely and unevenly on the chromosome. Compared with other chromosomes, most of the circRNAs to be identi ed were distributed on chromosome 1, chromosome 6 and chromosome 13, which is in keeping with the characteristics of circRNAs reported by others [17].In addition, we also detected circRNAs in the sex chromosomes (Fig. 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-69.7%) were exons followed by intergenic circRNAs (approximately 16.1-17.8%), and only a small portion (approximately 12.8-15.3%) were located in introns. The length analysis showed that most circRNAs in the three tissues were longer than 3000 bp (Fig. 1B). We further constructed Venn diagrams of the identi ed circRNAs in three tissues and at four time points (Fig. 1D-1E). We found that a total of 9563 circRNAs were coexpressed in the three tissues, but the number of circRNAs uniquely expressed in the liver was considerably greater than that in muscle and fat. At the same time, it was found that a total of 9202 circRNAs were coexpressed in the four time periods. However, the number of circRNAs identi ed at 90 d was the largest, and the number of circRNAs that were uniquely expressed at 90 d was notably greater than that at other time points. CircRNAs also showed tissue-speci c and time-speci c expression in skeletal muscle, fat, and liver during the development of Ningxiang pigs. between any two groups. A total of 1786 circRNAs were de ned as DECs among these comparison groups ( Fig. 2A).
Next, we attempted to evaluate the potential functions of differentially expressed circRNAs (DECs) between two closed time points. Compared with piglets (30 d after birth), among these DECs, 156 circRNAs were upregulated at 90d, 204 circRNAs were upregulated at 150 d, and 188 circRNAs were upregulated at 210 d (Fig. 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 the growth and development of skeletal muscle. KEGG pathway analysis showed that compared with the lactation period (30 d), the biological functions of these three closed time points were different. For example, in skeletal muscle, signi cantly upregulated pathways at 90d are primarily involved in such processes as propanoate metabolism, valine, leucine and isoleucine degradation, and starch and sucrose metabolism; the signi cantly upregulated pathways at 150dare primarily involved in such processes as ubiquitin-mediated proteolysis, cellular senescence, and starch and sucrose metabolism; the signi cantly upregulated pathways at 210d are primarily related to genetic information processing, including such processes as RNA degradation, ubiquitin-mediated proteolysis, and protein processing in the endoplasmic reticulum. However, the downregulated pathways of these three groups of closed time points are all related to lipid metabolism and transport, including such processes as fatty acid biosynthesis, primary bile acid biosynthesis, and unsaturated fatty acid and peroxisome biosynthesis(p < 0.05, Additional Fig. 1 and Additional Table   S3).Overall, studies of the host genes of DECs at various time points indicate that the circularization of these genes may be important for skeletal muscle function. Among these circRNAs, we also found that 55were identi ed as common DE genes throughout the muscle development process (Fig. 2C and Additional TableS4).These common DE genes play important roles in development. KEGG pathway analysis indicated that these common DE genes were mainly enriched in lysosomes, cellular senescence, the VEGF signaling pathway, and ABC transporters (Fig. 2D).
To study the role of circRNAs in the transformation of skeletal muscle function, in addition to 30 d vs 90 d, we added 90d vs 150 d and 150 d vs 210 d comparisons. The results represent circRNAs that were up-or downregulated only between two consecutive time points during muscle growth. Between 90d and 150 d, the signi cantly upregulated pathways were primarily involved in such processes as longevity regulation, EGFR tyrosine kinase inhibitor resistance, and lysine degradation, and the downregulated pathways were primarily involved in such processes as valine, leucine and isoleucine degradation, ubiquitinmediated proteolysis, and regulation of actin cytoskeleton. Compared with the results obtained at 150d, the pathways upregulated at 210d were determined to beprimarily related to rheumatoid arthritis, antigen processing and presentation, and cytokine-cytokine receptor interactions, and downregulated pathways were observed to be primarily related to such metabolic pathways as starch and sucrose metabolism, ketone body synthesis and degradation, and valine, leucine and isoleucine degradation (p < 0.05, Additional Fig. 1 and Additional Table S3).These up-and downregulated circRNAs at different time periods may be associated with regulating the initiation or termination of growth and/or physiological processes at speci c developmental stages.
Next, we focused on the changing rules of the 20 most abundant circRNAs in skeletal muscle. Notably, 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, and PRKAR1A) ( Fig. 2E and Additional Table S5).

Expression Pattern of circRNAs in the Adipose Tissue Growth and Development of NingxiangPigs
Adipose tissue is primarily composed of fatty acids (FAs) and triglycerides, which play important roles in the quality of meat [19]. However, the type and proportion of fatty acids (FAs) are also affected by such factors 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 (Fig. 3A).
Among these closed time points, we detected 274 circRNAs upregulated at 90 d, 392 circRNAs upregulated at 150 d, and 396 circRNAs upregulated at 210 d compared with 30 d (Fig. 3B and Additional Table S6).KEGG analysis showed that compared with the lactation period (30 d), the upregulated pathways at these three closed time points were signi cantly enriched in the AMPK signaling pathway. Thesehost genes primarily includePIK3R1, IGF1R, PPP2R3A, PPARG, RPS6KB1, ACACA, and PPP2R3A, which encode proteins associated with adipocyte development and lipid anabolism. However, the downregulated KEGG pathways of these three groups of closed time points are all related to such pathways as complement and coagulation cascades, fatty acid biosynthesis, and steroid hormone biosynthesis. In addition, between 90d and 150d, signi cantly upregulated pathways are primarily involved in such pathways as synthesis and degradation of ketone bodies, synthesis of basal transcription factors, and acute myeloid leukemia, and downregulated pathways are primarily involved in such pathways as complement and coagulation cascades, steroid hormone biosynthesis, and primary bile acid biosynthesis. Compared with 150d, the upregulated pathways at 210d are primarily related to such metabolic pathways as the citrate cycle (TCA cycle), lysine degradation and glutathione metabolism, and downregulated pathways are primarily related to such pathways as protein processing in the endoplasmic reticulum, starch and sucrose metabolism, and adherens junctions(p < 0.05, Additional Fig. 1 and Additional TableS3).
Compared to piglets at30 d, 256 circRNAs were identi ed as common DE genes during the entire development process of adipose tissue, indicating that the KEGG pathways were primarily enriched in the complement and coagulation cascades, lysine degradation, the PPAR signaling pathway, the TGF-beta signaling pathway, the AMPK signaling pathway, and the PI3K-Akt signaling pathway (Fig. 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),and6-phosphofructo-2-kinase/fructose-2,6-bisphosphatase 1 (PFKFB1) (Fig. 3E and Additional TableS8).

Spatiotemporal Dynamic Expression Pattern of circRNAs in the Liver of Ningxiang Pigs
The liver is the primary organ for triglyceride metabolism, and it plays a role as a regulator of energy metabolism together with adipose tissue. The liver 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 the liver, we performed expression pro ling of circRNAs in the liver at postnatal days 30, 90, 150, and 210.
This pro ling enabled the evaluation of dynamic changes in circRNA expression from lactation to fattening and the identi cation of circRNAs associated with fatty acid anabolism and transport.
We clustered the 3,362 differential circRNAs identi ed in the ve comparison groups (Fig. 4A). Between the two closed time points, we detected 598 circRNAs upregulated at 90 d, 580 circRNAs upregulated at 150 d, and 572 circRNAs upregulated at 210 d compared with 30 d (Fig. 4C and Additional TableS9). KEGG analysis showed that between the two closed time points in liver tissue, the signi cantly upregulated pathways were primarily enriched in lipid metabolism and amino acid metabolism, including such pathways as tryptophan metabolism, primary bile acid biosynthesis, valine, leucine and isoleucine degradation, and retinol metabolism. However, the signi cantly downregulated pathways were primarily related to signal transduction, such as the Notch signaling pathway, ECM-receptor interactions, and the MAPK signaling pathway (p < 0.05, Additional Fig. 1 and Additional TableS3).
In addition, between two consecutive time points, the number of differentially expressed circRNAs in the 90 d vs. 150 d and 150 d vs. 210 d groups was notably lower than that in the 30 d vs. 90 d groups (Fig. 4B). The KEGG pathway analysis also showed that the circRNAs expressed between 30d and 90 d participated in more biological pathways or more complex biological processes (p < 0.05, Additional Fig. 1 and Additional TableS3). These results indicate that during liver development, circRNAs may play important roles 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, notably exceeding the number of muscle and adipose tissues (Fig. 4D, Additional TableS10). KEGG pathway analysis also indicated that these common DE genes were primarily enriched in such pathways as complement and coagulation cascades, ECM-receptor interactions, tryptophan metabolism, the Notch signaling pathway, and steroid hormone biosynthesis. Next, we observed that several of the most abundant circRNAs in the liver originate from protein-coding genes that play pivotal roles in lipid synthesis and metabolism.(e.g.,Scaffold180_2562774_2570664,Scaffold155_4811844_4828350,and Chr01_143206410_143210729) ( Fig. 4E and Additional TableS11).

Construction of the circRNA-miRNA-mRNA coexpression networks through time-series analysis
According to the dynamic expression patterns of the circRNAs across the four developmental 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 (Fig. 5A-5C). Finally, we observed colored modules and only considered the largest module. Functional enrichment analysis showed that the largest module in the three tissues was enriched in a variety of biological processes. Some of these processes were related to signal transduction, cell growth and death, amino acid metabolism, and lipid metabolism(Additional Table S12).
Previous studies have shown that circRNAs act as miRNA sponges and indirectly regulate gene transcription. To explore whether there are ceRNAs with functional correlations or regulatory relationships between circRNAs that are strongly correlated with time points, we performed STEM analysis of miRNAs and mRNAs(AdditionalFigure2). In the circRNA STEM analysis, we selected the circRNAs in the colored modules to predict the miRNA and obtained the intersection with the miRNAs that are signi cant in the miRNA STEM analysis. Next, target analysis was performed with mRNAs that are signi cant in mRNA STEM analysis. We screened 2384 (fat), 2672 (muscle) and 9187 (liver) circRNA-miRNA-mRNA network pairs related to developmental time changes (Additional Table S13). We found that there were far more ceRNA network pairs in the liver than in muscle and fat. More importantly, we identi ed numerous ceRNA network pairs related tothe cell cycle, cell differentiation, fatty acid biosynthesis metabolism markers, such as CDK16, MYH3, ACSL3, ACSL5, APOA4, and FADS. (Fig. 5D-5F). This result 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
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 primers were designed for each circRNA, and both cDNA and gDNA (negative control) were employed as templates for PCR ampli cation ( Fig. 6 and Additional Table S14). Finally, among the 11circRNAs, 10 were successfully con rmed (90.90%), suggesting the reliability of our circRNA identi cation results. In addition, as an alternative visualization method, we employed 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 the nucleus and cytoplasm. Although signal strength cannot be utilized as a quantitative measure of circRNA abundance, we detected greater signal strength at 90 d and 150 d, which is consistent with the Chr01_143206410_143210729 sequencing data in muscle (Fig. 6C and Fig. 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].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] identi ed 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 signi cantly involved in muscle development and function. These researchers alsoconstructed the rst public S. scrofacircRNA database [17]. Morten et al. [25]pro led 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 identi ed a highly complex regulation pattern of thousands of circRNAs with a distinct spatiotemporal expression pro le, 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 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 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 identi ed 61,683 circRNAs in all evaluated tissues. The identi ed 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 identi ed 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 identi ed tissue-speci c circular RNAs, which provided important information regarding their functions. For example, we identi ed 13871, 6274 and 20646 tissue-speci c circRNAs in muscle, fat and liver, respectively (Fig. 1C). These ndings suggest that speci cally expressed circRNAs may play speci c roles at speci c times during the development of tissues or organs. This possibility is consistent with the ndings of other studies, indicating that circRNAs have tissue-/stage-speci c expression differences [25][26][27].
Subsequently, the identi cation of circRNAs by grouping statistics indicated that 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. The identi ed 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 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 regulatory mechanisms [30,31]. Among these genes, TNNI2 (fast-twitch skeletal muscle isoform, named TNI-fast) is a speci c protein of the fast muscle ber 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 bers [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 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 the properties of meat. Notably, pathway analysis of the genes giving rise to circRNAs peaking at 30d indicated a signi cant 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 signi cant 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 uctuations, 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 6phosphofructo-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 veri cation. 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 uctuations, 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 identi ed 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][53][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 signi cantly expressed in the STEM analysis of miRNA. Among the target mRNAs of miRNA-204, only 4 mRNAs (LALMA4, COL6A1, ABCD2, ELMOD3) were signi cantly 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-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 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 ndings indicate that the same circRNA may play different roles in fat development and function.
Interestingly, combined with STEM analysis, we identi ed 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.