To perform comprehensive profiling 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 adaptor-polluted reads were firstly 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 identified 61,683 unique circRNAs in all assessed biological tissues (Additional Table S1). Analysis of chromosome distribution indicates that the identified circRNAs are transcribed widely and unevenly on the chromosome. Compared with other chromosomes, most of the circRNAs identified 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 identified (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 identified 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 identified 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-specific expression and time-specific expression in the skeletal muscle, fat, and liver during Ningxiang pig development.
Table 1
Sequencing basic data of four different developmental stages of Ningxiang pig three organs/tissues
Terms
|
NX30d
|
NX90d
|
NX150d
|
NX210d
|
Muscle
|
Fat
|
Live
|
Muscle
|
Fat
|
Live
|
Muscle
|
Fat
|
Live
|
Muscle
|
Fat
|
Live
|
Raw reads number
|
65606883
|
57278618
|
56500489
|
50488816
|
53672934
|
55531984
|
53929255
|
45120232
|
45030983
|
56107588
|
58206111
|
51897636
|
Clean reads number
|
64137048
|
56517482
|
55574608
|
48897656
|
52575561
|
54755520
|
52912480
|
44451010
|
44381475
|
55237321
|
57266380
|
51069717
|
Clean reads rate
|
97.76%
|
98.67%
|
98.36%
|
96.85%
|
97.96%
|
98.60%
|
98.11%
|
98.52%
|
98.56%
|
98.45%
|
98.39%
|
98.40%
|
Clean Q30 bases rate
|
95.42%
|
95.02%
|
95.14%
|
95.35%
|
95.48%
|
95.55%
|
95.35%
|
95.48%
|
95.28%
|
95.63%
|
95.28%
|
95.52%
|
Mapped reads
|
120386916
|
106157588
|
102827570
|
91329628
|
97298423
|
103499271
|
99877320
|
81793231
|
83368071
|
104372240
|
105818225
|
95480783
|
Mapping rate
|
93.85%
|
93.90%
|
92.51%
|
93.38%
|
92.55%
|
94.51%
|
94.38%
|
92.00%
|
93.92%
|
94.48%
|
92.41%
|
93.48%
|
The values represent the reads and proportion that were compared to those in the Ningxiang pig reference genome (PRJNA531381) using the Hisat2 program.
|
3.2 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 profile 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 five 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 identified by the DEseq2 with P-value < 0.05 and had fold-change > 2 between any two groups. Total of 1786 circRNAs were defined 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 significantly upregulated pathways at 90d are mainly involved in Propanoate metabolism, Valine, leucine and isoleucine degradation, and Starch and sucrose metabolism, etc; The significantly upregulated pathways at 150d are mainly involved in Ubiquitin mediated proteolysis, Cellular senescence, and Starch and sucrose metabolism, etc; The significantly 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 identified 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, significantly 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 up- and downregulated circRNAs at different time periods may mean regulating the start/stop of growth and/or physiological processes at specific 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).
3.3 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 identified in the five 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 significantly 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, significantly 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 identified 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).
3.4 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 profiling of circRNAs in liver at postnatal days 30, 90, 150,and 210. This profiling allowed evaluation of dynamic changes in circRNA expression from lactation to fattening and identification of circRNAs associated with fatty acid anabolic and transport.
We clustered the 3,362 differential circRNAs identified in the five 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 significantly 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 significantly 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 identified 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 Time-series Analysis
According to their dynamic expression patterns across the four development stages, we found that all the identified circRNAs were classified into 4, 4 and 5 cluster profiles in the muscle, fat, and liver, respectively. Each tissue included at least 5 significantly enriched model profiles(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 figure 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 significant in the miRNA STEM analysis; Then target analysis was performed with mRNAs that are significant 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.
3.6 Verification of Identified circRNAs
In order to validate the circRNAs identified 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 amplification (Figure 6 and Additional Table S14). Finally, among the 11 circRNAs, 10 were successfully confirmed (90.90%), suggesting the reliability of our circRNAs identification result. In addition, as an alternative visualization method, we used fluorescence 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).