Identification and Characterization of CircRNAs of Two Pig Breeds as a New Biomarker in Metabolism-Related Diseases

Adipose tissue is a type of highly specialized loose connective tissue, which stores a high content of triglycerides and lipid-soluble substances. Surplus energy produced by the body during metabolism is generally stored in adipose tissue. If the capacity of adipose tissue to absorb and store lipids decreases, lipids may be then accumulated in non-adipose tissues, thus leading to obesity and imbalance in energy metabolism, and further leading to obesity-related diseases such as type II diabetes, atherosclerosis, insulin resistance, nonalcoholic fatty hepatitis, cardiovascular diseases, etc [1, 2].. In fact, health issues related to lipid metabolism disorders have received extensive attention, but the specific molecular mechanism remains unclear.

Circular RNA (circRNA) is a kind of endogenous circular molecule generated by alternative splicing. It forms a covalently closed continuous loop without 5’ and 3’ ends, thus it cannot be degraded by ribonuclease. In mammalian cells, circRNAs exist in large amounts and have properties such as good stability, conservativity and tissue specificity. It has been shown that circRNAs are related to Alzheimer’s disease [3], colorectal cancer and ovarian cancer [4], hepatocellular carcinoma [5], cardiovascular diseases [6], central nervous system disorders [7], etc., but their role in lipid metabolism remains unclear.

Recent research suggests that in eukaryotic cells circRNAs can work as miRNA sponges to regulate gene expression [8, 9]. Currently, there are evidences showing that miRNA can regulate the expression of genes associated with adipogenic differentiation and lipid metabolism [10, 11]. For example, in human preadipocytes and mouse 3T3-L1 cells, miR-143 promotes adipocyte differentiation by targeting MAPK1 [12]. Chen et al. reported that overexpression of miR-425-5p would inhibit the proliferation and differentiation of pig intramuscular preadipocytes and lead to down-regulated expression of PPARγ, FABP4 and FASN associated with fat formation [13]. In high fat-induced obese mice, miR-199a-3p regulates adipocyte differentiation by directly targeting 3’UTR of SCD, and it can change fatty acid composition by decreasing the proportion of unsaturated fatty acids [14]. However, it remains unknown whether these miRNAs involved in lipid metabolism are regulated by circRNAs.

Pigs are anatomically and physiologically similar to humans. Particularly, pigs and humans have similar cardiovascular, urinary, dermatological and digestive systems, and brown fat disappears in both after birth [15]. Therefore, pig, as an animal model of disease, has important clinical values. Large White pig is a typical lean meat breed with lower contents of subcutaneous fat and intramuscular fat than Laiwu pig. Because of the difference in fat deposition between the two pig breeds, circRNAs in their subcutaneous adipose tissues were sequenced in this paper and the expression patterns of the circRNAs were investigated. Bioinformatics methods were also used to systematically analyze circRNAs that can regulate lipid metabolism. The results may provide a theoretical basis for further study on lipid metabolism, and provide potential therapy targets for lipid metabolism-related diseases.

All procedures involving animals were approved by the Animal Care and Use Committee at Institute of Animal Sciences (Chinese Academy of Agricultural Sciences), where the experiments were conducted. All experiments were performed in accordance with the relevant guidelines and regulations set by the Ministry of Agriculture of the People’s Republic of China.

Experimental animals include two pig breeds: Large White pig and Laiwu pig (a local breed) with difference in fat deposition. Three female pigs of each breed were used. All pigs were raised under the same conditions in the same environment in Daqian Agriculture and Animal Husbandry Co., Ltd. (Laiwu City, China). They were fed according to the current nutritional needs (National Research Council, NRC, 1998). At slaughter age, their body weights were similar within species, and they were all healthy and in good body condition. After slaughter, the subcutaneous adipose tissue of the longissimus dorsi muscle of each animal was sampled.

