Comprehensive Analysis of circRNA Expression Profiles in Human Brown Adipose Tissue

Introduction

The excessive accumulation of adipose tissue is one of the leading factors of obesity.^1,2 Adipose tissues are typed as white adipose tissue (WAT), brown adipose tissue (BAT) and beige adipose tissue.³ The primary function of WAT is storing triglyceride, while BAT, in contrast to WAT, dissipates the energy of chemical bindings through uncoupling protein 1 (UCP1) that results in increased fatty acid oxidation and thermogenesis.^4,5 Biologically active BAT can promote thermogenesis through UCP1, and activation of BAT and induction of WAT browning accelerates glycolipid uptake, and is essential for improving lipid metabolism and insulin resistance.⁶ However, beige adipose tissue has the characteristics of both WAT and BAT, which can store triglyceride and promote thermogenesis after cold stimulation.^7,8 BAT has recently gained importance in metabolic research, partially attributed to its unique metabolic function.^9–11 Thus, new BAT biomarkers are required to be explored for clinical practices.

Circular RNAs (circRNAs), a special member of noncoding RNAs, are closed in structure, and formed by covalent bonding.^12–14 They are widely expressed in mammals.^14–16 Numerous studies indicated that circRNAs play a critical role in the development of various diseases.^17–21 Recent reports indicate that circRNAs are differentially expressed in mice BAT and WAT,^22,23 but the pattern of differentially expressed circRNAs in human BAT and WAT remains unclear.

In the present study, we analyzed the expression profiles of circRNA by using circRNA microarray to explore the differentially expressed circRNAs (DECs) in human BAT and WAT. Subsequently, six differentially expressed circRNAs were chosen for verifying the expression through qRT-PCR. Among them, four circRNAs were validated in human BAT and the biological functions were predicted by further bioinformatics analysis. Our study is the first comprehensive analysis of human BAT circRNA expression profiles, and identification of new biomarkers of BAT provides new strategies for the prevention and treatment of obesity and related disease.

Materials and Methods

Research Subjects

In this study, 30 patients with clavicle fractures and 3 patients with cholecystectomy were included. All patients were admitted in the Department of Orthopaedic Surgery and Department of Hepatobiliary Surgery of the Third Xiangya Hospital of Central South University from November 2021 to June 2022. All patients were excluded from diabetes, acute infection, tumor, and smoking. The adipose tissue samples were obtained from the supraclavicular fossa of 30 patients with clavicle fractures, but only 3 cases of BAT were confirmed with H&E staining and UCP1 expression. WAT samples were obtained from abdominal subcutaneous tissues of three patients with cholecystectomy surgery; those patients were of same gender, age, and weight with three BAT group patients.

The design of this research work was approved by the Ethics Review Committee of the Third Xiangya Hospital of Central South University, Changsha, China, and the research experiments were performed according to the guidelines outlined in the Declaration of Helsinki. All individual patients provided written informed consent.

H&E Staining

The paraffin-embedded adipose tissue was cut into small sections, dewaxed twice through xylene, and sections were sequentially soaked in alcohol at different concentrations followed by washing with distilled water. The sections were then stained with hematoxylin (Servicebio, Wuhan, China), washed in distilled water, dehydrated with 0.5% hydrochloric acid alcohol. Next, the sections were stained in an eosin staining solution (Servicebio, Wuhan, China), and dried after washing with distilled water. Lastly, neutral resin (Servicebio, Wuhan, China) was used to seal the sections, and photos were taken under a microscope (Leica DM IL LED Fluo, Germany).

Total RNA Extraction with RNA Labeling and Hybridization

