Integrative bioinformatics analysis to screen key genes and signalling pathways related to ferroptosis in obesity

ABSTRACT

Ferroptosis is closely associated with the development of disease in the body. However, there are few studies on ferroptosis-related genes (FRGs) in obesity. Therefore, key genes and signalling pathways related to ferroptosis in obesity were screened. Briefly, the RNA sequencing data of obesity and the non-obesity human samples and 259 FRGs were downloaded from GEO database and FerrDb database, respectively. The obesity-related module genes were firstly screened by weighted gene co-expression network analysis (WGCNA) and crossed with differentially expressed genes (DEGs) of obesity/normal samples and FRGs to obtain obesity-ferroptosis related (OFR) DEGs. Then, key genes were screened by PPI network. Next, the correlation of key genes and differential immune cells between obesity and normal samples were further explored by immune infiltration analysis. Finally, microRNA (miRNA)-messenger RNA (mRNA), transcription factor (TF)-mRNA networks and drug-gene interaction networks were constructed. As a result, 17 OFR DEGs were obtained, which mainly participated in processes such as lipid metabolism or adipocyte differentiation. The 4 key genes, STAT3, IL-6, PTGS2, and VEGFA, constituted the network. M2 macrophages, T cells CD8, mast cells activated, and T cells CD4 memory resting had significant differences between obesity and normal samples. Moreover, 51 miRNAs and 164 drugs were predicted for 4 key genes. All in all, this study has screened 4 FRGs, including IL-6, VEGFA, STAT3, and PTGS2, in obesity patients.

KEYWORDS:

• Obesity
• ferroptosis
• key genes
• immune infiltration
• bioinformatics

1. Introduction

Obesity is a global disease in which abnormal or excessive body fat can impair health [1]. Terribly, the incidence rate of obesity has increased year by year [2]. According to the available data, around 40% of the global population are overweight and 13% are obesity, figures which have almost tripled since 1970 [3]. Obesity not only affects leptin secretion or causes insulin resistance (IR), but also has a bearing on various chronic diseases including hypertension, diabetes mellitus, non-alcoholic fatty liver disease, cardiovascular and cerebrovascular diseases, chronic kidney diseases, respiratory problems, and cancer. It seriously threatens human health and causes heavy economic and medical burdens [4–10]. However, we have neither elucidated the interaction between obesity and these diseases nor discovered effective prevention and treatment methods. Therefore, exploring the in-depth mechanisms affecting obesity is necessary to delay the onset of chronic diseases and reduce the harm to patients.

Ferroptosis is a newly described form of cell death, which is characterized by excessive iron deposition and production of reactive oxygen species (ROS) by oxidized lipids, further promoting cell damage or death [11–13]. The main mechanism of ferroptosis is the depletion of glutathione and the reduction of glutathione peroxidase 4 (GPX4) activity. Lipid oxides cannot be metabolized through the glutathione reductase reaction catalysed by GPX4, and then oxidize lipids to produce ROS, thus promoting cell death. Ferroptosis in cells has significantly changed mitochondrial morphology [14]. Recently, numerous studies have reported that ferroptosis regulates the development of various diseases, such as cancer, neurodegenerative disorders, ischaemia/reperfusion, atherosclerosis, diabetes mellitus, alcoholic liver injury, non-alcoholic steatosis hepatitis (NASH), acute lung or kidney injury, and liver or kidney fibrosis [15,16].

According to several studies, 1) there is iron metabolism disorder or iron deposition in obesity. For example, high fat diet (HFD) was reported to increase the iron content in mouse adipose tissue (AT). The serum ferritin level in obesity patients was positively correlated with waist-to-hip ratio, body mass index (BMI) and visceral fat accumulation, and elevated serum ferritin could be detected in metabolic syndrome (MS) [17]. 2) Inhibition of iron level can improve the disorder of adipocyte and systemic metabolism [18]. It’s reported that deferoxamine could suppress the progression of obesity in diabetic conditions by improving oxidative stress and inflammatory cytokines in fat and diabetic white adipose tissue (WAT) in mice [19]. 3) In addition, the production of ROS in AT of obese mice increased selectively, while the expression of NADPH oxidase increased and antioxidant enzyme decreased. Fat accumulation is associated with systemic oxidative stress in both humans and mice [20]. 4) As an energy storage organ, brown adipose tissue (BAT) contains a large number of densely arranged iron-containing mitochondria, which increases the possibility of ROS and other harmful substances aggravating the damage of mitochondria and other cellular structures [18]. 5) It’s confirmed by clinical detection or animal tests that GPX4 expression could be reduced in a variety of obesity related diseases, such as obesity related liver and kidney injury [21], obesity cardiomyopathy, atherosclerosis, type 2 diabetes [13], and also, cell damage or death caused by ferroptosis could be observed. Inflammatory reaction, pathological changes or abnormal indicators were improved to varying degrees after treatment by inhibiting iron metabolism related genes or iron death inhibitors.

