Transcriptome of visceral adipose tissue identifies an inflammation-related ceRNA network that regulates obesity

Obesity is becoming an epidemic of widespread concern, but the underlying causes remain elusive. In this study, whole transcriptome RNA sequencing revealed differential profiles of noncoding (nc) RNAs and mRNAs in visceral adipose tissue from obese (BMI > 32.5 kg/m2) and lean (BMI < 20 kg/m2) individuals, with 1920 differentially expressed genes, 1466 long noncoding (lnc) RNAs, 122 micro (mi) RNAs, and 52 circular (circ) RNAs identified. Gene Set Enrichment Analysis, Gene Ontology analysis, Kyoto Encyclopedia of Genes and Genomes analysis revealed that these ncRNAs were involved in inflammation-related pathways that included cytokine–cytokine receptor interaction, the tumor necrosis factor and nuclear factor kappa B signaling pathways. The results indicated a critical role of inflammation in the pathogenesis of obesity. The network interaction of lncRNA, circRNA, and miRNA revealed a competing endogenous (ce) RNA network that was associated with inflammation. The ceRNA network included circORC5/miR-197-5p/TNFRSF10D and circNTRK2/miR-760/LAT, which were dysregulated in obese patients. In conclusion, this whole transcriptome study provided a pool of data that will be useful for identifying biomarkers of obesity and identified an obesity-associated ceRNA network that is regulated by circORC5 and circNTRK2.


Introduction
Obesity is increasing worldwide and is associated with metabolic diseases, such as diabetes, hypertension, cardiovascular disease, as well as an increased risk of malignancy [1,2]. Obesity is characterized by an imbalance between caloric intake and energy expenditure resulting in excessive accumulation of fat in visceral adipose tissue (VAT) [3]. In addition to energy storage, VAT is an endocrine organ, secreting adipose-derived cytokines, leptin, visfatin, and adiponectin [4]. Obesity is accompanied by hyperplasia and enlargement of adipocytes and increased secretion of multiple adipokines and proinflammatory cytokines, such as interleukin (IL)-1, IL-6, and tumor necrosis (TNF)-α that lead to chronic, low-grade local inflammation and insulin resistance [5,6]. Chronic inflammation and insulin resistance are the main characteristics of obesity, but details of mechanisms involved in either the inflammation or obesity remain to be uncovered [2,3,5].
Noncoding (nc) RNAs, including micro (mi) RNA, long noncoding (lnc) RNA, and circular (circ) RNA), are active in metabolic disease [7]. miRNAs are small endogenous ncRNAs with 20-22 bps. The major function of miRNAs is dependent on complementary interaction with the target mRNAs via complete or incomplete matching [7]. Several recent studies reported that dysregulated miRNAs enhanced chronic inflammation, insulin resistance, adipogenesis, and angiogenesis in VAT and thus promoted the process of obesity [8,9]. On the contrary, lncRNAs are transcripts that are more than 200 nt in length and lack an apparent proteincoding capacity [10,11]. For example, lncRNA-GAS5 regulates adipogenesis of mesenchymal stem cells by preventing 1 3 miR-18a from decreasing the expression of connective tissue growth factor (CTGF) [11]. CircRNAs are ncRNAs with a closed loop structure and reverse splicing of its 5′ and 3′ ends by covalent bonding, which is more stable than its linear RNA transcript [12,13]. SAMD4A, is a circRNA that is upregulated in VAT to activate the miR-138-5p/ZEH2 axis to promote the pathogenesis of obesity [13].
To identify regulatory mechanisms influencing obesity, we used whole transcriptome RNA sequencing (RNA-seq) to uncover the profiles of differentially expressed (DE) mRNAs and ncRNAs. The dysregulated profiles were then evaluated by Gene Set Enrichment Analysis (GSEA), Gene Ontology (GO) and Kyto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment to reveal the functions of the ncRNAs. The aim was to construct a competing endogenous (ce) RNA network that includes DE circRNA-miRNA-mRNA to explore the regulatory mechanisms of the ncRNAs and the potential relationship with obesity.