Total RNA materials were obtained from human BAT and WAT samples by using trizol (Invitrogen, Carlsbad, CA, USA). The total RNA was then treated with Rnase R (Epicentre, Madison, WI, USA), and the linear RNA was removed for enriching circRNAs. The circRNAs were amplified by random primer method (Arraystar Super RNA Labeling Kit, Arraystar, Rockville, MD, USA), and then transcribed to fluorescent cRNA. Subsequently, the labeled fluorescent cRNA was recovered (Qiagen, Hilden, Germany), and then subjected to activity and concentration measurement by using a spectrophotometer (NanoDrop ND-1000, NanoDrop Technologies, Wilmington, DE, USA). Next, 5-µL 10× blocker and 1-µL 25× fragmentation buffer was added to each labeled cRNA to make them fragmented and the mixture was kept at 60°C for 30 min, followed by the addition of 25-µL 2× hybridization buffer. Next, 50 μL hybridization solution was used to assemble onto circRNA expression chip slides. The slides were then incubated for 17 hours at 65°C. Finally, after cleaning and fixing the hybridized arrays, the results were produced by using Agilent scanner G2505C (Agilent Technologies, Santa Clara, CA, USA).

CircRNA Microarray Analysis

Raw data were extracted through Agilent Feature Extraction software (Agilent Technologies, Santa Clara, CA, USA), and then uploaded to the R software limma package for quantile normalization as well as subsequent data processing. By low-intensity filtering, we found at least three out of six circRNAs for further analyses. Volcano plot and fold change cutoff were combined to identify DEG with statistical significance between two groups or samples.