However, the relationship between ferroptosis and obesity has not been fully clarified. There are still many related diseases involved in obesity that need to expound the internal mechanism. Therefore, we aimed to excavate the relationship between ferroptosis and obesity. Briefly, bioinformatics analysis was applied to screen ferroptosis-related genes (FRGs) in obesity patients to determine the relevance of FRGs to obesity and their molecular mechanisms. Finally, this study provided a theoretical basis to deal with obesity or its related diseases.

2. Methods

2.1. Ethic statement

Animal experiments were performed under a project licence (No. kmmu20221557) granted by the Animal Ethics Committee of Kunming Medical University, in compliance with Chinese national guidelines for the care and use of animals. A protocol was prepared before the study without registration. Besides, we tried to use as few animals as possible to obtain maximum experimental data. The data of 12 mice used in this study were included in the results. Last but not the least, the mice were euthanized via overdose anaesthesia.

2.2. Data source

The GSE2508 and GSE25401 datasets were downloaded from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) database. This database is the largest and most comprehensive public gene expression data resource and covers gene expression data of almost all diseases. The GSE2508 dataset included sequencing data from 10 obesity human samples and 10 non-obesity human samples of subcutaneous adipose tissue (SAT) from the GPL8300 platform. The GSE25401 dataset contained sequencing data of 30 obesity samples and 26 non-obesity control samples from the GPL6244 platform. In addition, 259 FRGs were downloaded from the FerrDb database (http://www.zhounan.org/ferrdb/legacy/index.html). The FerrDb database is the first ferroptosis related database, providing information on genes, drugs and their interactive networks related to ferroptosis regulation. This study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

2.3. Screening for differentially expressed genes (DEGs) between obesity and control samples

In this study, the ‘limma’ R package [22] was used for differential analysis between obesity and non-obesity samples from the GSE2508 dataset. Firstly, the screening threshold for DEGs was set as |log2fold change (FC)| ≥ 0.5 and P < 0.05. The volcano and heat maps of DEGs were then plotted by the ‘ggplot2’ package and the ‘pheatmap’ package. In addition, the ‘ClusterProfiler’ package [23] was responsible for gene set enrichment analysis (GSEA) of the DEGs, and q < 0.25, P < 0.05, and false discovery rate (FDR) < 0.25 were defined as statistically significant pathways.

2.4. Screening for obesity-ferroptosis related (OFR) DEGs

Co-expression network analysis was performed on 10 obesity and 10 non-obesity samples in the GSE2508 dataset using the ‘weighted gene co-expression network analysis (WGCNA)’ package [24]. First of all, the data were log standardized, and then the samples were clustered to observe whether there were outliers or abnormal values. Secondly, the network was filtered for suitable soft thresholds based on the gene clustering of all samples to approximate a scale-free distribution. Next, the obtained topological matrices were clustered to calculate the adjacency and dissimilarity coefficients between genes. After that, the genes were divided into different modules with a minimum of 30 genes per module using the dynamic shear tree algorithm. Subsequently, a heat map of module-clinical phenotype relationships was drawn to assess the correlation of each module with both obesity and non-obesity clinical traits and screen the modules with the highest correlation to obesity. Next, intra-module analysis was performed. Notably, gene significance (GS) referred to the correlation of all gene expression profiles with this modular trait. Module membership (MM) was the absolute value of the correlation between the gene and the phenotypic trait. A correlation threshold of 0.7 was set to search for obesity-related module genes. Next, the DEGs between obesity and non-obesity samples, obesity-related module genes, and FRGs were intersected to obtain OFR DEGs. Lastly, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed on OFR DEGs using the ‘Cluster Profiler’ [23] package.

2.5. Construction of the protein-protein interaction (PPI) network

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database is a database of predicted PPIs. A PPI network was constructed in this paper using the STRING database for OFR DEGs, followed by visualization of the interactions network by ‘igraph’ package. Next, the Cytoscape (https://cytoscape.org/) plug-in Molecular COmplex Detection (MCODE) was used to identify significant modules. Besides, the Cytoscape plug-in cytoHubba was employed to explore the hub genes in the PPI network by studying the key nodes in the network.

2.6. Immune infiltration analysis

In this study, the immune infiltration of samples in the GSE2508 dataset were analysed by the Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) algorithm. The proportion of 22 immune cell types in each sample was observed first. Next, the proportion of these immune cell types in obesity and control samples was compared by Wilcoxon test. Finally, the correlation between key genes and differential immune cells was analysed by ‘psych’ package combined with Pearson method, with cor > 0.3, P < 0.05 as the threshold.

2.7. Expression level validation and methylation analysis of key genes

The expression levels of key genes in the GSE25401 and GSE2508 datasets were verified, respectively. After that, the methylation levels of key genes in obesity and normal samples were analysed using DiseaseMeth 2.0 (http://diseasemeth.edbc.org/). Notably, the methylation levels of all key genes were obtained and then compared by t-test.

2.8. Building of the miRNA-mRNA and TF-mRNA networks

The key genes were uploaded into the miRTarBase (http://mirwalk.umm.uni-heidelberg.de/) and TarBase databases (http://www.microrna.gr/tarbase) to predict their target microRNAs (miRNAs). Then, transcription factors (TFs) of the key genes were predicted using the ENCODE database. Ultimately, microRNA-messenger RNA (miRNA-mRNA) and TF-messenger RNA (TF-mRNA) relationship networks were constructed.

2.9. Construction of drug-gene interaction networks

The drug-gene interaction database (DGIdb) was applied to predict potential drugs or molecular compounds interacting with key genes, and the Cytoscape software to visualize the drug gene interaction network.

2.10. Validation of the key gene expression in vivo

2.10.1. Collection of animal tissues

A total of 12 male mice aged 26 weeks were purchased from Aniphe Biolaboratory Inc. Of them, 6 mice selected to establish obese mice models (weight 41.1–41.7 g) were fed with purified high fat feed (purchased from Jiangsu Medicience Biomedical Co., Ltd., No. MD12033) for 20 weeks. In addition, the control group comprised six C57BL/6 mice of the same age fed with a standard diet. The humidity, temperature, and light/dark cycle of the experimental environment were strictly guaranteed, and the experiment was carried out in strict accordance with the national standard for experimental animals. The SAT of the groin was dissected and quickly frozen in liquid nitrogen, then the RNA was extracted or the histology was prepared.

2.10.2. Quantitative reverse transcription polymerase chain reaction experiment

For the first thing, the inguinal SAT samples (50 mg) were collected from the HFD and control mice, respectively. Next, RNA was extracted using TRIzol (Ambion, Foster City, CA, USA) and reversely transcribed using SureScript-First-strand-cDNA-synthesis-kit (Servicebio, Wuhan, China) to obtain the reverse transcription product complementary DNA (cDNA). Later, 2 × Universal Blue SYBR Green qPCR Master Mix was used for PCR amplification on a CFX96 real-time quantitative fluorescence polymerase chain reaction (qRT-PCR) instrument. The primer sequences were shown in Table 1. Finally, t-test was adopted for comparison of the key gene expression between HFD mice and control mice.

Table 1. Primer sequences.

Primer	Sequence (5’-3’)
STAT3 F	TGCGGAGAAGCATTGTGAGTG
STAT3 R	ATTTGTTGGCGGGTCTGAAGT
IL-6 F	AGTTGTGCAATGGCAATTCTGA
IL-6 R	CTGGCTTTGTCTTTCTTGTTATCTT
PTGS2 F	TCTACAACAACTCCATCCTCCTG
PTGS2 R	GAAGCGTTTGCGGTACTCATT
VEGFA F	AACGATGAAGCCCTGGAGTG
VEGFA R	TCTGCGGATCTTGGACAAAC
GAPDH F	CCTTCCGTGTTCCTACCCC
GAPDH R	GCCCAAGATGCCCTTCAGT

STAT3, signal transducer and activator of transcription 3; IL-6, interleukin-6; PTGS2, prostaglandin endoperoxidase 2; VEGFA, vascular endothelial growth factor A; GAPDH, glyceraldehyde-3-phosphate dehydrogenase; F, forward; R, reverse.

2.11. Statistics analysis

Statistics analysis were performed through GraphPad Prism 5.0 software. Student’s paired t-test was used for comparison between the two groups. P < 0.05 was considered statistically significant.

3. Results

3.1. DEGs between obesity and non-obesity control samples

In the GSE2508 dataset, there were 436 DEGs between obesity and non-obesity samples, including 333 up-regulated differential genes and 103 down-regulated differential genes (Figure 1a,b, Table S1). In terms of KEGG pathway enrichment, the viral protein pathways, interleukin 17 (IL-17) signalling pathway, and chemokine signalling pathway, were positively correlated with the obesity sample group; while metabolic pathways were positively correlated with the non-obesity group (Figure 1c).

Figure 1. Differential genes expression of obesity samples in the GSE2508 dataset. (a) volcano plot of gene expression changes, red colour indicates upregulated genes and blue colour indicates downregulated genes. (b) heat map of gene expression changes, red colour indicates upregulated genes and blue colour indicates downregulated genes. (c-f) GSEA enrichment trend diagram of obesity samples and normal samples in viral protein interaction with cytokine and cytokine receptor (c), IL-17 signalling pathway (d), chemokine signalling pathway (e), metabolic pathways (f). GSEA, gene set enrichment analysis; IL, interleukin.

3.2. OFR DEGs enriched in 12 GO terms and 10 KEGG pathways

Due to acting as an outlier, GSM47334 was removed at the outlier position (Figure 2a). With a soft threshold of 7, the network was in a scale-free distribution (Figure 2b). The genes were clustered into 16 modules through the dynamic shear tree algorithm (Figure 2c). The black (cor = 0.83) and blue modules (cor = 0.76) had the highest positive correlation with obesity, and the brown module (cor = 0.85) had the highest negative correlation with obesity, so these three modules were selected for the next step of analysis (Figure 2d). The relationship between MM and GS of different modules showed a high correlation between the three modules and obesity (Figure 2e). A total of 17 OFR DEGs were obtained (Figure 3a,b) and enriched in 12 GO terms and 10 KEGG pathways, mainly including response to oxidative stress, regulation of angiogenesis, and regulation of lipid metabolic processes, and pathways such as the IL-17 signalling pathway (Figure 4a,b). The chord diagram and Cluego network diagram showed that PTGS2, IL-6, TNFAIP3, TXNRD1, SP1, and MAP3K5 were involved in response to oxidative stress; PTGS2, IL-6, CD44, and VEGFA were related to regulation of peptidyl-serine phosphorylation; PTGS2, TXNRD1, and SP1 participated in the regulation of lipid metabolic processes; PTGS2 and ABCC1 were involved in fatty acid derivative metabolic process (Figure 4c,d). In addition, PTGS2 participated in Kaposi sarcoma-associated herpesvirus infection, the Tumour Necrosis Factor (TNF) signalling pathway, the IL-17 signalling pathway, human cytomegalovirus infection and miRNAs in cancer; and MAP3K5 were involved in TNF signalling pathway (Figure 4e,f).

Figure 2. The co-expression modules analysis in the GSE2508 dataset. (a) clustering dendrogram of 20 samples. (b) left: the relationship between the scale-free fit index and various soft-thresholding powers; right: the relationship between the mean connectivity and various soft-thresholding powers. (c) clustering dendrogram of genes, various colours represent different modules. (d) the relationship of 2 traits and 16 modules. (e) the scatterplot describing the relationship between module membership and gene significance in black, blue, brown module.

Figure 3. Screening for OFR DEGs. (a-b) the Venn diagram for the intersection of obesity related DEGs and modular genes (a) and the intersection of DEGs highly correlated with obesity and ferroptosis related genes (b). DEGs, differentially expressed genes; OFR, obesity-ferroptosis related; FRGs, ferroptosis related genes.

Figure 4. The gene functional enrichment analysis. (a) GO enrichment. (b) KEGG enrichment. (c) chord diagram for the distribution of OFR DEGs in different GO pathways. (d) Cluego network diagram for the relationship between OFR DEGs and GO pathways. (e) chord diagram for the distribution of OFR DEGs in different KEGG pathways. (f) Cluego network diagram for the relationship between OFR DEGs and KEGG pathways. GO, gene Ontology; KEGG, Kyoto Encyclopedia of genes and Genomes; OFR DEGs, obesity-ferroptosis related differentially expressed genes; BP, biological progress; MF, molecular function; TNF, tumour necrosis factor; IL, interleukin; HIF, hypoxia inducible factor; EGFR, epidermal growth factor receptor.

3.3. The PPI network

The PPI network was composed of 16 nodes and 51 edges (Figure 5a), with 13 up-regulated genes and 3 down-regulated genes. Besides, the top four key genes in the PPI network were STAT3, IL-6, PTGS2, and VEGFA, and the sub-network of the four key genes was shown in Figure 5b.

Figure 5. The PPI analysis of key genes. (a) the PPI network. (b) top four hub genes. PPI, protein-protein interaction.

3.4. Immune infiltration analysis

The percentage of infiltration of immune cells in normal and obesity samples was shown in Figure 6a. Thereinto, T cells CD8 accounted for 41.69% in normal samples; besides, T cells CD8 (26.18% also served as the most represented cells in obesity samples. There were four immune cells, namely, M2 macrophages, T cells CD8, mast cells activated, and T cells CD4 memory resting, that differed between obesity and normal samples (Figure 6b). Furthermore, IL-6 was positively correlated with mast cells activated and M2 macrophages. PTGS2 showed a positive correlation with M2 macrophages, T cells CD4 memory resting, and mast cells activated, and a negative correlation with T cells CD8. STAT3 was positively correlated with mast cells activated and negatively correlated with T cells CD8. VEGFA had a positive correlation with T cells CD8 and a negative correlation with mast cells activated. Overall, these key genes may directly or indirectly influence immune cells to participate in the development of the disease (Figure 6c).

Figure 6. The gene functional enrichment analysis. (a) GO enrichment. (b) KEGG enrichment. (c) chord diagram for the distribution of OFR DEGs in different GO pathways. (d) Cluego network diagram for the relationship between OFR DEGs and GO pathways. (e) chord diagram for the distribution of OFR DEGs in different KEGG pathways. (f) Cluego network diagram for the relationship between OFR DEGs and KEGG pathways. GO, gene Ontology; KEGG, Kyoto Encyclopedia of genes and Genomes; OFR DEGs, obesity-ferroptosis related differentially expressed genes; BP, biological progress; MF, molecular function; TNF, tumour necrosis factor; IL, interleukin; HIF, hypoxia inducible factor; EGFR, epidermal growth factor receptor.

3.5. Expression level validation and methylation analysis of key genes

Compared with the normal samples, the testing and validation sets exhibited up-regulated expression of STAT3, IL-6, and PTGS2 and markedly down-regulated expression of VEGFA (Figure 7a,b). In addition, the methylation levels of IL-6, PTGS2, and VEGFA were different between obesity and non-obesity samples (Figure 7c).

Figure 7. Immune infiltration analysis. (a) percentage of infiltration of immune cells in normal and obesity samples. (b) GSE2508-GPL8300 immune cell distribution box line. (c) correlation map between key genes (IL-6, PTGS2, STAT3, and VEGFA) and immune cells.

*, P<0.05; **, P<0.01 ***, P<0.001. NK, natural killer; ns, not significant; IL, interleukin.

3.6. Networks of miRNA-mRNA and TF-mRNA

The related miRNA-mRNA and TF-mRNA of key genes were showed in Table 2. In the miRNA-mRNA network, STAT3 was regulated by 17 miRNAs including hsa-miR-130b-3p; IL-6 was regulated by hsd-miR-149-5P; PTGS2 was regulated by Hsa-miR-26a-5p, hsa-miR-26B-5P, hsa-miR-143-3P and hsa-miR-144-3P; VEGFA was regulated by hsa-miR-15a-5p and 29 other miRNAs (Figure 8a). The TF-mRNA network indicated that E2F1 and ELK4 regulated PTGS2; EP300 regulated STAT3 and IL-6; and CTCF, RAD21, ZBTB7A, and CCNT2 regulated VEGFA (Figure 8b).

Figure 8. Expression level of key genes (IL-6, PTGS2, STAT3, and VEGFA) in the testing set and verification set. (a-b) the expression levels of key genes were verified in GSE25401 dataset (testing set) (a) and GSE2508 dataset (verification set) (b). (c) methylation level of key genes (IL-6, PTGS2, STAT3, and VEGFA) between obesity and normal samples analyzed by DiseaseMeth 2.0. *, P < 0.05; **, P < 0.01; ***, P < 0.001. IL, interleukin.

Table 2. The related miRNA-mRNA and TF-mRNA of key genes.

Genes	miRNAs-mRNA	TFs-mRNA
STAT3	hsa-miR-130b-3p, hsa-miR-125b-5p, hsa-miR-519d-3p, hsa-let-7a-5p, hsa-let-7c-5p, hsa-let-7a-5p, hsa-miR-21-5p, hsa-miR-106b-5p, hsa-miR-29b-3p, hsa-miR-20a-5p, hsa-miR-20b-5p, hsa-miR-125a-5p, hsa-miR-181a-5p, hsa-miR-29a-3p, hsa-miR-93-5p, hsa-miR-106a-5p, hsa-miR-17-5p	POU2F2, STAT1, JUND, EP300
IL-6	hsa-miR-149-5p	EP300, CTCF, RAD21
PTGS2	hsa-miR-26a-5p, hsa-miR-143-3p, hsa-miR-144-3p, hsa-miR-26b-5p	E2F1, ELK4
VEGFA	hsa-miR-29a-3p, hsa-miR-29b-3p, hsa-miR-29d-3p, hsa-miR-106a-5p, hsa-miR-106b-5p, hsa-miR-20b-5p, hsa-miR-93-5p, hsa-miR-125a-5p, hsa-miR-15a-5p, hsa-miR-17-5p, hsa-miR-424-5p, hsa-miR-140-5p, hsa-miR-206, hsa-miR-361-5p, hsa-miR-1-3p hse-miR-302d-3p, hsa-miR-6838-5p, hsa-miR-205-5p, hsa-miR-199a-5p, hsa-miR-195-5p, hsa-miR-156-5 g, hsa-miR-16-5p, hsa-miR-200b-3p, hsa-miR-200c-3p, hsa-miR-373-3p, hsa-miR-497-5p, hsa-miR-299-5p, hsa-miR-372-3p, hsa-miR-429	CTCF, RAD21, CCNTS, 2BTB7A

STAT3, signal transducer and activator of transcription 3; IL-6, interleukin-6; PTGS2, prostaglandin endoperoxidase 2; VEGFA, vascular endothelial growth factor A.

3.7. Drug-gene interaction networks

The pharmacogenetic interaction network was shown in Figure 9. A total of 83 drugs were predicted by PTGS2, including naproxen and rofecoxib. VEGFA predicted 38 drugs including aflibercept. IL-6 predicted 25 drugs, which included siltuximab, infliximab, etanercept, and adalimumab. A total of 23 drugs, including cucurbitacin E and bardoxolone methyl, were predicted by STAT3.

Figure 9. Construction of miRNA-mRNA and TF-mRNA networks by the miRtarbase and TarBase databases. (a) gene-miRNA-mRNA network; (b) gene-TF-mRNA network. miRNA, microRNA; mRNA, messenger RNA; TF, transcription factor.

3.8. Validation of the expression of key genes in vivo

The qRT-PCR assay showed that STAT3 and PTGS2 were expressed more highly in HFD mice, and VEGFA was expressed lower in HFD mice, compared with control mice (P < 0.05). Besides, the IL-6 levels were lower in HFD mice than in control mice, but the difference was not significant (Figure 10a-d).

Figure 10. qRT-PCR assay. (a-d) the IL-6 (a), PTGS2 (b), STAT3 (c) and VEGFA (d) levels of SAT samples in mice in NC group and HFD group were detected by qRT-PCR. *, P < 0.05; **, P < 0.01. IL, interleukin; GAPDH, glyceraldehyde-3-phosphate dehydrogenase; ns, not significant; SAT, subcutaneous adipose tissue; NC, normal control; HFD, high fat diet; qRT-PCR, quantitative real-time polymerase chain reaction.

4. Discussion

In this study, the correlation between FRGs and obesity was revealed. Firstly, 17 OFR DEGs were screened using the WGCNA software package. The analysis of GO and KEGG enrichment revealed that the biological processes associated with differential genes were mainly related to oxidative stress, vascular formation, lipid metabolism, and adipocyte differentiation. After construction of the PPI network, the top four key genes (STAT3, IL-6, PTGS2, and VEGFA) were obtained. Besides, the qRT-PCR assay confirmed that STAT3 and PTGS2 were more highly expressed in HFD mice while VEGFA was more lowly expressed in HFD mice compared with control mice.

Ferroptosis is an iron-dependent regulated cell death caused by iron overload and reactive oxygen species-dependent excessive accumulation of lipid peroxidation. Obesity is directly linked to changes in iron storage. Iron overload in AT is not only one of the characteristics of obesity or IR, but also can affect mitochondrial function and adiponectin production. Increased extracellular free iron stimulates the release of ROS, further promoting IR [25,26]. The above characteristics of adipocytes give us reason to believe that iron overload may function in obesity and its related diseases by inducing ferroptosis. In this study, the analyses of GO and KEGG enrichment revealed that the biological processes associated with differential genes were mainly related to oxidative stress, vascular formation, lipid metabolism, and adipocyte differentiation. Reduction of lipid synthesis often causes free fatty acid overload, thereby inducing mitochondria reactive oxygen species (ROS) generation and endoplasmic reticulum (ER) stress. FFAs stimulate nicotinamide adenosine dinucleotide phosphate oxidase (NOXs) to generate more ROS, and then leads to ferroptosis through lipid peroxidation [27]. Iron can impact adipocyte lipid handling in several ways, including lipid accumulation during adipogenesis, lipid release during lipolysis, and lipid peroxidation [25]. Moreover, iron treatment of primary adipocytes is sufficient to increase lipolysis [28] and reduce insulin-mediated glucose uptake via changes in adipokine secretion and adipocyte IR [29]. Ferroptosis can induce angiogenesis, which contributes to the expansion of healthy fat pads, allowing AT to show a phenotype more alike that of BAT [30].

After construction of the PPI network, the top four key genes (STAT3, IL-6, PTGS2, and VEGFA) were obtained. To be specific, STAT3 activity was significantly elevated in visceral adipose tissue (VAT) in mice fed with HFD, compared with mice fed with a low-fat diet (LFD). Some studies have pointed out that functional ablation of Stat3 in T cells reduces diet-induced obesity, improves insulin sensitivity and glucose tolerance, and suppresses VAT inflammation [31,32]. STAT3 plays a role in obesity and its related diseases by regulating T cells mediated macrophages phenotypic conversion. In female breast cancer, fatty acid oxidation driven by STAT3 reduces glycolysis in CD8⁺ T-effector (TEFF) cells and inhibits its anti-tumour effect, which is crucial for obesity promoting tumour growth [33]. In addition, STAT3 mediates ferroptosis through binding to consensus response elements in the SLC7A11, GPX4, and FTH1 gene promoters [34]. Genetic inhibition of STAT3 activity triggers ferroptosis through lipid peroxidation and Fe2+ accumulation in gastric cancer cells, which demonstrates significant anti-tumour effects in gastric cancer cell xenograft model [34]. In the miRNA-mRNA network, STAT3 was regulated by 17 miRNAs such as hsa-miR-125a-5p, hsa-miR-125b-5p, hsa-miR-130b-3p, hsa-miR-181a-5p, and hsa-let-7a-5p. Generally, miR-125a was suppressed in the livers of diet-induced obesity (DIO) mice [35], the negative feedback of which regulated target gene STAT3. Such regulation would affect gluconeogenesis, glycogen synthesis, and adipogenesis in mouse livers [36]. In addition to directly affecting obesity or other obesity related diseases, miRNAs also exert a non-negligible function mediated by STAT3 in obesity.

PTGS2, also known as cyclooxygenase 2 (COX-2), participates in various physiological and pathological reactions and disease processes [37]. Up-regulated expression of COX-2 was reported in adipocytes of rats fed with HFD. Iron overload induces the production and peroxidation of lipid ROS through significantly up-regulating the expression of PTGS2, which corresponded to the molecular mechanism of ferroptosis [38]. Furthermore, adipocyte COX-2 is closely related to obesity-related AT inflammation. On the one hand, COX-2 can activate nuclear factor-κB (NF-κB) and Hypoxia-inducible Factor-1α (HIF-1α)-mediated inflammatory pathways through PGE2-EP3 mediated activation. On the other hand, COX-2-dependent CD74 expression and macrophage migration inhibitory factor (MIF) release can mediate M1 macrophage polarization, thereby inducing AT inflammation or systemic IR [39,40]. Therefore, it can be speculated that the expression of COX-2 is different in obesity patients and may have different impact on IR and related diseases. Moreover, Zhao et al. claimed that administration of adipose tissue macrophage exosomes (ATM Exos) induced higher levels of PTGS2, lipid ROS and cardiac injury. The possible mechanism is that miR-140-5p in ATM Exos targets the Solute Carrier Family 7 Member 11 (SLC7A11) to inhibit GSH synthesis, induce iron death and lead to heart injury [41]. In the miRNA-mRNA network, PTGS2 was regulated by four miRNAs including hsa-miR-26b-5p, hsa-miR-143-3p, hsa-miR-144-3p, and hsa-miR-26a-5p. However, there are no studies focusing on the correlation between PTGS2 and these four miRNAs in obesity or ferroptosis. Hence, their correlation deserves a further exploration.

VEGFA is mainly involved in regulating angiogenesis during homoeostasis and diseases [42]. Overexpression of VEGFA in WAT during the initial phase of HFD feeding using a doxycycline (Dox)-induced mouse model revealed that VEGFA angiogenesis contributed to the expansion of healthy fat pads, causing AT to show a phenotype more alike that of BAT. When the binding of VEGF-A to VEGF receptor 2 was suppressed in the same stage, the mouse presented with clear metabolic disorder and decrease in insulin sensitivity. However, the same VEGFA-VEGFR2 blockade led to improved insulin sensitivity and overall metabolic health. In short, these findings highlighted how angiogenesis difference in different stages of development of obesity affect the overall AT pathophysiological process [43,44]. A recent research revealed weight loss of COVID-19 patients caused by increasing of VEGF dependent-angiogenesis in the WAT [45]. Possibly, the effect of VEGF activation on weight loss can be reversely explored, namely, the adipose atrophy and body weight loss of patients may be promoted by increasing the expression of VEGFA and the number of microvessels in AT and promoting the browning phenotype of WAT. In the miRNA-mRNA network, VEGFA was regulated by 29 miRNAs including hsa-miR-205-5p, hsa-miR-15a-5p, hsa-miR-195-5p, hsa-miR-29c-3p, hsa-miR-361-5p, and hsa-miR-15b-5p, they are also closely related to obesity, IR, and glucose and lipid metabolism. In obese cases, miR-15a was down-regulated [46]. HFD can induce NF-κB and then activate miR-15a, which can down-regulate VEGFA and inhibit angiogenesis in mesenchymal stem cells (MSCs) [47]. Besides, in non-small cell lung cancer (NSCLC), the expression of miRNA-29c and VEGFA show a negative correlation. Down-regulation of miRNA-29c affects the expression of PI3K/Akt signalling pathway-related proteins, targeting VEGFA to regulate the progression of NSCLC, which suggests the possibility of involvement in the formation of obesity-related tumours [48].

IL-6, as a pleiotropic cytokine characterized with both proinflammatory and anti-inflammatory effects, may be up-regulated in obesity patients [49]. IL-6 stimulates lipolysis and is released from skeletal muscle during exercise, while exercise reduced visceral adipose tissue mass, this effect of exercise was abolished in the presence of IL-6 blockade [25,50]. A large amount of data have demonstrated that the biological effects of IL-6 vary with the source site in vivo, where adipocyte derived IL-6 can stimulate the accumulation of atm, whereas myeloid-derived IL-6 can inhibit the polarization of M1 macrophages and improve glucose and insulin tolerance [51]. IL-6 exposure caused cartilage cell ferroptosis by inducing cellular oxidative stress and disturbing iron homoeostasis [52]. The current understanding to the regulatory function of IL-6 exposure in ferroptosis is limited. In the miRNA-mRNA network, IL-6 was only regulated by hsa-miR-149-5p. In HFD mice SAT, up-regulated expression of miR-149-3p might aggravate glucose intolerance and liver steatosis [53]. There is evidence that miR-149 inhibits inflammatory cytokines such as IL-6, indirectly reduces visfatin, and then plays an anti-inflammatory role [54]. In the qRT-PCR verification, the HFD mice exhibited lower IL-6 levels than the control mice, but the difference was not significant. It needs further research to verify whether IL-6 is increased in patients with obesity.

The TF-mRNA network indicated that various TFs could regulate key genes. The E2F1 related TNF-superfamily paracrine loop may lead to AT dysfunction in obesity patients [55]. MiR-494-3p/JunD is associated with obesity-related metabolic cardiomyopathy [56]. EP300 is a major regulator of obesity gene disorder [57]. The expression of STAT1 in WAT is negatively correlated with serum blood glucose levels. The absence of STAT1 in adipocytes can reduce inflammation and improve insulin sensitivity [58].

The pharmacogenetic interaction network of the four key genes was further explored in this study. For the prediction of therapeutic drugs, experts had proposed a strategy of ‘reducing chronic inflammation’. Small molecular or monoclonal antibodies, such as etanercept, infliximab, adalimumab and salicylate, targeting inflammation to restore insulin sensitivity and β-cell function, plays a role in improving metabolic disorder and reducing cardiovascular and cerebrovascular risks in obesity patients. At present, though the effectiveness, safety and combined drug regimen of small molecular or monoclonal antibodies need to be further verified, this research field has broad prospects [59]. Similarly, the COX-2 inhibitor naproxen can effectively reduce blood lipids, change AT gene expression and plasma bile acid profile [60]. Rofecoxib can reduce PGF2α and 8-isoprostaglandin levels in urine of obese rats and the expression of proinflammatory chemokine mRNA in renal microvessel, ultimately reducing glomerular injury [61]. Curcumin has been shown to inhibit a variety of cell signalling pathways (including NF-κB, STAT3, Nrf2, ROS and COX-2) and present with activity in a variety of chronic diseases including tumour, diabetes and obesity [62]. Moreover, cucurbitacin E can reduce adipocyte hypertrophy, visceral obesity and adipogenesis gene expression by inhibiting JAK-STAT in VAT. To improve adipose tissue dysfunction by reducing leptin, TNF-α and increasing adiponectin can also be used as one of the new methods to treat metabolic diseases in the future [63]. Bardoxolone Methyl can prevent HFD-induced macrophage infiltration and inflammation, and reduce IL-6, STAT3 protein levels and TNF-α mRNA expression. Moreover, the development of insulin resistance and liver steatosis induced by long-term HFD in mice can be prevented by regulating the molecular mechanisms involved in insulin signal transduction, lipid metabolism or liver inflammation [64]. In the mouse model of brain oedema caused by middle cerebral artery occlusion, the expression of VEGFA and VEGFR2 increased in obese mice, with a more significant increase in the death rate, infarct volume, swelling and blood-brain barrier damage. The application of aflibercept reduces the increased permeability of vascular endothelial induced by VEGFA, thereby selectively reducing brain oedema after stroke in obese mice [65].

Based on the data of this study, we can further explore following different ideas to guide the diagnosis and treatment of obesity and its related diseases: 1) The adipose atrophy and body weight loss can be promoted by increasing the expression of VEGFA and promoting the browning phenotype of WAT; 2) Reducing the secretion of PTGS2 and IL-6 to decrease the release of ROS and inhibit the process of ferroptosis may delay or improve the obesity-related organ damage; 3) Suppressing PTGS2, IL-6 and STAT3 in obesity can restore anti-tumour immunity or block the growth, proliferation or metastasis of tumour cells. However, there are some limitations. For example, the data comes from public platform resources, and the accuracy is limited due to no external verification. Therefore, the potential mechanisms of key genes, miRNA or TF regulated by key genes, and related immune cells and drugs in the occurrence, development, or prevention of obesity need to be further explored. In the future, it is necessary to subject these parts to specific clinical experiments, animal models, or cell models for in-depth verification.

5. Conclusion

In this study, 4 key genes (STAT3, IL-6, PTGS2, and VEGFA) associated with ferroptosis in obesity were screened, and confirmed by qRT-PCR. Briefly, immune cells associated with key genes, such as T cells CD4 memory resting and M2 macrophages may participate the development of obesity. In the miRNA-mRNA network, STAT3 is regulated by 17 miRNAs including hsa-miR-125a-5p and hsa-miR-125b-5p. VEGFA is regulated by 29 miRNAs including hsa-miR-205-5p and hsa-miR-15a-5p. IL-6 is only regulated by hsa-miR-149-5p. PTGS2 is regulated by four miRNAs including hsa-miR-26b-5p, hsa-miR-143-3p, hsa-miR-144-3p, and hsa-miR-26a-5p. In the TF-mRNA network, E2F1 and ELK4 regulate PTGS2, EP300 regulates STAT3 and IL-6, and CTCF, RAD21, ZBTB7A, and CCNT2 regulate VEGFA. Finally, 83 drugs are predicted by PTGS2, including naproxen and rofecoxib. There are 38 drugs including aflibercept predicted by VEGFA, 25 drugs predicted by IL-6 including siltuximab, infliximab, etanercept, and adalimumab, and 23 drugs including cucurbitacin E, bardoxolone methyl predicted by STAT3. Undoubtedly, this study deserves a further exploration.

Author’s contributions

(I) Conception and design: MK Li, LQ Ma; (II) Administrative support: None; (III) Provision of study materials or patients: MK Li; (IV) Collection and assembly of data: MK Li, C Xing; (V) Data analysis and interpretation: MK Li, C Xing; (VI) Manuscript writing: All authors; (VII) Final approval of manuscript: All authors.

Availability of data and materials

The datasets used in this article are available from corresponding author on reasonable request.

Ethics approval and consent to participate

Animal experiments were performed under a project licence (No. kmmu20221557) granted by the Animal Ethics Committee of Kunming Medical University, in compliance with Chinese national guidelines for the care and use of animals.

Acknowledgments

We thank the funding source. The funders had no role in the study design, data collection, data analyses, interpretation, or writing of the manuscript. We thank our team members for their helpful discussions.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/21623945.2023.2264442.

Additional information

Funding

This work was supported by a grant from the National Natural Science Foundation of China (NSFC; No. 82160117).