Integrative analysis of circRNAs, miRNAs, and mRNAs profiles to reveal ceRNAs networks in chicken intramuscular and abdominal adipogenesis

Background Tissue-specific fat deposition is regulated by a series of complex regulatory mechanisms. Reports indicate that epigenetic regulators, such as circular RNAs (circRNAs), are crucial in diseases progression, animal development, metabolism, and adipogenesis. In this study, to assess the functional roles of circRNAs in adipogenesis and tissue-specific fat deposition, we comprehensively analyzed the Ribo-Zero RNA-Seq and miRNAs data during chicken intramuscular and abdominal adipogenic differentiation. Results circRNAs and miRNAs profiles during chicken adipogenic differentiation were found in adipocytes derived from various adipose tissues. It was also discovered that high levels of downregulated miRNAs potentially promote adipogenesis by activating their target genes which are associated with fatty acid metabolism and adipogenic differentiation. Through analysis of the correlation between the expression levels of circRNAs and adipogenic genes, as well as the dynamic expression patterns of circRNAs during adipogenic differentiation, several candidate circRNAs were identified. Moreover, competing endogenous RNA (ceRNAs) networks were constructed during chicken intramuscular and abdominal adipogenesis by combining miRNAs with mRNAs data. Several candidate circRNAs potentially influence adipogenesis by regulating miRNAs via PPAR and fatty acid metabolism-related pathways were identified, such as circLCLAT1, circFNDC3AL, circCLEC19A and circARMH1. Conclusion In conclusion, our findings reveal that circRNAs and the circRNA-miRNAs-mRNAs-ceRNAs network may play important roles in chicken adipocytes differentiation and tissue-specific fat deposition.


Background
Different resources and organizational microenvironment dictate various physiological and biochemical characteristics of adipose tissues [1][2][3]. Since preadipocytes derived from various tissues exhibit distinct proliferation speed and adipogenic potential, it is believed that transcriptional and post-transcriptional regulation mechanisms cause stage-and tissue-specific fat deposition in animals [4][5][6].
In this study, comprehensive circRNAs atlas during chicken adipogenic differentiation across various tissuesderived preadipocytes were provided. Future, we identified differentially expressed (DE) circRNAs, miRNAs and their dynamic expression patterns during adipocytes differentiation. Then, circRNAs and RNA-Seq data were combined to construct circRNAs-miRNAs-mRNAs ceRNA network involved in chicken adipogenesis. This study provides an in-depth understanding of epigenetic underlying regulation mechanisms of adipogenesis and tissues-specific fat deposition in birds.

Overview of circRNAs data by Ribo-zero RNA-Seq in chicken preadipocytes and adipocytes
Here, we analyzed 960 million raw paired-end (PE) reads generated by Ribo-Zero RNA-Seq data. After filtering with low-quality reads, we retained a total of 910 million clean reads. Unique clean reads of 88-94% were mapped to the chicken genome 5 (galGal5) ( Table 1). It was found that circRNAs were distributed in most of the chicken chromosomes. The expression levels of cir-cRNAs in different tissue-derived adipocytes from the adipose tissue were identified then analyzed to find a better understanding of the cirRNAs in chicken adipocytes. As shown in Fig. 1a, intramuscular preadipocytes (IM_Pre) and adipocytes (IM_Ad), abdominal preadipocytes (Ab_Pre), and adipocytes (Ab_Ad) were fell into two clusters, suggesting that the differences within the adipose tissues are potentially larger than that of adipocytes. Moreover, the Pearson correlation coefficient of circRNA expression between the cell samples within groups waved from 0.86 to 0.93, suggesting a reliable consistency of cell samples within groups. The bar plot showed that no apparent difference in global expression levels of circRNAs across different groups (Fig. 1b). To classify the length distribution of circRNAs, the resource of all circRNAs in different gene regions were compared, including exon, intergenic, and intron regions. Eventually, 787, 710, 980, 1083, 1406, 1426, 1447 and 1004 circRNAs were identified in eight groups respectively (Fig. 1c). Results showed that the majority of circRNAs (65%) were generated from exon regions, followed by intergenic regions (27%) and intron regions (8%) (Fig.  1d). As shown in Fig. 1d, most of the intersect circRNAs (60%) shared from different software were 200-1000 nt in length. Moreover, a considerable number of circRNAs which longer than 3000 nt were identified in the present study (Fig. 1d). Besides, most of circRNAs (94%) were found more than two exons (Fig. 1e). The circRNAs were distributed across 31 chromosomes (chromosomes 1-28, 33, Z, and W) of chicken with chromosomes 1, 2, 3 and 21 carrying the largest share (Fig. 1f, Fig. 2).