Gene Ontology (GO) enrichment was conducted by using the GO seq R package. Encyclopedia of Genes and Genomes (KEGG) pathway enrichment was achieved by DAVID online tool (https://david.ncifcrf.gov/).

QRT-PCR Amplification of Partial circRNAs

Human BAT and WAT were analyzed based on raw data and validated using qRT-PCR. cDNA was reverse transcribed by Revert Aid TMM-MuLV reverse transcriptase (Invitrogen, Carlsbad, CA, USA). qRT-PCR was performed by using a thermocycler (Applied Biosystems, New York, NY, USA). The primer sequences were designed by circPrimer1.2 (http://www.bioinf.com.cn/) and are listed in Supplementary Table 1.

The circRNA–microRNA–mRNA Network Prediction

For exploring the biological function of circRNAs, the circRNA–microRNA binding sites were predicted by online software miRanda and TargetScan. The DECs within all the comparisons were marked with circRNA–microRNA interaction information. The miRNA and mRNA interaction prediction was performed by using the software, DIANA TOOLS. And the circRNA–microRNA–mRNA interaction was constructed and visualized by Cytoscape 2.8.2.¹⁵

Statistics Analysis

The data analyses between two groups were performed with paired t‐test (GraphPad Prism 8.0, GraphPad, CA, USA), where P<0.05 was considered statistically significant.

Results

Differential Expression Profile of circRNAs in BAT and WAT

H&E staining was used to verify phenotypic characters of patients’ sample from which BAT or WAT obtained. WAT had a larger cell size and uni-locular lipid droplets, while in BAT there was a smaller cell size and multi-locular lipid droplet (Figure 1A). In addition, UCP1, the molecular marker of BAT, was significantly upregulated in BAT compared to WAT (Figure 1B). Next, DECs between BAT and WAT were detected with a human circRNA microarray. The box plot showed that the Log2 ratio distribution was comparable in all of the samples (Figure 2A). Thus, the results suggested that circRNAs in BAT had an expression pattern that varied largely from WAT. The hierarchical clustering results indicated that the circRNAs expression profiles of BAT were apparently different from that of WAT (Figure 2B). The Scatter plot also showed a clear discrepancy in circRNAs expression between BAT and WAT (Figure 2C). The volcano plot showed that DECs between BAT and WAT have statistical significance (multiplicity >2.0, P<0.05) (Figure 2D). Our analyses suggested that circRNA expression patterns are different in BAT compared to WAT.

Figure 1 H&E staining and UCP1 expression in BAT and WAT. (A) BAT and WAT were stained with H&E, scale bar 200μm. (B) UCP1 mRNA expression of BAT and WAT. n=3, **p<0.01.

Figure 2 DECs in human BAT and WAT. (A) The box plot displayed the normalized intensity distribution of all data. (B) Hierarchical clustering maps of DECs were shown in BAT and WAT. Red represents relatively high expression; green represents relatively low expression. (C) Scatter plots were used to assess changes in circRNAs expression between BAT and WAT. CircRNAs above the top green line and below the bottom green line indicate more than a 1.5-fold variation in circRNAs between the two comparison samples. (D) The volcano plot showed statistically significant DECs between BAT and WAT (fold change >2.0, P<0.05). (E) Classification and distribution of DECs in human chromosomes.

Based on the above results, we found that there were 353 DECs in BAT compared to WAT, of which 152 circRNAs were upregulated and 201 circRNAs were downregulated. The leading 10 upregulated and downregulated circRNAs are shown in Table 1 and Table 2, respectively. The categories and distribution of these DECs in human chromosomes were identified indicating that most of these circRNAs are located on chromosome 1–4 (Figure 2E).

Table 1 Upregulated circRNAs in Human BAT Vs WAT

Table 2 Downregulated circRNAs in Human BAT Vs WAT

Bioinformatic Analysis of DECs

We performed GO and KEGG analyses for DECs according to their host genes. In upregulated circRNAs, the GO analysis indicated that DECs in BAT were linked to organonitrogen compound catabolic process, intracellular anatomical structure, and ion binding, and so on (Figure 3A). In the downregulated circRNAs, DECs of GO analysis were enriched in the mRNA metabolic process, cytoplasm and RNA binding, and so on (Figure 3B). In upregulated circRNAs, KEGG analysis displayed that the most noticeable enriched pathways were thermogenesis, focal adhesion, and cell cycle (Figure 3C). However, in downregulated circRNAs, the most significantly enriched pathways were found propanoate metabolism, protein processing in the endoplasmic reticulum, and carbon metabolism (Figure 3D).

Figure 3 GO and pathway analysis of DECs in BAT and WAT. GO analysis of upregulated (A) and downregulated (B) DECs in BAT and WAT. Pathway analysis of upregulated (C) and downregulated (D) DECs in BAT and WAT.

Validation of Expression of DECs

The four circRNAs with most upregulated expression and two circRNAs most downregulated expression according to fold change were verified with qRT-PCR analysis, and the partial circRNAs did not show statistical difference in circRNA microarray result. Our results showed that hsa_circ_0006168, hsa_circ_26337, and hsa_circ_0007507 expressions are significantly elevated in BAT compared to WAT. In contrast, the expression of hsa_circ_0030162 was significantly declined, which is consistent with circRNAs microarray expression profiles (Figure 4A–D). However, the level of hsa_circ_0004789 and hsa_circ_0059662 had no noticeable difference between the two groups (Figure 4E and F). So, we propose that hsa_circ_0006168, hsa_circ_26337, hsa_circ_0007507, and hsa_circ_0030162 are novel biomarkers for human BAT.

Figure 4 Validated of DECs by qRT-PCR. The expression of hsa_circ_0006168 (A), hsa_circ_0026337 (B), hsa_circ_000750 (C), hsa_circ_0030162 (D), hsa_circ_0004789 (E) and hsa_circ_0059665 (F) in BAT and WAT. n=3, *p<0.05. Results are expressed as the means.

Prediction of circRNAs Targeting microRNA–mRNA Network

To further evaluate the potential function of interested circRNAs, we predicted the potential circRNA-microRNA network by TargetScan and miRanda. According to the highest prediction score, our results showed five microRNAs interacting with hsa_circ_0006168, hsa_circ_26337, hsa_circ_0007507, and hsa_circ_0030162, respectively. For instance, hsa-miR-100-3p, hsa-miR-125b-5p, hsa-miR-133a-5p, hsa-miR-143-5p and hsa-miR-148b-5p might be the target microRNAs of hsa_circ_0006168, while hsa-miR-452-3p, hsa-miR-548ag, hsa-miR-1296-5p, hsa-miR-504-5p and hsa-miR-664a-5p were predicted as the targets of hsa_circ_0007507. The interaction network of circRNA–microRNA–mRNA is summarized in Figure 5.

Figure 5 Prediction of the circRNA–microRNA–mRNA interaction network. The target microRNAs of hsa_circ_0006168 (A), hsa_circ_0026337 (B), hsa_circ_000750 (C), and hsa_circ_0030162 (D) were predicted by Arraystar’s home-made software and target mRNA of microRNA was predicted by DIANA TOOLS software. circRNA–microRNA–mRNA interaction network was visualized and constructed by Cytoscape 2.8.2.

Discussion

Adipose tissues are important targets for treating metabolic disorders, like obesity, and BAT and WAT play different roles in regulating metabolism.^24,25 CircRNAs were reported to be involved in the development and progression of various metabolic diseases.^26,27 Xu et al found that overexpression of circRNA-CDR1as is associated with diabetes improved insulin levels in mice.²⁸ Li et al reported that circScd1 overexpression promotes steatosis in non-alcoholic fatty liver.²⁹ Previous reports also showed that circRNAs are abundantly expressed in BAT and WAT of mice,^22,23 however, there is no clear report on circRNA expression profiles in human BAT and WAT. So, here we aimed to explore the DECs between BAT and WAT, and attempted to identify some novel biomarkers for metabolic diseases.

In GO and KEGG analysis of DECs, and we tried to predict the functions of circRNAs in human BAT. The results indicated that DECs are associated with cellular proteolytic metabolic processes. The most significantly enriched pathways include thermogenesis focal adhesion cell cycle, and propanoate metabolism, etc, in particular, some circRNAs are associated with the thermogenesis pathway, which is consistent with the molecular characteristics of BAT. Moreover, the cell cycle is a critical pathway to participate in the adipogenesis of BAT.³⁰

In addition, we verified six DECs by qRT-PCR, and found that hsa_circ_0006168, hsa_circ_26337, hsa_circ_0007507, and hsa_circ_0030162 expression patterns are consistent with the circRNA microarray data. hsa_circ_0006168 has been reported to participate in regulating the proliferation and invasion of glioblastoma and esophageal cancer cells.^31–33 In addition, hsa_circRNA_0007507 and its host gene, RASA1, are upregulated in gastric cancer.³⁴

Usually circRNAs act as a sponge for microRNAs.^14,35,36 Thus, we predicted the potential association network of circRNAs with its binding target microRNAs and mRNAs. For instance, miR-125b-5p, miR-133a-5p, and miR-143-5p were predicted to interact with hsa_circ_0006168. Interestingly, miR-125b-5p is negatively regulated by UCP1 expression and BAT formation in humans and mice. Furthermore, miR-125b inhibited Brite adipose formation.³⁷ Knockout of miR-133a improved BAT thermogenic gene expression including PRDM16, UCP1 and CIDEA, and promoted BAT adipogenesis.^38–40 Knockout miR-143 elevated thermogenesis and lipolysis in BAT and decrease the BAT formation.⁴¹ miR-143 was also reported to be downregulated in BAT than WAT in case of both of pre- and mature adipocytes.⁴² Although we did not find studies on other’s predictions on microRNAs of hsa_circ_26337, hsa_circ_0007507, and hsa_circ_003016, we speculate that hsa_circ_0006168 most likely plays a vital role in BAT function.

This study had some limitations too. Firstly, the sample size was small; although we collected 30 cases of adipose tissue from the patient’s supraclavicular fossa, only 3 cases of BAT were verified by H&E staining and UCP1 expression. Secondly, only a small number of circRNAs were validated. So further investigations are required into these perspectives in future research works.

Conclusion

In this study, by using circRNA microarray, we identified several DECs including hsa_circ_0006168, hsa_circ_26337, hsa_circ_0007507, and hsa_circ_0030162 between human BAT and WAT. Our findings indicated that these circRNAs biomarkers are associated with adipose tissue formation, as well as lipolysis, and they are identified as potential therapeutic targets for treating obesity and related diseases. The regulatory mechanisms of these circRNAs will be further elucidated in future studies.