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 circRNAs 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|22327655 was significantly downregulated in mature adipocytes when compared with preadipocytes in both IMF and AbF groups. Z:78636651|78640537 (circCDKN2A) was considerably downregulated in AbAd group, while upregulated in the IMAd group. Additionally, 7:22754779|22756470 (circTNS1) was only found in intramuscular fat preadipocytes and adipocytes, 21:6563897|6564555 (circHSPG2), 3:66699169|66707019 (circSLC22A16), and 16:242878|309387 (circKIFC1L) were specifically expressed in intramuscular adipocytes. 3:79378778|79415580 (circBCKDHB), 28:1095851|1096498 and 1:83410785|83437136 were specifically expressed in abdominal adipocytes.
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, circCLEC19A, 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 circRNAs 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 circRNAs data is reliable.
Expression characteristics of miRNAs in chicken intramuscular and abdominal adipogenesis
To comprehensively analyze the characteristics of miRNAs during tissue-specific adipogenesis, miRNAs expression levels were compared across various groups. A total of 80 DE miRNAs (58 upregulated and 22 downregulated) were identified between AbPre and AbAd groups, 225 DE miRNAs (85 downregulated and 140 DE miRNAs upregulated) between AbPre and IMPre groups, 206 DE miRNAs (98 downregulated and 108 upregulated) between AbAd and IMAd groups and 111 DE miRNAs (89 upregulated and 22 upregulated) between IMPre and IMAd groups respectively (Fig. 6A) (Table S2). Among them, 81 (66 up-regulated, 15 down-regulated) and 50 (34 up-regulated, 16 down-regulated) specific-DE miRNAs were identified in IM and Ab groups during adipogenic differentiation respectively (Fig. 6B). Moreover, it was observed that most of DE miRNAs were significantly downregulated in mature adipocytes, including miR-130b-5p, miR-148a-5p, miR-15c-5p, miR-16c-5p, miR-30a-3p (Fold change > 2, Q_value < 0.001). However, miR-146b-5p was significantly upregulated in mature adipocytes compared to preadipocytes (Fold change > 4, Q_value < 0.001). The heatmaps of DE miRNAs across various groups are displayed in Fig. 6C.
Construction of circRNAs-miRNAs-mRNAs Network
Notably, 31 DE miRNAs were shared between IM and Ab groups during adipogenesis, including miR-15c-5p, miR-206, miR-148a-5p, miR-128-1-5p, miR-30a-3p, miR-27b-5p, miR-92-5p and three novel miRNAs (not shown) (Fig. 7A). To further investigate the functional roles of miRNAs in chicken adipogenesis, the target genes of DE miRNAs were predicted through bioinformatics analyses. KEGG pathway analysis revealed that their target genes were highly enriched in the PPAR signaling pathway, metabolic pathways, fatty acid metabolism pathway (Fig. 7B). Further, through integrating analysis of DE cirRNAs, miRNAs and mRNAs data, we constructed a ceRNAs network during chicken adipogenesis (Fig. 7C). It was noted that circFNDC3AL, circHSPG2, circLCLAT1 and circCLEC19A were the hub cirRNAs of the ceRNAs network (Fig. 7C). Furthermore, KEGG pathway analysis showed that the ceRNAs network regulates target genes mainly via fatty acid metabolism, focal adhesion and ECM-receptor interaction pathways (Fig. 7D).
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), 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 circARMH1 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 four 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. 8D). In contrast, the expression levels of circCLEC19A were dramatically decreased after two 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 circLCLAT1 potentially regulates animal adipogenesis or lipid metabolism.
Validation of circRNA-miRNA-mRNAs ceRNA networks
To verify the circRNA-miRNA-mRNAs regulation networks during intramuscular adipogenic differentiation, we analyzed their expression levels in intramuscular preadipocytes and mature adipocytes. Bioinformatics analysis revealed that cirFNDC3AL exhibit many MREs (miRNAs response elements), such as miR-130a-5p, miR-15c-5p, miR-18a/b-3p (Fig. 9A). Results from RNA-Seq and RT-qPCR showed that cirFNDC3AL and its downstream genes were significantly upregulated in differentiated-mature intramuscular adipocytes, compared to intramuscular preadipocytes (Fig. 9C, 9D, 9F). On the other hand, several adipogenic differentiation-related miRNAs were significantly downregulated in adipocytes, compared to intramuscular preadipocytes, including miR-222, miR-456-3p, miR-15C-5p, miR-130b-5p, miR-130a-5p (Fig. 9B, 9E). Additionally, it was worth noting that the target genes of circFNDC3AL were mainly associated with lipid metabolism and adipogenesis, such as PPARG, PLIN2, CPT1A, INSIG1. In addition, bioinformatic analysis and qRT-PCR results found that miR-146b-5p and miR-147 potentially regulate RUNXT1 and FADS2, while miR-34a-5p and miR-30e-5p might regulate PDGFRA, IGF2BP3, FNDC3B, MYH9 (Fig. S2). circLCLAT1 potentially influence chicken adipogenesis by sponging these miRNAs, thus regulating target genes. The findings suggested that circLCLAT1-miRNAs (miR-146b-5p, miR-147, miR-30e-5p, and miR-34a-5p)-mRNAs pathway showed opposite expression patterns during intramuscular adipogenic differentiation (Fig. 10). CircLCLAT1 and its downstream genes were significantly downregulated in differentiated-mature intramuscular adipocytes (Fig. 10C, 10D and 10F), while miR-146b-5p and miR-147 were significantly upregulated in mature adipocytes compared to preadipocytes (Fig. 10B, 10E).