Identification of differentially expressed (DE) circRNAs during adipogenic differentiation across various groups
To explore functional roles of circRNAs in tissues-specific adipogenic differentiation, we analyzed and compared the expression levels of circRNAs between different groups. Eventually, we identified 17, 14, 20, and 9 differentially expressed (DE) circRNAs across different groups (IMPre_vs_IMAd, AbPre_vs_AbAd, AbPre_vs_IMPre and AbAd_vs_IMAb) (Fig. 3, Table 2). Most of DE circRNAs were specificially found in preadipocytes or adipocytes. In addition, several DE circRNAs showed different expression patterns during adipogenic differentiation between AbF and IMF groups (Table  S1). Among them, 7:22323550|22,327,655 was significantly downregulated in mature adipocytes when compared with preadipocytes in both IMF and AbF groups. Z:78636651|78,

Functional characterization of DE circRNAs
To investigate the potential functions of DE circRNAs, GO terms analysis of host genes of DE circRNAs were performed. We found that the molecular function (MF) of gene ontology (GO) analysis during adipogenesis were primarily enriched in lipid binding, phospholipase activity, and lipase activity, whereas cell component (CC) were mainly enriched in regulation of autophagy, glycerophospholipid biosynthetic process, phosphatidylinositol-mediated signaling and inositol lipid-mediated signaling (Fig. 4).

Validation of circRNAs
To confirm the circRNAs in chicken adipocytes, five DE circRNAs including circFNDC3AL, circHSPG2, cir-cCLEC19A, circARMH1, and circLCLAT1 were randomly selected for validation using PCR and Sanger sequencing ( Fig. 5a). Divergent and convergent primers were designed to amplify circRNAs back-spliced junction (BSJ) site and linear mRNAs. As showed in Fig. 5b, BSJ sites of circRNAs were amplified and validated by sanger sequencing. The amplified PCR products using convergent primers of cir-cRNAs were distinct both in cDNA and gDNA samples, whereas the divergent primers only amplified the cDNA samples. PCR products amplified using the divergent primers were further identified through Sanger sequencing (Fig. 5b). The sequencing results corresponded to Ribo-zero RNA-Seq data, which suggested that the cir-cRNAs data is reliable.

Candidate circRNAs related to adipogenic differentiation
Further, correlation analysis between circRNAs and their target genes based on the ceRNAs network was performed to investigate the crucial roles of circRNAs in the adipogenesis of chicken. As shown in Fig. 8a, the expression levels of circFNDC3AL and circHSPG2 were significantly and positively related to PPAR signaling pathway genes in intramuscular adipogenesis (r > 0.90, p < 0.01), such as peroxisome proliferator-activated receptor γ (PPARγ), fatty acid binding protein 4 (FABP4), Fig. 2 The circos plot of circRNAs in chicken adipocytes. The ring from the outer layer to the inner layer represent represents AbPre1, AbPre1, AbAd1, AbAd2, IMPre1, IMPre2, IMAd1and IMAd2 groups, respectively perilipin 2 (PLIN2), CD36 molecule (CD36). In contrast, the expression level of circLCLAT1 was significantly and negatively related to PPAR signaling pathway genes and positively related to lysocardiolipin acyltransferase 1 (LCLAT1) and fatty acid desaturase 2 (FADS2) (r > 0.95, p < 0.01) (Fig. 8a). Moreover, the expression levels of cir-cARMH1 were significantly and negatively related to those of genes in abdominal adipocytes (r > 0.92, p < 0.01), while circCLEC19A was significantly and positively related to those of genes (r > 0.90, p < 0.01) (Fig.  8b). To further investigate the potential roles of circRNAs in the fat tissue-specific deposition, the expression of cirRNAs was analyzed in different tissues. Our results demonstrated that the expression of these four candidate circRNAs was highly expressed in the abdominal fat tissue of chicken (Fig. 8c). In addition, the dynamistic expression patterns of candidate circRNAs were detected during tissue-specific adipogenetic differentiation. Our results showed that the expression levels of circFNDC3AL were decreased after 4 days of induction, while dramatically upregulated in the late stage of intramuscular adipogenesis (Fig. 8d). The expression levels of circCLEC19A were significantly increased during the late stage of abdominal adipogenesis (Fig. 8 d).
In contrast, the expression levels of circCLEC19A were dramatically decreased after 2 days' of adipogenetic induction, while slightly increased in the late stage of intramuscular adipogenesis. In addition, the expression levels of circARMH1 were decreased during the abdominal adipogenesis (Fig. 8d). Interestingly, we noticed that chicken circLCLAT1 was highly conservative among humans, mice and pigs (Fig. S1), implying that cir-cLCLAT1 potentially regulates animal adipogenesis or lipid metabolism.

