Assessment of CircRNA Expression Profiles and Potential Functions in Brown Adipogenesis

Brown adipose tissue (BAT) is specialized for energy expenditure, thus a better understanding of the regulators influencing BAT development could provide novel strategies to defense obesity. Many protein-coding genes, miRNAs, and lncRNAs have been investigated in BAT development, however, the expression patterns and functions of circRNA in brown adipogenesis have not been reported yet. This study determined the circRNA expression profiles across brown adipogenesis (proliferation, early differentiated, and fully differentiated stages) by RNA-seq. We identified 3,869 circRNAs and 36.9% of them were novel. We found the biogenesis of circRNA was significantly related to linear mRNA transcription, meanwhile, almost 70% of circRNAs were generated by alternative back-splicing. Next, we examined the cell-specific and differentiation stage-specific expression of circRNAs. Compared to white adipocytes, nearly 30% of them were specifically expressed in brown adipocytes. Further, time-series expression analysis showed circRNAs were dynamically expressed, and 117 differential expression circRNAs (DECs) in brown adipogenesis were identified, with 77 upregulated and 40 downregulated. Experimental validation showed the identified circRNAs could be successfully amplified and the expression levels detected by RNA-seq were reliable. For the potential functions of the circRNAs, GO analysis suggested that the decreased circRNAs were enriched in cell proliferation terms, while the increased circRNAs were enriched in development and thermogenic terms. Bioinformatics predictions showed that DECs contained numerous binding sites of functional miRNAs. More interestingly, most of the circRNAs contained multiple binding sites for the same miRNA, indicating that they may facilitate functions by acting as microRNA sponges. Collectively, we characterized the circRNA expression profiles during brown adipogenesis and provide numerous novel circRNAs candidates for future brown adipogenesis regulating studies.

Introduction

In recent years, the number of obese people is increasing rapidly and becoming pandemic. As obesity is highly associated with metabolic syndromes, such as diabetes, cardiovascular diseases, and even cancer, it severely threatens public health (Blüher, 2019). Obesity develops when the energy intake is more than the energy expenditure. Human adipose tissue mainly includes white adipose tissue (WAT) and BAT. The WAT, which contains large lipid droplets, is the main place of energy storage, whereas BAT contains plenty of mitochondria and converts chemical energy into heat to maintain body temperature (Betz and Enerbäck, 2018). It had been considered that BAT only exists in infants, however recent reports detected functional BAT in adult humans (Leitner et al., 2017). The BAT activity positively correlates with human metabolic rate and reduces energy storage, thus strategies that enhance BAT development and increase BAT activity are considered as possible approaches to combat obesity (Bhatt et al., 2017).

Adipogenesis is the process through which preadipocytes differentiate into adipocytes. Previous studies have identified numerous protein regulators of brown adipogenesis, e.g. Pparγ, Prdm16, and Pgc1α (Shapira and Seale, 2019). In recent years, a lot of studies have been focusing on noncoding RNAs. It is known that many miRNAs (Gharanei et al., 2020) and lncRNA (Squillaro et al., 2020) regulate brown adipogenesis and thermogenesis. Recently, circRNA emerged as a new type of noncoding RNA. CircRNAs are covalently linked non-coding RNAs with neither 5′ caps nor 3’ polyadenylated tails, thus it is more stable than linear RNAs (Li et al., 2018a). In 1979, circRNA was first observed in Hela cells by electron microscopy (Hsu and Coca-Prados, 1979). At that time, circRNA was considered as rare byproducts or spliced intermediates. However, with the routine use of high-throughput sequencing technology, circRNAs are expressed in various organisms (Memczak et al., 2015). CircRNAs are highly conserved between species and the expression patterns are tissue-specific and developmental stage-specific. Recent work suggests that circRNAs play vital roles in various biological activities, such as normal tissue development, pathological processes, and even disease progression. CircRNAs exert their function in several ways. First, circRNAs can bind miRNAs and inhibit their functions. Second, circRNAs can act as protein sponges to promote protein interaction or directly affect protein functions. Third, circRNAs can bind to RNA Pol II complex and act as cis-regulators of transcription. Fourth, a few circRNAs are reported to encode functional proteins (Zhang et al., 2019a).

Recent work indicated that circRNAs are involved in adipose tissue development and obesity. We previously reported that 3,771 circRNAs were expressed in white adipogenesis (Zhang et al., 2021a). Arcinas et al. detected thousands of circRNAs in mouse and human WAT. They identified circTshz2-1 and circArhgap5-2 were essential for adipogenesis (Arcinas et al., 2019). Liu et al. found circSAMD4A (hsa_circ_0004846) was significantly upregulated in obese people. CircSAMD4A affected pre-adipocytes differentiation by acting as the miR-138-5p sponge and regulated EZH2 expression (Liu et al., 2020). Wang et al. discovered that 9,311 circRNAs were expressed in duck pre-adipocytes and adipocytes. Circ-PLXNA1 could regulate adipogenesis by binding to miR-214 and affecting CTNNB1expression (Wang et al., 2020). Zhang et al. found circNrxn2 could promote WAT browning by binding to miR-103 and enhance FGF10 expression levels (Zhang et al., 2019b). Recently, circRNAs were reported to be differentially expressed during brown to white adipose tissue transformation in goats (Zhang et al., 2021b). However, the global expression patterns and functions of circRNAs in brown adipogenesis have not been reported.

BAT primary stromal vascular fraction (SVF) cells which contain a lot of primary BAT preadipocytes are widely used to study adipogenesis in vitro (Shan et al., 2016; Oguri et al., 2020). In the present study, we focus on determining the expression profiles and potential roles of circRNAs in brown adipogenesis. The primary SVF cells were isolated from mice interscapular BAT and circRNAs were detected during adipogenesis by RNA-seq. A great number of novel circRNAs were identified and their expression patterns were characterized. Then, we determined the tissue-specific and differentiation stage-specific circRNAs. We also demonstrated that the differential expression circRNAs (DECs) may potentially be involved in regulating brown adipogenesis by acting as miRNA sponges.