Total RNAs from each sample were isolated with TRIzol reagent (Invitrogen Life Technologies, Carlsbad, USA) according to the manufacturer’s instructions. The purity and concentration of the total RNAs were determined with NanoDrop ND-1000. RNA integrity was tested by denaturing agarose gel electrophoresis.

Sequencing libraries were generated using a TruSeq Stranded Total RNA LT Kit with Ribo-Zero TM Gold (RS-122-2301) for Illumina according to the manufacturer’s recommendations. Ribo-Zero rRNA Removal Kit was used to remove ribosomal RNAs from the total extracted RNAs. RNase R was used to remove linear RNAs. Then, breaking reagent was added to break circRNAs into fragments. With the circRNA fragments as templates, random hexanucleotide primers were used to synthesize first-strand cDNAs. Subsequently, second-strand cDNAs were synthesized with second-strand cDNA synthesis mix. When second-strand cDNAs were synthesized, dUTP was used instead of dTTP and different adapters were ligated to the 5′ and 3′ ends of the second strand (The Index Primer were shown in Table 1). UNG enzyme was used to digest the cDNA strand containing dUTP, and only the cDNA strand with different adapters was kept. cDNAs were purified using cDNA purification kit. After purification, end repair was performed, poly (A) tail was added to cDNA and sequencing adapter was ligated. Subsequently, size selection was performed. Then, PCR amplification was performed. The quality of the library was examined with Agilent 2100 Bioanalyzer. After ensuring the quality of the library, Illumina HiSeqTM 2500 (or other sequencing equipment) was used for sequencing.

Paired-end sequencing was adopted. After the sequencing data were obtained, NGSQCToolkit (v2.3.3) was used for quality control and adapter trimming [16]. On this basis, low-quality bases and N-bases were trimmed from the reads. Finally, high-quality clean reads were obtained for subsequent analysis. BWA [17] was used to align the clean reads to reference gene sequences. CIRI was used to identify circRNAs. The alignment results (SAM format) were scanned to search for paired chiastic clipping (PCC) and paired-end mapping (PEM) signals, as well as GT-AG splicing signals. All the sequences with junction sites were realigned to reference gene sequences using dynamic programming algorithm to ensure the reliability of the identified circRNAs. DEseq2 [18] was used for differential expression analysis of circRNAs. Comparing treatment group and reference group, |fold Change| ≥2 and p-value < 0.05 were used to indicate significant difference in circRNA expression. Then, the number of up/down-regulated circRNAs was obtained.