Sample collection
Clinical samples were collected from patients diagnosed with obesity or chronic appendicitis at our hospital between January and May 2021. Patients with a body mass index (BMI) ≥ 32.5 kg/m 2 were included in an obese group and those with a BMI ≤ 20 kg/m 2 were included in a normalweight group. About 100 mg of VAT (Fig. S1) was obtained during surgery, histologically confirmed, immediately frozen in liquid nitrogen, and stored at − 80 °C until RNA extraction. Thirty patients were enrolled in an obesity group and 30 in a normal-weight group. Differences in BMI, abdomen circumference, systolic pressure, diastolic pressure, fasting blood glucose, cholesterol, and triglycerides between the two groups were significant (all P < 0.001, Table S1). The study was approved by the ethics committee of our hospital (2021-scientific-367) and all patients gave their informed consent.

Whole transcriptome RNA-seq
Six VAT samples from obese patients and four samples from normal-weight patients were selected for RNA-seq, which was performed by Biomarker Technologies (Beijing, China). Total RNA was extracted with TRIzol (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions. The RNA quality and purity were determined with a Nanodrop 2000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). The cDNA libraries for circRNA, lncRNA, and small RNA were constructed separately. Briefly, 1.5 μg RNA samples were used as input material for rRNA removal using Ribo-Zero rRNA Removal Kits (Epicentre, Madison, WI, USA). Sequencing libraries were generated with Next Ultra Directional RNA Library Prep kits (Illumina, New England Biolabs, USA) following the manufacturer's instructions. First-strand cDNA was synthesized using random hexamer primer and reverse transcriptase. After adenylation of the 3′ ends of DNA fragments, Nebnext hairpin loop structure adapters were ligated to prepare for hybridization. The library fragments were purified with AMPure XP Beads (Beckman Coulter, Beverly, MA, USA) to select insert fragments of 150-200 bp in length. The cDNA was treated with 3 μL uracil DNA glycosylase and DNA glycosylase-lyase endonuclease VIII (USER) enzyme (New England Biolabs, USA) at 37 °C for 15 min before PCR. PCR was performed with Phusion High-Fidelity DNA polymerase, universal PCR primers and Index (X) Primer. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform at Biomarker Technologies (Beijing, China), and paired-end reads were generated.
The transcriptome was assembled using Cufflinks v2.1.1 [14] and Scripture [15], based on the reads mapped to the reference genome (GRCh38). Differential expression analysis of the two groups was performed using the DESeq2 R package v.1.10.1 (https:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ DESeq2. html). RNAs found by RNA-seq with an adjusted P value < 0.05 and an absolute value of log 2 fold change of > 1 were considered as DE.