Materials and Methods

BAT SVF Cells Culture

Primary BAT SVF cells were isolated from C57BL/6J background mice as described (Shan et al., 2016). Briefly, the interscapular brown fat was dissected from 8-week-old mice. Next, the brown fat was cut to fine pieces and digested with 1.5 mg/ml collagenase type I (catalog number SCR103, Sigma-Aldrich) at 37°C water bath for 50 min. The tissue was vortexed every 10 min. Then, the digestion was filtered through 100-μm and 70-μm cell strainers. The digestion was centrifuged at 400 g for 6 min to enrich SVF cells. The cells were seeded into a cell culture plate with growth medium (DMEM and 10% fetal bovine serum). To induce adipogenic differentiation, the growth medium was supplemented with 2.85 mM recombinant human insulin (catalog number I8830, Solarbio), 0.3 mM dexamethasone (catalog number D8040, Solarbio), and 0.63 mM 3-isobutyl-methylxanthine (catalog number I7018, Sigma-Aldrich). After 96 h, the medium was switched to growth medium, supplemented with 10 nM triiodothyronine (catalog number T6397, Sigma-Aldrich), and 200 nM insulin to induce mature adipocytes. For Oil red O staining, cells were fixed with 10% formaldehyde for 5 min and stained with staining solutions (catalog number G1262, Solarbio) according to the instructions.

RNA-Sequencing

Total RNA of the BAT adipocytes was collected using Trizol (catalog number 15596026, Thermo Fisher Scientific). The total RNA was treated with Epicentre Ribozero rRNA Removal kit (catalog number RZH1046, Epicentre) and RNase R (catalog number RNR07250, Epicentre) to remove ribosomal RNA and linear RNA. Then, the sequencing libraries were prepared by NEBNext Ultra Directional RNA Library Prep Kit (catalog number E7760S, NEB, USA) following the manufacturer’s instructions. At last, the libraries were sequenced on an illumine platform and 150 bp paired reads were generated.

Computational Analyses of CircRNA

CircRNAs were detected by CIRI2 as previously described (Zhang et al., 2018). Then the sequence data were mapped to mm9 with BWA-MEM algorithm (Li, 2013). CircRNA candidates mapped by at least two reads in both replicated were kept for subsequent study. Then, the DECs were identified with the likelihood ratio test from R package DESeq2 (version 1.10.1) (Hiraike et al., 2017). The DECs were clustered based on their expression patterns by the degPatterns function from the DEGreport package (Pantano, 2021).

RT-PCR