Gene ontology (GO) (http://www.geneontology.org/) is an international standard for gene function classification. It classifies gene functions into three categories including molecular function, biological process and cellular component. Pathway enrichment analysis can help determine the main metabolic pathways and signaling pathways in which differentially expressed genes are involved. Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg) database [19] is a main public database for metabolic analysis and regulatory network research. To further explore the main biological functions of differentially expressed genes, we determined the GO terms and signaling pathways in which differentially expressed genes were enriched on the basis of hypergeometric distribution. If q-value ≤ 0.05, the enrichment was significant.

Total RNAs were reverse transcribed to synthesize cDNAs using PrimeScript RT Reagent Kit (Perfect Real Time; TaKaRa, Osaka, Japan), which were then subject to quantitative real-time PCR (qRT-PCR) analysis on an Applied Biosystems 7500 Fast Real-Time PCR System (Roche, Basel, Switzerland; v 2.0.5) with SYBR Green qPCR SMix (ROX; Roche). qRT-PCR results were used for verifying the gene expression levels revealed by above sequencing analysis. Six circRNAs (circRNA_26852, circRNA_15067, circRNA_25555, circRNA_11897, circRNA_23437 and circRNA_14707), three miRNAs (ssc-miR-486, ssc-miR-27a and ssc-miR-27b-3p) and three mRNAs (ABHD5, SMAD7 and SCD) were selected for amplification by specific divergent primers (Table 2). MiRNA reverse primers are universal primers that come with the QIAGEN kit. With glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene, 5S, Beta-actin (ACTB) as internal references of circRNA, miRNA, mRNA, respectively, 2^-ΔΔCt method was used to calculate the relative expression levels of gene between samples. t-test was performed to statistically analyze the relative expression levels.

miRanda software [20] was used to perform miRNA target prediction for circRNAs. In combination with our previous results on miRNAs and mRNAs in subcutaneous adipose tissues of Large White pig and Laiwu pig [21, 22], miRanda [20] and Targetscan [23] were used to predict the target genes of differentially expressed miRNAs. Differentially expressed mRNAs among miRNA target genes were selected and Cytoscape [24] was used to construct the circRNA-miRNA-mRNA interaction network. ClueGO [25] and KOBAS [26] were used to perform GO and KEGG enrichment analyses, respectively, of the differentially expressed genes related to the differentially expressed circRNAs.

All data were presented as mean values ± standard deviations. When comparisons were made, a Student’s t-test was performed. p < 0.05 was considered as statistically significant.

Laiwu pig and Large White pig have difference in fat deposition. By comparing the expression of circRNAs between two pig breeds, we identified circRNAs that might potentially regulate fat deposition in them and further explored the underlying molecular mechanism. A total of 29, 763 circRNAs were identified in the subcutaneous adipose tissues of Laiwu pig and Large White pig. Differential expression analysis was performed (Laiwu pig vs. Large White pig, denoted as L_PX vs. D_PX). With |fold change| ≥ 2.0 and p-value < 0.05 as criterion (Fig. 1), 275 differentially expressed circRNAs (70 up-regulated circRNAs and 205 down-regulated circRNAs) were obtained. The up-regulated circRNAs included 3 exonic circRNAs, 17 intergenic circRNAs, 49 sense-overlapping circRNAs and 1 antisense circRNAs (Fig. 2A), while the down-regulated circRNAs included 7 exonic circRNAs, 21 intergenic circRNAs, 174 sense-overlapping circRNAs and 3 antisense circRNAs (Fig. 2B).

According to our previous results on miRNAs and mRNAs in subcutaneous adipose tissues of Large White pig and Laiwu pig [21, 22] as well as the target miRNAs of circRNAs, six differentially expressed circRNAs were further verified by qRT-PCR. The results show that in subcutaneous adipose tissue of Laiwu pig, circRNA_26852, ssc-miR-486 and SCD were significantly up-regulated, while circRNA_15067, circRNA_23437, circRNA_14707, circRNA_11897, ssc-miR-27a, ssc-miR-27b-3p, ABHD5 and SMAD7 were significantly down-regulated. circRNA_25555 tended to be down-regulated but the difference was not statistically significant (Fig. 3).

CircRNAs can work as miRNA sponges and regulate miRNAs. Analysis of circRNA-miRNA interaction can help reveal the functions and functioning mechanisms of circRNAs as miRNA sponges. miRNA-circRNA interaction pairs were sorted according to p-values and the top 300 miRNA-circRNA interaction pairs with smaller p-values were selected for the construction of circRNA-miRNA interaction network (Fig. 4).

The expression levels of circRNAs determined by qRT-PCR were consistent with those determined by sequencing analysis. circRNA_26852 and circRNA_11897 that were up-regulated and down-regulated, respectively, in the subcutaneous adipose tissue of Laiwu pig were selected for subsequent analysis. ssc-miR-874 and ssc-miR-486 that can bind circRNA_26852 as well as ssc-miR-27a and ssc-miR-27b-3p that can bind circRNA_11897 (Table 3) were selected according to bioinformatics analysis and our previous results [21, 22] on miRNAs and mRNAs in the subcutaneous adipose tissue of two pig breeds. Subsequently, circRNA-miRNA-mRNA interaction networks were constructed (Fig. 5). Four miRNAs (ssc-miR-874, ssc-miR-486, ssc-miR-27a, and ssc-miR-27b-3p) in Fig. 5 were differentially expressed in the subcutaneous adipose tissues of Large White pig and Laiwu pig (Table 4). Other miRNAs in Fig. 5 were not differentially expressed but might also be potential targets of circRNA_26852 and circRNA_11897. All the mRNAs in Fig. 5 were differentially expressed in the subcutaneous adipose tissues of Large White pig and Laiwu pig [21]. qRT-PCR results of some genes including SCD, ABHD5 and SMAD7 were consistent with our previous sequencing results of differentially expressed mRNAs.

In order to deeply explore the functions of target genes, GO and KEGG pathway enrichment analyses were performed. The target genes of circRNA_26852 are mainly enriched in biological processes related to fat deposition and lipid metabolism, such as regulation of protein kinase D signaling, regulation of triglyceride catabolic process, negative regulation of lipid storage, phosphatidic acid biosynthetic process, etc. (Fig. 6A). Pathway enrichment analysis indicates that the target genes of circRNA_26852 are enriched in biosynthesis of unsaturated fatty acids, fatty acid metabolism, PPAR signaling pathway, TGF-beta signaling pathway, and disease-related signaling pathways (such as Insulin signaling pathway, AGE-RAGE signaling pathway in diabetic complications, Huntington’s disease, etc.) (Fig. 6B). The target genes of circRNA_11897 are also enriched in biological processes related to lipid metabolism such as fatty acid biosynthetic process, MAPK cascade reaction, ERK1 and ERK2 cascade reaction, cell proliferation, etc. (Fig. 7A). Pathway enrichment analysis shows that the target genes of circRNA_11897 are enriched in PPAR signaling pathway, AMPK signaling pathway, non-alcoholic fatty liver disease (NAFLD), Alzheimer’s disease, pathways in cancer, etc. (Fig. 7B). Differentially expressed miRNAs might target these differentially expressed genes in the pathways and then change their signal transduction. Since circRNAs can work as miRNA sponges, circRNAs might inhibit the effect of miRNAs on the target genes and then regulate fat deposition and metabolism-related diseases.

This study analyzed, for the first time, the expression of circRNAs in the subcutaneous adipose tissues of Large White pig and Laiwu pig. The aim was to identify potential circRNAs related to adipogenic differentiation and lipid metabolism. A total of 275 differentially expressed circRNAs were identified in the subcutaneous adipose tissue. Among them, 70 were up-regulated and 205 were down-regulated in the subcutaneous adipose tissue of Laiwu pig. The expressions of six circRNAs were further validated by qRT-PCR analysis. miRNA target prediction was performed for these differentially expressed circRNAs. The host genes of differentially expressed circRNAs and target genes of up-regulated circRNA_26852 and down-regulated circRNA_11897 were subject to GO and KEGG enrichment analysis. The regulatory role of circRNAs in fat deposition and lipid metabolism was explored from the perspective of host genes and competing endogenous RNAs. Our results may provide potential therapy targets for obesity and metabolic syndrome.

CircRNAs can absorb miRNAs associated with adipogenic differentiation and lipid metabolism. For the differentially expressed circRNAs, miRanda was used to detect some differentially expressed miRNAs that can bind them. It was found that ssc-miR-874 and ssc-miR-486 are targets of circRNA_26852, and ssc-miR-27a and ssc-miR-27b-3p are targets of circRNA_11897. Guo et al. reported that the expression of miR-874 is down-regulated in the adipose tissue of Japanese Black Wagyu cattle compared with that in Holstein cattle, thus miR-874 might regulate fat deposition in cattle [27]. In the subcutaneous adipose tissue of Laiwu pig, the expression of miR-874 was up-regulated. Japanese Black Wagyu cattle and Laiwu pig both have strong ability to deposit fat in their bodies, but the expression levels of miR-874 in them are different. It is inferred that this might be related to species and location. Some research indicates that in the subcutaneous adipose tissue of HIV-infected patients the expression of some miRNAs is abnormal and among them the expression of miR-874 is up-regulated [28]. Then, the research infers that changes in the expression levels of these miRNAs might lead to lipid metabolism disorders related to HIV. In addition, miR-874 might promote insulin synthesis [29]. Therefore, it is inferred that miR-874 is related to fat deposition and fat metabolism. Prats-Puig et al. found that the expression of miR-486 is improved in the plasma of obese children before adolescence and it is related to fat content, adiponectin, etc [30].. The expression of miR-486 in these obese children is consistent with that in Laiwu pig. Some report shows that miR-486 can inhibit the transcription factor, i.e., forkhead box O1 (Fox o1) [31], and plays a role in insulin functioning and triglyceride metabolism [32]. Consequently, it is inferred that miR-486 might be involved in lipid metabolism and associated with diabetes and obesity. Since ssc-miR-874 and ssc-miR-486 are targets of circRNA_26852, circRNA_26852 might participate in adipogenic differentiation and lipid metabolism. Wang et al. found that overexpression of miR-27a in pig adipocytes could promote lipolysis and the production of more glycerin and fatty acids [33]. Some report indicates that miR-27a can inhibit adipocyte differentiation by targeting PPARγ [34]. Song et al. found that miR-27b-3p can regulate fatty acid metabolism by targeting PPARγ and then play a role in the maturation and development of pig oocytes [35]. circRNA_11897 can bind miR-27a and miR-27b-3p, which might lead to the up-regulation of the target genes of miR-27a and miR-27b-3p and then the regulation of adipogenic differentiation and lipid metabolism. In sum, it is inferred that circRNAs may work as competing endogenous RNAs and regulate lipid metabolism, which may provide potential prevention targets for metabolism-related diseases.

Competing endogenous RNAs (ceRNAs) are RNAs with miRNA recognition elements. ceRNAs and mRNAs can competitively bind the same miRNA [36], thereby regulating each other. This mechanism is called “ceRNA” hypothesis. The target genes of circRNAs were analyzed to explore in depth the role of circRNAs in fat deposition and fat metabolism. GO and KEGG enrichment analyses reveal that the target genes of circRNA_26852 are enriched in biological processes related to fat deposition and lipid metabolism such as regulation of triglyceride catabolic process, negative regulation of lipid storage, phosphatidic acid biosynthetic process, etc. Pathway enrichment analysis shows that the target genes of circRNA_26852 are enriched in PPAR signaling pathway, TGF-beta signaling pathway, and disease-related signaling pathways (such as insulin signaling pathway AGE-RAGE signaling pathway in diabetic complications, Huntington’s disease). Zhang et al. reported that miR-140-3p can regulate adipocyte differentiation by targeting TGF-β1 [37]. In addition, TGF-β1 can phosphorylate SMAD3. The activated SMAD3 can bind C/EBPs to inhibit its transcriptional activity and thus play a role in the early stage of adipocyte differentiation [38]. PPARs are members of receptor superfamily and they include three subtypes: PPARα, PPARβ and PPARγ. PPARα gene can promote the hydrolysis of triglycerides into free fatty acids, activate β oxidation in liver and muscle, and decrease the synthesis of fatty acids and glycerin [39]. PPARγ is a key transcription factor regulating the transformation of mesenchymal stem cells into adipocytes [40]. The target genes of circRNA_11897 are also enriched in biological processes related to lipid metabolism such as fatty acid biosynthetic process, MAPK cascade reaction, ERK1 and ERK2 cascade reaction, cell proliferation, etc. MAPK, as serine/threonine protein kinase, can phosphorylate many protein kinases and nuclear transcription factors. MAPK cascade reaction is thus one of the important signaling pathways in cell. Four MAPK signaling pathways have been identified in mammalian cells: p38 MAPK pathway, ERK1/2 signaling pathway, ERK5/BMK1 pathway and JNK pathway. Zhao et al. found that activated p38 MAPK can promote fat synthesis in bone marrow mesenchymal stem cells [41]. In early stage of adipogenic differentiation, activated ERK signaling pathway can promote the expression of C/EBP and PPARγ and promote adipocyte differentiation [42]. In later stage, ERK1/2 phosphorylates and inactivate PPARγ, thus adipocyte differentiation is inhibited [43]. Therefore, the regulatory role of ERK1/2 in adipocyte differentiation is determined by the accurate time of ERK phosphorylation [44]. It is thus inferred that circRNA might regulate the expression of related genes in signaling pathway by binding miRNA. In this way, circRNA plays a regulatory role in fat deposition and lipid metabolism.

Analysis of circRNA-miRNA-mRNA interaction network shows that separating the synthetic ligands of α-β-hydrolase domain containing protein 5 (ABHD5, target gene of ssc-miR-486) from perilipin-1 (PLIN1) can directly activate adipose triglyceride lipase (ATGL) and promote lipolysis [45]. Our previous study shows that the expression of ABHD5 is down-regulated in the subcutaneous adipose tissue of Laiwu pig [21], whereas the expression of ssc-miR-486 is up-regulated in the subcutaneous adipose tissue of Laiwu pig [22]. Since circRNA_26852 can work as miR-486 sponge, it is inferred that circRNA_26852–miR-486–ABDH5 axis can regulate lipid metabolism. SCD gene is the target gene of ssc-miR-27b-3p. SCD is a key enzyme in lipogenesis and it can catalyze the transformation of saturated fatty acid (SFA) into monounsaturated fatty acid (MUFA) and play an important role in biosynthesis of fatty acid [46]. There is a lack of SCD1 gene in mice, thus feeding the mice with sugar cannot induce the synthesis and secretion of triglycerides [47]. The expression of SCD is up-regulated in the subcutaneous adipose tissue of Laiwu pig [21], while the expression of ssc-miR-27b-3p is down-regulated [22]. Bioinformatics analysis shows that circRNA_11897 can work as miR-27b-3p sponge. It is inferred that circRNA_11897 might play a role in regulating adipogenic differentiation and lipid metabolism. SMAD is a kind of intracellular protein and it can change extracellular signal into intracellular growth factor receptor, which then activate the transcription of downstream genes in nucleus [48, 49]. SMAD7 can inhibit the activation of SMAD2 and SMAD3 by TGFβ [50], and participate in adipocyte differentiation [51]. The expression of SMAD7 gene, as the target gene of ssc-miR-874, is down-regulated in the subcutaneous adipose tissue of Laiwu pig [21], while the expression of ssc-miR-874 is up-regulated [22]. Since circRNA_26852 can work as ssc-miR-874 sponge, it is inferred that circRNA_26852 might regulate the expression of target genes by binding ssc-miR-874 and further regulate adipocyte differentiation and lipid metabolism.

In sum, this study used circRNA sequencing and bioinformatics technique to analyze, for the first time, the expression of circRNAs in the subcutaneous adipose tissues of Large White pig and Laiwu pig. CircRNA-miRNA target prediction was performed and circRNA-miRNA-mRNA interaction network was constructed. By GO and KEGG pathway analyses, we identified pathways enriched in adipogenic differentiation and lipid metabolism, as well as disease-related pathways. It is inferred that circRNAs might regulate adipogenic differentiation and lipid metabolism. This research provides a theoretical basis for further study on fat deposition mechanism and provides potential therapy targets for metabolism-related diseases.

X.Y.M. conceived and designed the study and wrote the paper. A.L. performed the experiment, analyzed data and wrote the paper. W.L.H., X.X.Z. and L.L.X performed the experiment and interpreted the data. All the authors read and approved the final manuscript. This work was supported by a grant from the Agricultural Science and Technology Innovation Program (ASTIPIAS05), the National Basic Research and Development Program of China (973 Program) (No. 2015CB943100) and the Basic Research Fund for Central Public Research Institutes of CAAS (No. 2013ywf-zd-2).

The authors declare that they have no competing interests.