Discussion
Preadipocytes-derived from abdominal fat (AbF) have previously been reported to exhibit higher adipogenic differentiation potential than the ones derived from intramuscular fat (IMF) in chickens [26]. Adipogenesis is controlled by a series of complex regulatory mechanisms, including transcription level, post-transcription level, among others [27][28][29]. Noncoding RNAs have notable effects on adipogenesis and lipid metabolism [30,31]. Several studies about noncoding RNAs on adipogenesis were geared towards mammals [11,12,30,31]. However, information on circular RNAs (circRNAs) in poultry adipogenesis is limited.
In this work, we analyzed Ribo-zero RNA-Seq and corresponding miRNAs data during chicken tissue-specific adipogenic differentiation. Three software (CIRI2, CIR-Cexplore2 and find_circ) were used to identify cir-cRNAs. Although fewer circRNAs were obtained with this way, it improved the accuracy rate of our findings. A total of 8843 circRNAs were identified in preadipocytes and adipocytes. This suggested that several circRNAs are expressed in chicken preadipocytes and adipocytes. Of note, the characteristics of circRNAs based on length, exon number, and genomic distribution in chicken were in line with previous reports [19,32]. Previous studies revealed that the expression level of cir-cRNAs was tissues-and stage-specific [33,34]. Abdominal and intramuscular fat deposition resulted from a variety of complex molecular mechanisms. Notably, tissue-specific adipogenesis is potentially regulated by circRNAs [25]. In the present study, we identified DE circRNAs across various tissue and cell types. Also, previous reports identified several DE circRNAs in chicken, such as circCELC19A which was suggested to play functional roles in follicle development [32]. In addition, cir-cFNDC3AL and circTNS1 were identified in embryonic muscle tissues of chicken [19].
Previously, we reported that miRNAs could be crucial regulators in chicken adipogenesis [8,35]. Here, we systematically analyzed and compared the expression patterns of miRNAs during intramuscular and abdominal adipogenesis. It was revealed that several miRNAs were significantly downregulated in mature adipocytes. On the other hand, many genes associated with lipid metabolism and adipogenesis were upregulated during the adipogenic differentiation across various tissues, implying that the downregulated miRNAs potentially promote adipogenesis through their target genes in chicken. Generally, it is believed that circRNAs could function as a molecular sponge of miRNAs thereby regulating their target genes [36]. By integrating the analysis of circRNAs, miRNAs, and mRNAs data, hub ceRNAs networks were constructed during chicken adipogenic differentiation. Notably, the hub genes in the ceRNAs network were enriched in fatty acid metabolism, PPAR signaling pathway, and p53 signaling pathway. Furthermore, the downstream genes in the ceRNA pathways were strongly related to candidate circRNAs in the present study (Fig. 7c and Fig. 8a), suggesting that these circRNAs might play functional roles during adipogenesis. RT-qPCR results of the circFNDC3AL-miRNAs-mRNAs pathway corroborated the results of bioinformatics analysis. Interestingly, fibronectin type III domain-containing protein 3B (FNDC3B) and 5 (FNDC5) regulate white fat browning and adipogenesis [37][38][39][40]. Besides, perlecan (HSPG2, Heparan Sulfate Proteoglycan 2) was found to regulate obesity and lipid deposition [41]. Thus, we speculated that cir-cFNDC3AL and circHSPG2 potentially regulate chicken adipogenesis. However, in-depth studies on the functions of chicken circFNDC3AL and circHSPG2 on adipogenic differentiation are essential. Previous research reported that runtrelated transcription factor 1 (RUNX1T1) [42,43], myosin, heavy chain 9, non-muscle (MYH9), fatty acid desaturase 2 (FADS2), insulin-like growth factor 2 mRNA-binding protein 1 (IGF2BP1) were associated with adipogenesis and fat deposition [44][45][46][47]. We found that circLCLAT1 might influence adipogenic differentiation by regulating downstream genes (RUNX1T1, FADS2, MYH9, IGF2BP3, and PDGFRA) through four miRNAs (miR-34a-5p, miR-30e-5p, miR-146b-5p, and miR-147). Nonetheless, future in-depth studies should be undertaken to explore the molecular regulation mechanism of circFNDC3AL and circLCLAT1/miRNAs/ mRNAs pathway in chicken adipogenesis.

Conclusions
In summary, we comprehensively identified expression profiles of circRNAs in adipocytes derived from various adipose tissues in chicken. By integrating circRNAs, miR-NAs, and mRNAs data, we constructed ceRNAs networks that regulated chicken adipogenesis. Over all, this study offers insights into poultry tissue-specific adipogenesis in poultry and reveals novel clues to studying circRNAs on adipogenic differentiation and fat deposition in poultry. Fig. 7 The construction of circRNAs-miRNAs-mRNAs ceRNAs networks during chicken adipogenesis. a Heatmap of 28 DE miRNAs across various groups. b The pathway enrichment analysis of target genes of miRNAs. c The circRNAs-media ceRNAs networks during chicken adipogenesis. The cube represents circRNAs, circles represent miRNAs and arrows represent genes. Sea-blue color represents a downregulated trend, while dark-red and brown-yellow color represents an upregulated trend. d The pathway enrichment analysis of downstream genes in ceRNAs networks