GSEA, GO and KEGG enrichment analysis
Gene sets were downloaded from the MsigDB database [16], and broad hallmarks and KEGG symbols were included for GSEA, and pathways with corrected P values of < 0.05 were considered significantly enriched. GO and KEGG enrichment analysis were conducted by DAVID (https:// david. ncifc rf. gov) and KOBAS (http:// kobas. cbi. pku. edu. cn/ kobas3). GO includes three items, biologic processes (BPs), molecular function (MF), and cellular components (CCs). The DE mRNAs that were targets of DE miRNAs, DE lncRNAs, or DE circRNAs were enriched by GSEA, GO or KEGG analysis. The top ten enriched GO items and KEGG pathways were selected and visualized by SangerBox (http:// sange rbox. com).

Protein-protein interaction (PPI) analysis
STRING (https:// www. string-db. org) was used to construct PPI networks that showed the gene interactions involved in selected KEGG pathways. An interaction score > 0.4 was used as the extraction cutoff for the PPI pairs. CytoHubba (https:// apps. cytos cape. org/ apps/ cytoh ubba) was used to select the top ten hub genes in the PPI network. The relative results were visualized by Cytoscape v3.7.0 [17].

ceRNA network
lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA networks were constructed using current hypothetical models of ceRNA regulation and the RNA-seq results. DE miRNA targets of DE circRNAs were predicted by miRanda (http:// www. mirba se. org). The target DE genes of DE miR-NAs were predicted by the miRBase database (http:// www. mirba se. org). Finally, three dysregulated circRNAs were selected for inclusion in the ceRNA network, and the results were visualized by Cytoscape v3.7.0.

Quantitative real-time (qRT)-PCR validation
A total of 60 VAT samples (including 30 in obesity group and 30 in normal group) were selected for qRT-PCR validation. Total RNA was extracted from VAT with TRIzol (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions and then reverse transcribed to cDNA with FastKing RT Kits and gDNase (Tiangen Biotech, Beijing, China). qRT-PCR was performed with SuperReal PreMix Plus (SYBR Green; Tiangen Biotech, Beijing, China) on an ABI 7500 RT-PCR system (Applied Biosystems, Foster city, CA, USA) following the protocol provided by the manufacturer. Relative RNA expression was calculated by the 2 −ΔΔCt method, the mean CT value in normal group was used as the reference, and 18 s rRNA (for circRNA), GAPDH (for mRNA), and U6 (for miRNA) were used as internal controls. The primers that were used are listed in Table S2. SPSS 25.0 (IBM Corp., Armonk, NY, USA) was used for the statistical analysis and GraphPadPrism7.0 (GraphPad Software, La Jolla, CA, USA) was used to draw the figures. Normally distributed variables were reported as means ± standard deviation and categorical data were reported as numbers and percentages (%). Between-group differences were compared by paired t-tests or Chi-square tests. P values < 0.05 were considered statistically significant.
GO enrichment confirmed that inflammatory responses were the most affected BPs (Fig. 4A). The enriched MFs and CCs identified in the GO analysis are listed in Fig. 4A. KEGG pathway enrichment analysis found that cytokine-cytokine receptor interaction; the PI3K-Akt, chemokine, TNF, and NF-kappa B signaling pathways were the most likely to be associated with obesity (Fig. 4B). GO and KEGG integrative analysis to screen the key genes associated with obesity defined 4658 genes that were involved in pathways related cancer cell proliferation, such as the PI3K-Akt, NF-kappa B, and TNF signaling pathways, and cytokine-cytokine receptor interaction (Fig. S2A, B).

Genes involved in three key obesity-associated pathways
We selected three important inflammation-related KEGG pathways, which were "cytokine-cytokine receptor interactions", "TNF signaling pathway" and "NF-kappa B signaling pathway", and the enriched genes in these three pathways were further compared with DE target genes of circRNA and miRNA, the common genes were then selected for further PPI analysis. A total of 52 common genes in obesity-associated pathways were chosen, 29 involved in cytokine-cytokine receptor interactions, 11 in the TNF signaling pathway, and 12 in the NF-kappa B signaling pathway ( Table 2). The analysis suggested that obesity was a disease associated with chronic inflammation. The PPI network for the genes involved in the pathways was constructed with STRING online tools (Fig. 5A-C).

CeRNA mechanisms associated with obesity
In order to further explore the potential regulatory involvement of DE circRNAs in the pathogenesis of obesity, we constructed three ceRNA networks that were most likely associated with inflammation. According to the baselines of co-expression relationships between ceRNAs and bioinformatic prediction of the target miRNA and mRNA of cir-cRNA, the top 3 DE circRNA did not show any significance or the co-expression relationships were not available; as a result, three other circRNAs with a significant co-expression relationship with its ceRNAs were selected for further analysis (hsa_circ_0003000, hsa_circ_0000545, and hsa_ circ_0139142). In addition, hsa_circ_0003000 was targeted by three miRNAs, miR-1226-3p, miR-197-5p, miR-6780-5p (Fig. 6A). The miRNAs regulated 123 DE genes. The overlapping genes in the three selected KEGG pathways and PPI network included LTB, CSF3R and LAT, They were in the ceRNA network and regulated by has-circ-0000545 and hascirc-0139142. has-circ-003000 sponged miR-1226-3p and miR-6780a-5p as well as miR-197-5p. The targeted genes included XCR1, CARD11, and TNFRSF10D (Fig. 6B).

Validation of obesity-related ncRNAs in the ceRNA network
We performed qRT-PCR in human samples to validate the DE findings. Has_circ-0003000 and has_circ-0139142 were significantly downregulated in VAT from obese compared with normal-weight individuals (Fig. 7A). Consistent with the ceRNA network, miR-760 and miR-197-5p were targeted by hsa_circ-0003000 and has_circ-0139142, which were significantly upregulated in obese VAT (Fig. 7B). The expression of LAT and TNFRSF10D mRNA was significantly downregulated in obese VAT. Those genes were key genes for hsa_circ-0003000 and hsa_circ-0139142 regulation in the ceRNA networks (Fig. 7C). Our analysis confirmed that the inflammation-related ceRNA network had important roles in regulating obesity. Because the host genes of circ-0003000 and circ-0139142 were ORC5 and NTRK2 we name them circORC5 and circNTRK2, respectively.

Discussion
Accumulating evidence confirms that excessive energy storage in VAT leads to dysregulated secretion of adipokines and cytokines and state of chronic, low-grade inflammation and insulin resistance associated with obesity [18,19]. In this study, the VAT transcriptome revealed that DE genes and DE ncRNAs have an important role in the pathogenesis of obesity by modulating the process of inflammation and adipogenesis.
Numerous studies have illustrated the key role of miR-NAs in obesity [20]. The expression of miR-28-3p [21], miR-125-5p [22], and miR-103-5p [23] was shown to be downregulated, and the expression of miR-27 [24], miR-142-3p [24], miR-486, and miR-806 [25] upregulated in the VAT of obese patients. The DE miRNAs were shown to be involved in obesity-associated biological processes, including adipogenesis, angiogenesis, inflammation, insulin resistance, and macrophage polarization [21,23]. Our RNA-seq analysis screened novel DE miRNAs that were associated with obesity and might eventually serve as biomarkers or targets for the treatment and diagnosis of obesity.
We identified both DE lncRNAs and circRNAs that regulated the pathogenesis of obesity [26]. A ceRNA model of circRNA or lncRNA has been associated the inhibition function of miRNA downregulating mRNA [26,27]. Previous studies revealed that most lncRNAs and circRNAs function as ceRNAs to regulate gene expression and biological processes [28,29]. LncRNA/MEG3 was found to be downregulated during adipogenesis, which significantly promoted adipogenesis by modulating a miR-140-5P/PPAR-γ axis [28]. CircRNA is emerging as an interesting area because of its specific structure and function in regulating obesity [30]. At present, a total of 52 DE circRNAs have been identified, and the involvement of most of them in regulating obesity has not been investigated.
We used bioinformatic analysis to investigate the roles DE circRNAs active in a ceRNA network to regulate genes Fig. 1 The high-throughput sequencing results between Obesity and Normal group. A The general expression profiles of mRNAs, lncR-NAs, miRNAs and circRNAs; B Differentially expressed mRNAs, lncRNAs, miRNAs and circRNAs between Obesity and Normal group; C The relative location information of RNAs in chromosome, the outermost circle represents chromosome location and then mRNAs, lncRNAs, circRNAs and miRNAs from outside to inside, respectively. The expression profiles of up-and downregulated RNAs are marked red and blue, respectively, and the height of each line represents differential expression profiles [log(FDR)]; D Different RNA expression ratio in related chromosomes, each bar with different color represents different RNA, and the X-axis represents 23 chromosomes, and the Y-axis represents RNA ratio (the number of RNA in related chromosome/the total number of RNA), the higher the bar, the higher the RNA ratio ◂ associated with chronic inflammation by regulating IL-1, IL-6, and TNF-alpha [31,32]. Inflammatory factors were increased in patients with obesity [33,34]. Conversely, IL-13 and IL-33, which are anti-inflammatory cytokines, were decreased in obesity [35]. A high-fat diet with IL-33 injection improved insulin sensitivity in an animal study, which indicated that anti-inflammatory treatment could decrease insulin resistance in obesity [35,36]. In this study, KEGG enrichment showed that pathways including cytokine-cytokine receptor interactions, and TNF and NF-kappa B pathway signaling, regulated obesity-related inflammation and insulin resistance [37][38][39]. CXCL6, IL-6, IL-1B, TNFAIP3, and CARD11 have been identified as inflammatory factors involved in glucose intolerance and insulin resistance [40]. Our study demonstrated that those genes regulated the PPI and ceRNA network associated with obesity. CXCL6 was reported to be upregulated in high-fat diet-induced obese mice and promoted cell proliferation and immune responses in addition to inflammation [38,39]. In our study, CARD11 was identified as the target of circORC5 (hsa_circ_000300), which is a member of caspase recruitment domain family and promotes inflammation and insulin resistance by activated NF-kappa B signaling [40].

Conflict of interest
The authors declare no conflict of interest.

Ethical approval
The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Ethics Committee of Beijing Chaoyang Hospital, Capital Medical University (2021-scientific-367).
Informed consent Informed consent was obtained from all subjects involved in the study, written informed consent has been obtained from the patient(s) to publish this paper.