The cDNA was obtained by using random primers and a reverse transcription kit (catalog number RR037A, Takara). CircRNA sequences were amplified by divergent primers which were designed by Primer3 (https://primer3.ut.ee/) to cover the back-splicing sites. PCR products were extracted by agarose gel DNA extraction kit (catalog number 9762, Takara). Then the purified DNA was sent to perform Sanger sequencing. Then qPCR was performed to examine the expression levels of circRNAs on a LightCycler 96 system (Roche, Germany) using TB Green Premix Ex II (catalog number RR820A, Takara) according to the instructions. The relative expression levels were normalized to 18S.

CircRNA Expression Patterns Analysis

The expression patterns of circRNAs were classified by the STEM program with No normalization/add 0 option (Ernst and Bar-Joseph, 2006). GO annotations were obtained from the STEM program with the default option.

CircRNA-miRNA Networks Construction

The circRNA sequences were obtained from package circPrimer (Zhong et al., 2018)and the potential miRNA binding sites were predicted by miRanda (version 3.3a) (John et al., 2004). To obtain more rigorous prediction results, the energy threshold was set to −7 kcal/mol and the score threshold was set to 150, while the other parameters were set to default values. Then, the circRNA-miRNA interaction networks were constructed by Cytoscape (version 3.8.2) (Shannon et al., 2003).

Results

Identification of circRNAs in brown Adipogenesis

To detect circRNAs transcripts in brown adipogenesis, we cultured primary BAT SVF cells and induced them to differentiate. As shown in Figure 1A, the BAT SVF cells were in fibroblast-like spindle shape on day 0 (D0, proliferation stage), whereas the cells round-up on day 4 (D4, early differentiated stage) and were filled with large lipid droplets on day 8 (D8, fully differentiated stage) as indicated by Oil red O staining. Total RNA was collected on D0, D4, and D8 post differentiation. The circRNAs were enriched and then RNA-seq was performed. CIRI2, which is based on back splicing alignment, was used to detect circRNAs. Then the high confident circRNAs with at least 2 reads were detected in both 2 replicates were selected for subsequent analysis. In this way, a total of 3,869 distinct circRNA candidates were detected, indicating that a larger number of circRNAs were expressed in brown adipogenesis (Supplementary Table S1). Notably, up to 36.9% of the circRNAs were novel and not annotated in the circBase database (http://www.circbase.org/). We found that only 32.3% of circRNAs continue to be expressed, while 547, 435, and 696 circRNAs were specifically expressed in proliferation, early differentiated, and fully differentiated stages respectively (Figure 1B). We found that 48 circRNAs were not aligned to the annotated genome sequence, the other 3,821 circRNAs were derived from 2,046 gene locus. Based on our previous study, 20,968 genes were transcribed to mRNA during brown adipogenesis (with a minimum of 4 reads, GEO accession number GSE173710). Thus, as many as 9.61% of genes could produce circRNAs.

FIGURE 1

FIGURE 1. Identification and characterization of circRNAs in brown adipogenesis. (A) Representative pictures of BAT pre-adipocytes differentiation. Scale bar = 20 µm. (B) The number of CircRNAs identified during BAT pre-adipocytes differentiation. (C) The number of circRNA based on back-splicing regions in each adipogenesis stage. (D) The number of circRNAs derived from each chromosome. (E) Correlation of circRNA and mRNA number from the same chromosome. (F) Principal Component Analysis (PCA) of circRNA expression during BAT and WAT pre-adipocytes differentiation. (G) Venn diagrams of circRNAs expressed in BAT and WAT adipogenesis.

According to the back-splicing regions, circRNAs are classified into exonic, intronic, and intergenic circRNAs. We found the ratio of the three types of circRNAs were similar across brown adipogenesis, in which exonic circRNAs were identified as the major type which accounts for more than 93% of all circRNAs, then followed by intronic circRNAs and intergenic circRNAs (Figure 1C). Next, we examined the chromosome distribution of circRNAs. It is noticed that circRNAs were unevenly distributed. Chromosome 2 gave rise to 311 circRNAs which is the greatest, while chromosome X produced only 78 circRNAs (Figure 1D). We suspected that the number of circRNA was related to the length of the chromosome. Thus, we calculated the association between them, and a significant correlation was found (R = 0.51, p < 0.05). It is reported that both circRNAs and liner mRNAs are sliced from pre-RNAs, we further examined the association between circRNA and linear mRNA number from the same chromosome (R = 0.69, p < 0.001, Figure1E). These results suggested that the biogenesis of circRNA was related to linear mRNA transcription.

As mentioned earlier, BAT and WAT are functionally distinct adipose tissue. To compare the circRNA expression profiles, we collected the WAT SVF circRNA expression data from our previous study (GEO accession number GSE178502), in which the same mice were used as the current study and WAT SVF were isolated from inguinal WAT. Principal component analysis (PCA) analysis showed that samples derived from different tissues located far away, indicating BAT and WAT circRNA expression profiles were largely different (Figure 1F). Then we examined the tissue-specific circRNAs. As shown in (Figure 1G), 2,725 circRNAs were co-expressed in both adipose tissues, while 1,144 and 986 circRNAs were specifically expressed in BAT and WAT, respectively. Thus, almost 30% of the circRNAs were specific to BAT or WAT.

Alternative Back-Splicing of CircRNAs

In most cases, a host gene only produced one circRNA, however, these circRNAs only account for 31.8% of all the circRNAs expressed in BAT. The other circRNAs came from alternative back splicing (Figure 2A). Strikingly, some host genes even produced more than 10 circRNA isoforms, for example, Atrophin 2, which acts as transcriptional co-repressors and plays key roles during embryogenesis (Zoltewicz et al., 2004), produced up to 20 distinct circRNAs. The above results suggest that alternative back-splicing is an important source of circRNA diversity.

FIGURE 2

FIGURE 2. Alternative back-splicing of circular RNAs. (A) Distribution of circRNAs variants among genes. (B) Generation model of circRNA variants derived from the same acceptor in Nedd4. Splicing acceptor is indicated with ▼, while splicing donor is indicated with▽. (C) Generation model of circRNA variants derived from same acceptor or same donor in Nifa. (D–F) Expression profile of circRNA isoforms for Nedd4, Glis3, and Nifa across brown adipogenesis.

As reported, one gene can make several circRNAs by selectively using different splice donors or splice acceptors (Liang et al., 2017). For example, Nedd4 is an E3 ubiquitin ligase that is involved in adipogenesis by regulating Pparγ stability (Yao et al., 2019). Nedd4 could produce three distinct circRNA isoforms by using the same splice acceptor and different splice donors (Figure 2B). Another example is Glis3, which is a transcription factor and plays a key role in neonatal diabetes, type 1 and type 2 diabetes (Shan et al., 2016). Glis3 produced three circRNA isoforms (circ_0000942, circ_0007671, and circ_0007670) that also share the same acceptor. Another example is Nifa, which co-localizes with Pparγ and transcriptionally induces the brown fat gene expression during brown adipocyte differentiation (Hiraike et al., 2017). As shown in Figure 2C, Nifa produced four circRNA isoforms (circ_Nifa_1, circ_Nifa_2, circ_Nifa_3, and circ_0011879). The alternative splicing was very complex. circNifa1 and circNifa3 shared an acceptor, while circ_0011879 and circNifa2 shared another acceptor. Meanwhile, circ_0011879 and circNifa1 shared a donor.

Next, we explored the expression of these circRNA isoforms. For circRNAs derived from Nedd4, we found circNedd4_1 and circNedd4_2 were the dominant isoforms and the expression patterns of the three circRNA isoforms were completely different. CircNedd4_2 continually increased, while circNedd4_2 decreased on D4. CircNedd4_3 began to increase on D8 (Figure 2D). For the circRNAs derived from Glis3, circ_0000942 increased from D0 to D4, and then start to decrease. Circ_0007671 and circ_0007671 showed similar expression profiles, which kept unchanged until D4 and then decreased to the lowest expression level on D8 (Figure 2E). For the circRNAs derived from Nifa, circ_Nifa_1 and circ_0011879 were the dominant isoforms. Then expression patterns of them were similar, which showed an increase on D4 and then decreased (Figure 2F). In summary, alternative back-splicing expands the diversity of circRNAs, and expression among the back-splicing variants is different.

Expression Patterns of CircRNAs in BAT Adipogenesis

As mentioned earlier, BAT preadipocytes undergo a tremendous change of morphology during adipogenesis. To examine the overall circRNA expression levels, we did a boxplot and found that the levels of circRNAs during brown adipogenesis were comparable to each other (Figure 3A). To analyze the variation of circRNA expression profiles, we performed PCA. As shown in Figure 3B, the samples derived from the same group were located near each other, indicating the high repeatability of the results. As expected, the samples came from D4 and D8 located nearer than the D0, indicating the circRNA expression patterns were similar in differentiation stages, whereas they were largely different with the proliferation stage. Hierarchical clustering heatmap showed the correlation between different samples ranged from 0.88 to 0.98 (Figure 3C). Replicates in each stage were highly correlated, meanwhile, the D4 and D8 samples were clustered together, which was consistent with the PCA results.

FIGURE 3

FIGURE 3. The expression patterns of circRNA during brown adipogenesis. (A) Relative expression levels of circRNAs during brown adipogenesis. (B) PCA plot of circRNA expression during brown adipogenesis. (C) Correlation of circRNA expression calculated by DESeq2 rlog-normalized RNA-seq data. (D) Time series analysis of circRNA expression patterns. The number on the top left indicated the expression profile number. (E–H) GO terms enriched for the host genes of indicated circRNA clusters.

To get an overview of circRNA expression patterns in brown adipogenesis, we carried out a time series analysis by STEM. 14 clusters were significantly enriched (with colored background, Figure 3D), which belong to 6 groups (with the same color). When we examined the GO enrichment results, we found green and yellow clusters were interesting, which correspond to the up and down regulated circRNAs. GO analysis for their parental gene showed the green cluster was associated with brown adipogenesis (Supplementary Table S2). For example, most genes in cluster 40 were associated with lipid metabolisms, such as sterol metabolic process, cholesterol metabolic process, and lipid catabolic process (Figure 3E). The genes in cluster 42 were enriched with GO terms of cold-induced thermogenesis, adaptive thermogenesis, and regulation of fat cell differentiation (Figure 3F). While in the yellow group (clusters 35 and 44), most genes were associated with cell mitosis. For instance, regulation of cell division, mitotic cell cycle process, and mitotic sister chromatid segregation (Figures 3G,H). These results were consistent with the process of brown adipogenesis, in which pre-adipocytes exit the cell cycle, enter the differentiation process, and obtain the function of thermogenesis.

Differentially Expressed CircRNAs in brown Adipogenesis

To examine the DECs in the time-course of brown adipogenesis, we performed the Likelihood ratio test by DESeq2, and the cut-off was set as adjusted p < 0.05. As a result, 117 DECs were identified, with 77 upregulated and 40 downregulated (Supplementary Table S3). Notably, much more up-regulated circRNAs were found in the D8 group than in the D0 and D4 groups (Figure 4A). The volcano plot showed log 2 fold change ranged from −7.5 to 10 (Figure 4B). The top upregulated circRNA was circ_Acss3_1, which could not be detected on D0, then it increased to 479 circCPM on D4 and reached an extremely high level on D8 (2428 circCPM).

FIGURE 4

FIGURE 4. Dynamic differentially expressed circRNAs in brown adipogenesis. (A) Heatmap of DECs across different stages of brown adipogenesis. (B) Volcano plot showing circRNA abundance across brown adipogenesis. The green color presents the significant changed DECs (padj <0.05). (C) Expression patterns of the DECs analyzed by degPatterns function from the DEGreport package (minc = 5).

It has been suggested that circRNA may positively or negatively regulate host gene expression (Shao et al., 2021). To explore the relationship between DECs and their host genes, we collected mRNA expression profiles in brown adipogenesis from our previous study (GEO accession number GSE173710). Then we calculated the correlation of them by Pearson correlation test. We found 40 circRNA-mRNA pairs were significantly correlated (p < 0.05, Supplementary Table S4). Interestingly, all of them showed positive correlations, which are consistent with our previous study in white adipogenesis (Zhang et al., 2021a). These results suggested that these DECs may potentially regulate their host genes expression in brown adipogenesis.

It is reported that some of the circRNAs showed highly evolutionary conservation between humans and mice (Jeck et al., 2013). To evaluate the conserved circRNAs, we obtained the human circRNA sequences from circBase (Glažar et al., 2014). Then, alignments were conducted using the NCBI-BLAST-2.11.0+ (https://www.ncbi.nlm.nih.gov/books/NBK131777/) to identify the regions of the mouse circRNAs that corresponded to the human circRNAs. Among the 117 mouse DECs, 85 of them aligned to human circRNAs (E value <10⁻⁵, Supplementary Table S5). In addition, it should be noted that the number of the conserved circRNAs may be underestimated, as circRNAs are tissue-specific and their expression profile in human BAT is not reported yet. These results indicated that a high ratio of the DECs is conserved between mice and humans.

According to the previous study, the circSAMD4A (derived from host gene SAMD4A, also known as SAMD4) regulated white pre-adipocytes differentiation (Liu et al., 2020). In the current study, we also found the circRNA was differentially expressed in brown adipogenesis, suggesting that it may also be involved in brown adipogenesis. Another reported circRNA derived from Arhgap5 (Arcinas et al., 2019) was also detected in our study, but it was not significantly differentially expressed in brown adipogenesis.

Then, we classified the genes that exhibited similar change across adipogenesis by using DEGreport which uses a hierarchical clustering approach based on pair-wise correlations. The DECs were classified into four groups (Figure 4C; Supplementary Table S3). The genes included in each group ranged from 7 to 63. The largest one was group1, which contained 63 circRNAs with a continuously increasing trend during adipogenesis. Then, it was group 3, which contained 33 circRNAs and kept on decreasing in adipogenesis. Group 2 contained 14 circRNAs which showed a sharp increase after D4. Group 4 contained only 7 circRNAs which were found to increase on D4 and then declined on D8. Although GO analysis showed no terms were significantly enriched, we identified many DECs parental genes which are essential in adipogenesis and brown adipocyte bioenergetics, such as Pparγ, Zbtb16, Sik2, Snrk, Nfia, and Tead1.

Validation of circRNAs

To determine the authenticity of the RNA-seq results, seven of the DECs were selected randomly to detect the back-splicing sites by Sanger sequencing. The circRNA sequence and genome locus were obtained by circPrimer (version 2.0). Then, we performed RT-PCR using divergent primers (version 4.1.0, Supplementary Table S6) and all the circRNAs were successfully amplified. We collected the corresponding DNA products and performed Sanger sequencing, as shown in Figures 5A,B, the expected back-splicing sites were detected. Then, the relative circRNAs expression levels were detected by qPCR. Compared to RNA-seq results, the expression trend of most circRNAs was similar. (Figure 5C). Further, we found the qPCR and RNA-seq results were significantly correlated (R = 0.580, p < 0.01, Figure 5D).

FIGURE 5

FIGURE 5. Verification of circRNAs. Representative circRNA generation models and Sanger sequencing results of circRNAs (A) circ_0000529, (B) circ_0001447. (C) Relative circRNA expression levels were detected by qPCR. (D) Correlation analysis of RNA-seq results and qPCR results (n = 3).

Prediction of microRNA-circRNA Interaction

Recently, it has been reported that circRNAs can act as miRNA sponges and repress miRNA functions (Hansen et al., 2013). As miRNAs control a large set of biological processes in adipogenesis, circRNAs may be involved in the process through miRNAs (Zhong et al., 2019). To validate the hypothesis, we chose the top 10 abundant DECs to predict the potential circRNA-miRNA interactions with miRanda. We found all of them were predicted to contain miRNA binding sites, except circ_0001511. To reduce false positives, we only kept the miRNAs that could be detected in brown adipogenesis (GEO accession number GSE45499) (Chen et al., 2013). In the end, we predicted 49 miRNAs may interact with the DECs. Cytoscape was used to construct the miRNA-circRNA interaction networks (Figure 6; Supplementary Table S7). Notably, each circRNA contained several different miRNA binding sites and a total of 97 miRNA binding sites were identified. More interestingly, most of the circRNAs contained multiple binding sites for the same miRNA, for example, circ_001396 was predicted to contain as many as 8 binding sites for mir-466f-3p, which greatly increased the possibility of their interaction.

FIGURE 6

FIGURE 6. The predicted circRNA-miRNA interaction networks. The potential miRNA binding sites of circRNAs were predicted by miRanda. The edge number represents the number of miRNA binding sites in the indicated circRNA sequences. The miRNAs with brown color indicate that these miRNAs are associated with adipogenesis or metabolic diseases according to HMDD (version 3.2) or previous study.

The Human microRNA Disease Database (HMDD) collects experiment-supported miRNAs associated with disease (Huang et al., 2019). We further explored the potential biological functions of the identified miRNAs. According to HMDD (version 3.2) and published data, 20 of the 49 miRNAs were revealed to be associated with adipogenesis or metabolic diseases, including obesity, mitochondrial metabolism disease, type 2 diabetes, and non-alcoholic fatty liver disease (Figure 6; Supplementary Table S8). For example, let-7, miR-337, miR-503, and miR-182 were identified in the networks. They were reported to play regulatory roles in adipogenesis or thermogenesis. In addition, miR-15a and miR-30e may be involved in mitochondrial metabolism disease and obesity, respectively. In summary, the DECs may potentially interact with these miRNAs to exert their regulatory functions.

Discussion

Previous studies had shown that a lot of circRNAs were expressed in adipose tissue. During Yak adipocytes differentiation, 7,203 circRNA were detected (Zhang et al., 2020). In human WAT, 6,925 circRNAs were detected, while in mouse WAT 2,380 circRNAs were detected (Arcinas et al., 2019). We previously showed 3,771 circRNA were detected during mouse white adipogenesis. However, the circRNA expression in BAT has not been reported. In the present work, the expression profiles of circRNAs were investigated during brown adipogenesis. As many as 3,671circRNAs were identified and 36.9% of them were not annotated yet. This is the first study that reports the expression of circRNAs expression profiles in brown adipogenesis. Compare to our previous work studying circRNA in white adipogenesis, we found 1,144 of the circRNAs were specifically expressed in BAT. Further, we showed many circRNA were differentialy expressed during brown adipogenesis. These results are consistent with the idea that circRNAs are tissue-specific and developmental stage specific.

The functions of many circRNAs remain unclear at present. However, accumulated reports suggest they can regulate parental gene expression through various mechanisms (Shao et al., 2021). As both circRNAs and linear RNAs are generated from pre-RNA, they compete with each other for splicing, thus decrease the expression levels of parental genes (Ashwal-Fluss et al., 2014). On the other hand, circRNAs can increase parental gene expression. For example, circEIF3J and circPAIP2 bind to RNA Pol II and U1 snRNA, then activate parental genes transcription. Circ-Sirt1 can compete with miR-132/212 to bind Sirt1, leading to the enhancement of Sirt1 expression (Kong et al., 2019). CircRNA sisR-4 promoted parental gene expression by activating enhancer (Tay and Pek, 2017). Based on the idea that circRNAs may regulate parental genes expression, we predicted their potential functions. According to the annotation of the parental genes, downregulated circRNAs were enriched in GO terms related to cell proliferation and cell cycle, whereas the up-regulated circRNAs were enriched in cell differentiation and thermogenesis. These results were consistent with the activity of brown adipogenesis. At the early stage of brown adipogenesis, preadipocytes undergo a post confluent mitosis and exit the cell cycle. Then the pre-adipocytes are committed to adipocytes. A batch of adipose-related genes began to accumulate (Ntambi and Young-Cheul, 2000). In addition, we found many parental genes of the DECs participated in adipogenesis modulation. For example, Pparγ is indispensable for adipogenesis (Brun et al., 1996), which produced five circRNA isoforms and one of them showed a continuously increasing expression. Circ_0001335 derived from Nsd2 which inhibits H3K27me3 and increases expression of C/EBPα and Pparγ, thus promoting adipogenesis (Brun et al., 1996). Fkbp5 and Fndc3b are essential for adipogenesis (Muraoka et al., 2009; Li et al., 2018b), both can also produce circRNAs. Meanwhile, some of the DECs parental genes are involved in regulating brown adipocyte bioenergetics. Parental gene Zbtb16 is a transcription factor that is induced upon cold exposure. It can increase the expression of BAT marker genes and β-oxidation genes (Plaisier et al., 2012). Another circRNA parental gene was Sik2. Sik2/TORC2 signaling cascade regulates PGC-1α and UCP-1 gene expression in BAT (Muraoka et al., 2009). In addition, several other circRNA parental genes, such as Snrk (Li et al., 2018b), Nfia (Hiraike et al., 2017), and Tead1 (Tharp et al., 2018) are also involved in BAT-specific genes expression and bioenergetics. In summary, we find many parental genes are key adipogenesis regulators.

It is reported that circRNA can post-transcriptional modulate gene expression by binding miRNA. Several circRNAs have been reported to affect adipogenesis. In obese individuals, CircSAMD4A expression increased. It could bind with miR-138-5p and induce adipogenesis (Yang et al., 2011; Liu et al., 2020). In bone marrow mesenchymal stem cells, CDR1as increases adipogenesis by binding with miR-7-5p (Chen et al., 2020). In the current study, we identified that 9 of the top 10 abundant DCEs contained miRNA binding sites. We found circ_0001017 was predicted to contain several Let-7 binding sites. Let-7 is one of the well-studied miRNAs which inhibits pre-adipocytes clonal expansion and terminal differentiation via targeting HMGA2 (Sun et al., 2009). MiR-337-3p can inhibit TWIST1 and promote transcription of genes participating in brown fat metabolism (Vonhögen et al., 2020). In circ_0000819, two miR-337-3p binding sites were identified. In addition, we found that circ_0001017 may interact with miR-503 which was reported to regulate brown adipogenesis by inhibiting BMPR1a (Man et al., 2020). Thus, these DECs may be involved in brown adipogenesis through post-transcriptional regulation by binding with miRNAs.

Conclusion

Collectively, we characterized circRNA expression patterns during brown adipogenesis. We identified a large number of novel circRNAs and found alternative back splicing is an important source of circRNA diversity. We found many of the circRNAs are brown adipogenesis and differentiation stage specific. We also predicted that adipose circRNAs may regulate adipogenesis by acting as microRNA sponges. These novel circRNAs may potentially serve as new regulators of BAT development.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The animal study was reviewed and approved by the Xinyang Normal University Animal Care and Use Committee.

Author Contributions

PZ designed the experiments and wrote the manuscript, PZ, MS, CD, and XC performed the experiments, PZ, ZC, YX, CL, and HX analyzed the data. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (31601167, 31972537), Major Science and Technology Projects of Henan Province (201300111200), the Central Plains Technological Innovation Leading Talents Project of Henan Province (194200510022), Hainan Provincial Natural Science Foundation of China (320CXTD445), Department of Education in Henan Province (21A230017), the Nanhu Scholars Program of XYNU.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Lei Wang, Feng Xing for their help in data analysis.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.769690/full#supplementary-material

References

Arcinas, C., Tan, W., Fang, W., Desai, T. P., Teh, D. C. S., Degirmenci, U., et al. (2019). Adipose Circular RNAs Exhibit Dynamic Regulation in Obesity and Functional Role in Adipogenesis. Nat. Metab. 1, 688–703. doi:10.1038/s42255-019-0078-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashwal-Fluss, R., Meyer, M., Pamudurti, N. R., Ivanov, A., Bartok, O., Hanan, M., et al. (2014). circRNA Biogenesis Competes with Pre-mRNA Splicing. Mol. Cel 56, 55–66. doi:10.1016/j.molcel.2014.08.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Betz, M. J., and Enerbäck, S. (2018). Targeting Thermogenesis in Brown Fat and Muscle to Treat Obesity and Metabolic Disease. Nat. Rev. Endocrinol. 14, 77–87. doi:10.1038/nrendo.2017.132

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhatt, P. S., Dhillo, W. S., and Salem, V. (2017). Human Brown Adipose Tissue - Function and Therapeutic Potential in Metabolic Disease. Curr. Opin. Pharmacol. 37, 1–9. doi:10.1016/j.coph.2017.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Blüher, M. (2019). Obesity: Global Epidemiology and Pathogenesis. Nat. Rev. Endocrinol. 15, 288–298. doi:10.1038/s41574-019-0176-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Brun, R. P., Tontonoz, P., Forman, B. M., Ellis, R., Chen, J., Evans, R. M., et al. (1996). Differential Activation of Adipogenesis by Multiple PPAR Isoforms. Genes Development 10, 974–984. doi:10.1101/gad.10.8.974

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, G., Wang, Q., Li, Z., Yang, Q., Liu, Y., Du, Z., et al. (2020). Circular RNA CDR1as Promotes Adipogenic and Suppresses Osteogenic Differentiation of BMSCs in Steroid-Induced Osteonecrosis of the Femoral Head. Bone 133, 115258. doi:10.1016/j.bone.2020.115258

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Siegel, F., Kipschull, S., Haas, B., Fröhlich, H., Meister, G., et al. (2013). miR-155 Regulates Differentiation of Brown and Beige Adipocytes via a Bistable Circuit. Nat. Commun. 4, 1769. doi:10.1038/ncomms2742

PubMed Abstract | CrossRef Full Text | Google Scholar

Ernst, J., and Bar-Joseph, Z. (2006). STEM: a Tool for the Analysis of Short Time Series Gene Expression Data. BMC Bioinformatics 7, 191. doi:10.1186/1471-2105-7-191

PubMed Abstract | CrossRef Full Text | Google Scholar

Gharanei, S., Shabir, K., Brown, J. E., Weickert, M. O., Barber, T. M., Kyrou, I., et al. (2020). Regulatory microRNAs in Brown, Brite and white Adipose Tissue. Cells 9, 2489. doi:10.3390/cells9112489

PubMed Abstract | CrossRef Full Text | Google Scholar

Glažar, P., Papavasileiou, P., and Rajewsky, N. (2014). circBase: a Database for Circular RNAs. RNA 20, 1666–1670. doi:10.1261/rna.043687.113

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansen, T. B., Jensen, T. I., Clausen, B. H., Bramsen, J. B., Finsen, B., Damgaard, C. K., et al. (2013). Natural RNA Circles Function as Efficient microRNA Sponges. Nature 495, 384–388. doi:10.1038/nature11993

PubMed Abstract | CrossRef Full Text | Google Scholar

Hiraike, Y., Waki, H., Yu, J., Nakamura, M., Miyake, K., Nagano, G., et al. (2017). NFIA Co-localizes with PPARγ and Transcriptionally Controls the Brown Fat Gene Program. Nat. Cel Biol. 19, 1081–1092. doi:10.1038/ncb3590

PubMed Abstract | CrossRef Full Text | Google Scholar

Hsu, M.-T., and Coca-Prados, M. (1979). Electron Microscopic Evidence for the Circular Form of RNA in the Cytoplasm of Eukaryotic Cells. Nature 280, 339–340. doi:10.1038/280339a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, Z., Shi, J., Gao, Y., Cui, C., Zhang, S., Li, J., et al. (2019). HMDD v3.0: a Database for Experimentally Supported Human microRNA-Disease Associations. Nucleic Acids Res. 47, D1013–D1017. doi:10.1093/nar/gky1010

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeck, W. R., Sorrentino, J. A., Wang, K., Slevin, M. K., Burd, C. E., Liu, J., et al. (2013). Circular RNAs Are Abundant, Conserved, and Associated with ALU Repeats. RNA 19, 141–157. doi:10.1261/rna.035667.112

PubMed Abstract | CrossRef Full Text | Google Scholar

John, B., Enright, A. J., Aravin, A., Tuschl, T., Sander, C., Marks, D. S., et al. (2004). Human microRNA Targets. Plos Biol. 2, e363. doi:10.1371/journal.pbio.0020363

PubMed Abstract | CrossRef Full Text | Google Scholar

Kong, P., Yu, Y., Wang, L., Dou, Y.-Q., Zhang, X.-H., Cui, Y., et al. (2019). circ-Sirt1 Controls NF-Κb Activation via Sequence-specific Interaction and Enhancement of SIRT1 Expression by Binding to miR-132/212 in Vascular Smooth Muscle Cells. Nucleic Acids Res. 47, 3580–3593. doi:10.1093/nar/gkz141

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H. (2013). Aligning Sequence Reads, Clone Sequences and Assembly Contigs with BWA-MEM, arXiv preprint arXiv:1303.3997.

Google Scholar

Leitner, B. P., Huang, S., Brychta, R. J., Duckworth, C. J., Baskin, A. S., McGehee, S., et al. (2017). Mapping of Human Brown Adipose Tissue in Lean and Obese Young Men. Proc. Natl. Acad. Sci. USA 114, 8649–8654. doi:10.1073/pnas.1705287114

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Feng, B., Nie, Y., Jiao, P., Lin, X., Huang, M., et al. (2018). Sucrose Nonfermenting-Related Kinase Regulates Both Adipose Inflammation and Energy Homeostasis in Mice and Humans. Diabetes 67, 400–411. doi:10.2337/db17-0745

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Yang, L., and Chen, L.-L. (2018). The Biogenesis, Functions, and Challenges of Circular RNAs. Mol. Cel 71, 428–442. doi:10.1016/j.molcel.2018.06.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, G., Yang, Y., Niu, G., Tang, Z., and Li, K. (2017). Genome-wide Profiling of Sus scrofa Circular RNAs across Nine Organs and Three Developmental Stages. DNA Res. 24, 523–535. doi:10.1093/dnares/dsx022

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Liu, H., Li, Y., Mao, R., Yang, H., Zhang, Y., et al. (2020). Circular RNA SAMD4A Controls Adipogenesis in Obesity through the miR-138-5p/EZH2 axis. Theranostics 10, 4705–4719. doi:10.7150/thno.42417

PubMed Abstract | CrossRef Full Text | Google Scholar

Man, X. F., Hu, N., Tan, S. W., Tang, H. N., Guo, Y., Tang, C. Y., et al. (2020). Insulin Receptor Substrate‐1 Inhibits High‐fat Diet‐induced Obesity by browning of white Adipose Tissue through miR‐503. FASEB j. 34, 12308–12323. doi:10.1096/fj.201903283rr

PubMed Abstract | CrossRef Full Text | Google Scholar

Memczak, S., Papavasileiou, P., Peters, O., and Rajewsky, N. (2015). Identification and Characterization of Circular RNAs as a New Class of Putative Biomarkers in Human Blood. PLoS ONE 10, e0141214. doi:10.1371/journal.pone.0141214

PubMed Abstract | CrossRef Full Text | Google Scholar

Muraoka, M., Fukushima, A., Viengchareun, S., Lombès, M., Kishi, F., Miyauchi, A., et al. (2009). Involvement of SIK2/TORC2 Signaling cascade in the Regulation of Insulin-Induced PGC-1α and UCP-1 Gene Expression in Brown Adipocytes. Am. J. Physiology-Endocrinology Metab. 296, E1430–E1439. doi:10.1152/ajpendo.00024.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Ntambi, J. M., and Young-Cheul, K. (2000). Adipocyte Differentiation and Gene Expression. J. Nutr. 130, 3122S–3126S. doi:10.1093/jn/130.12.3122s

CrossRef Full Text | Google Scholar

Oguri, Y., Shinoda, K., Kim, H., Alba, D. L., Bolus, W. R., Wang, Q., et al. (2020). CD81 Controls Beige Fat Progenitor Cell Growth and Energy Balance via FAK Signaling. Cell 182, 563–577. e20. doi:10.1016/j.cell.2020.06.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Pantano, L. (2021). DEGreport: Report of DEG Analysis, R package version 1.28.0.

Google Scholar

Plaisier, C. L., Bennett, B. J., He, A., Guan, B., Lusis, A. J., Reue, K., et al. (2012). Zbtb16 Has a Role in Brown Adipocyte Bioenergetics. Nutr. Diabetes 2, e46. doi:10.1038/nutd.2012.21

PubMed Abstract | CrossRef Full Text | Google Scholar

Shan, T., Zhang, P., Jiang, Q., Xiong, Y., Wang, Y., and Kuang, S. (2016). Adipocyte-specific Deletion of mTOR Inhibits Adipose Tissue Development and Causes Insulin Resistance in Mice. Diabetologia 59, 1995–2004. doi:10.1007/s00125-016-4006-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res. 13, 2498–2504. doi:10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao, T., Pan, Y.-h., and Xiong, X.-d. (2021). Circular RNA: an Important Player with Multiple Facets to Regulate its Parental Gene Expression. Mol. Ther. - Nucleic Acids 23, 369–376. doi:10.1016/j.omtn.2020.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Shapira, S. N., and Seale, P. (2019). Transcriptional Control of Brown and Beige Fat Development and Function. Obesity 27, 13–21. doi:10.1002/oby.22334

PubMed Abstract | CrossRef Full Text | Google Scholar

Squillaro, T., Peluso, G., Galderisi, U., and Di Bernardo, G. (2020). Long Non-coding RNAs in Regulation of Adipogenesis and Adipose Tissue Function. Elife 9, e59053. doi:10.7554/eLife.59053

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, T., Fu, M., Bookout, A. L., Kliewer, S. A., and Mangelsdorf, D. J. (2009). MicroRNA Let-7 Regulates 3T3-L1 Adipogenesis. Mol. Endocrinol. 23, 925–931. doi:10.1210/me.2008-0298

PubMed Abstract | CrossRef Full Text | Google Scholar

Tay, M. L.-I., and Pek, J. W. (2017). Maternally Inherited Stable Intronic Sequence RNA Triggers a Self-Reinforcing Feedback Loop during Development. Curr. Biol. 27, 1062–1067. doi:10.1016/j.cub.2017.02.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Tharp, K. M., Kang, M. S., Timblin, G. A., Dempersmier, J., Dempsey, G. E., Zushin, P.-J. H., et al. (2018). Actomyosin-Mediated Tension Orchestrates Uncoupled Respiration in Adipose Tissues. Cel Metab. 27, 602–615. e4. doi:10.1016/j.cmet.2018.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Vonhögen, I. G. C., el Azzouzi, H., Olieslagers, S., Vasilevich, A., de Boer, J., Tinahones, F. J., et al. (2020). MiR-337-3p Promotes Adipocyte Browning by Inhibiting TWIST1. Cells 9, 1056. doi:10.3390/cells9041056

CrossRef Full Text | Google Scholar

Wang, L., Liang, W., Wang, S., Wang, Z., Bai, H., Jiang, Y., et al. (2020). Circular RNA Expression Profiling Reveals that Circ-PLXNA1 Functions in Duck Adipocyte Differentiation. PLoS ONE 15, e0236069. doi:10.1371/journal.pone.0236069

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Z., Bian, C., Zhou, H., Huang, S., Wang, S., Liao, L., et al. (2011). MicroRNA Hsa-miR-138 Inhibits Adipogenic Differentiation of Human Adipose Tissue-Derived Mesenchymal Stem Cells through Adenovirus EID-1. Stem Cell Development 20, 259–267. doi:10.1089/scd.2010.0072

CrossRef Full Text | Google Scholar

Yao, Q., Liu, J., Xiao, L., and Wang, N. (2019). Sonic Hedgehog Signaling Instigates High-Fat Diet-Induced Insulin Resistance by Targeting PPARγ Stability. J. Biol. Chem. 294, 3284–3293. doi:10.1074/jbc.ra118.004411

CrossRef Full Text | Google Scholar

Zhang, P., Chao, Z., Zhang, R., Ding, R., Wang, Y., Wu, W., et al. (2019). Circular RNA Regulation of Myogenesis. Cells 8. doi:10.3390/cells8080885

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, P.-p., Han, Q., Sheng, M.-x., Du, C.-y., Wang, Y.-l., Cheng, X.-f., et al. (2021). Identification of Circular RNA Expression Profiles in White Adipocytes and Their Roles in Adipogenesis. Front. Physiol. 12. doi:10.3389/fphys.2021.728208

CrossRef Full Text | Google Scholar

Zhang, P., Xu, H., Li, R., Wu, W., Chao, Z., Li, C., et al. (2018). Assessment of Myoblast Circular RNA Dynamics and its Correlation with miRNA during Myogenic Differentiation. Int. J. Biochem. Cel Biol. 99, 211–218. doi:10.1016/j.biocel.2018.04.016

CrossRef Full Text | Google Scholar

Zhang, T., Zhang, Z., Xia, T., Liu, C., and Sun, C. (2019). circNrxn2 Promoted WAT Browning via Sponging miR-103 to Relieve its Inhibition of FGF10 in HFD Mice. Mol. Ther. - Nucleic Acids 17, 551–562. doi:10.1016/j.omtn.2019.06.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Zhan, S., Yang, S., Zhong, T., Guo, J., Cao, J., et al. (2021). Dynamic Expression Profiles of Circular RNAs during Brown to White Adipose Tissue Transformation in Goats (Capra hircus). Animals 11, 1351. doi:10.3390/ani11051351

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Guo, X., Pei, J., Chu, M., Ding, X., Wu, X., et al. (2020). CircRNA Expression Profile during Yak Adipocyte Differentiation and Screen Potential circRNAs for Adipocyte Differentiation. Genes 11, 414. doi:10.3390/genes11040414

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, Q., Huang, J., Wei, J., and Wu, R. (2019). Circular RNA CDR1as Sponges miR-7-5p to Enhance E2F3 Stability and Promote the Growth of Nasopharyngeal Carcinoma. Cancer Cel Int 19, 252. doi:10.1186/s12935-019-0959-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, S., Wang, J., Zhang, Q., Xu, H., and Feng, J. (2018). CircPrimer: a Software for Annotating circRNAs and Determining the Specificity of circRNA Primers. BMC Bioinformatics 19, 292. doi:10.1186/s12859-018-2304-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Zoltewicz, J. S., Stewart, N. J., Leung, R., and Peterson, A. S. (2004). Atrophin 2 Recruits Histone Deacetylase and Is Required for the Function of Multiple Signaling Centers during Mouse Embryogenesis. Development 131, 3–14. doi:10.1242/dev.00908

PubMed Abstract | CrossRef Full Text | Google Scholar