Tissues and cell samples collection
The experimental birds were obtained from the Avian Farm of Henan Agriculture University (Zhengzhou, Henan, China). Here, five types of tissues (heart, liver, breast muscle, abdominal, and kidney tissues) were collected and stored at − 80°C until use. Under deep anesthesia, the birds were euthanized by intramuscular injection of pentobarbital (Sigma, St. Louis, MO, USA) (40 mg/kg). Intramuscular and abdominal preadipocytes were isolated from the breast muscle and abdominal adipose tissue of 2 weeks old-chicken following the method described by Zhang et al [8,26]. Pollution-free treatment of animal carcasses was performed after the experiment. In vitro adipogenic differentiation model of chicken preadipocytes constructed according to the method by Zhang et al [26]. After 80-90% confluence, cells were exposed to the MDIO differentiation medium. Cell samples were harvested at 0, 2, 4, 6, 8, and 10 days after adipogenic induction.

Identification of differentially expressed (DE) circRNAs and host genes functional annotation
DE circRNAs were analyzed using the edgeR package [52], with the threshold of p-value < 0.05 and |fold change| ≥ 2. Host genes of circRNAs were used for predicting the potential functions of circRNAs. KOBAS 3.0 (http://kobas.cbi.pku.edu.cn/kobas3/?t=1) online software was used for Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses.
Genomic DNA (gDNA) and total RNA extraction, cDNA synthesis and quantitative real-time PCR (qRT-PCR) Genomic DNA (gDNA) and total RNA were extracted from tissues/cell samples using animal genomic DNA kit (Tiangen, China) and RNAiso plus (Takara, Dalian, China) following to the manufacturers' instructions respectively. cDNA was synthesized with a PrimeScriptTM RT reagent Kit with gDNA Eraser (Takara, Dalian, China). Primers used for qRT-PCR were designed by pri-mer3plus (Supplementary Table S3). qRT-PCR was performed with the SYBR® Premix Ex TaqTM II kit (Takara) on a LightCycler® 96 Real-Time PCR system. GAPDH and U6 were used as the reference genes of cir-cRNAs, mRNAs, and miRNAs, respectively. The 2 −ΔΔCt method was used to calculate gene expression levels.

Validation of circRNAs
To validate the circRNAs related to chicken adipogenesis, PCR and Sanger sequencing assays were performed. Divergent primer was designed to amplify the junction sites of circRNA, whereas a convergent primer was designed to amplify the linear mRNA (Supplementary  Table S3). Random primers were used for cDNA synthesis. Then, the cDNA of intramuscular and abdominal preadipocytes and adipocytes were mixed to a cDNA pool. The PCR products were sequenced by Sangon Biotech Co. Ltd. (Shanghai, China).

The integrative analysis of circRNAs, miRNAs, and mRNAs
Since circRNAs can act as sponges of miRNA, the potential miRNAs binding sites in chicken adipogenesis were analyzed. The four groups of circRNA data including IMPre, IMAd, AbPre, and AbAd were analyzed during chicken adipogenesis were analyzed, and the corresponding four groups of miRNA data during adipogenesis were descripted and analyzed in our previous studies [35,53]. miRNAs with |fold change| > 2 and FDR < 0.05 were regarded as differentially expressed miRNAs. TargetScan (http://www.targetscan.org/vert_72/), miRanda (http:// www.microrna.org/microrna/home.do), and RNAhybrid software (https://bibiserv.cebitec.uni-bielefeld.de/rnahybrid) were used to predict the miRNA sites in mRNAs and circRNAs with default parameters, and investigate putative interactions between miRNAs and mRNAs or cir-cRNAs. Furthermore, the target genes with opposite expression trends to miRNAs were regarded as candidate target genes. Based on the co-expression of DE circRNAs and miRNAs (|Pearson's correlation coefficient| > 0.8 and p < 0.05), the circRNAs-miRNAs-mRNAs ceRNA networks were constructed and visualized using Cytoscape (version 3.5.0) (http://www.cytoscape.org/) [54].

Statistical analysis
All statistical data were analyzed using SPSS 22.0 (SPSS Inc., Chicago, IL, USA). All results were presented as mean ± S.E.M, a two-tailed t-test to analyzed the data. p < 0.05 and ** p < 0.01 were considered significant and extremely significant respectively.