Identification of biomarkers, pathways and potential therapeutic agents for white adipocyte insulin resistance using bioinformatics analysis

ABSTRACT For the better understanding of insulin resistance (IR), the molecular biomarkers in IR white adipocytes and its potential mechanism, we downloaded two mRNA expression profiles from Gene Expression Omnibus (GEO). The white adipocyte samples in two databases were collected from the human omental adipose tissue of IR obese (IRO) subjects and insulin-sensitive obese (ISO) subjects, respectively. We identified 86 differentially expressed genes (DEGs) between the IRO and ISO subjects using limma package in R software. Gene Set Enrichment Analysis (GSEA) provided evidence that the most gene sets enriched in kidney mesenchyme development in the ISO subjects, as compared with the IRO subjects. The Gene Ontology (GO) analysis indicated that the most significantly enriched in cellular response to interferon-gamma. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis revealed that the DEGs were most significantly enriched in cytokine-cytokine receptor interaction. Protein–Protein Interaction (PPI) network was performed with the STRING, and the top 10 hub genes were identified with the Cytohubba. CMap analysis found several small molecular compounds to reverse the altered DEGs, including dropropizine, aceclofenac, melatonin, and so on. Our outputs could empower the novel potential targets to treat omental white adipocyte insulin resistance, diabetes, and diabetes-related diseases.


Introduction
Insulin resistance and insulin resistance-related complication have become important causes of mortality and morbidity through the world. Many studies have ascribed insulin resistance and diabetes to obesity [1,2]. Obesity is broadly characterized as an expansion of white adipose tissue mass to reserve the excessive energy in the form of triglycerides. During the recent decades, white adipose tissue has been emerged as a metabolic regulator for its secreting adipokines including proinflammatory or anti-inflammatory factors [3]. One of the main reasons of dysfunction of white adipose tissue to the impaired suppression of lipolysis in the presence of high insulin levels, is white adipose insulin resistance that plays a critical role in the pathophysiology of diabetes, non-alcoholic fatty liver disease, diabetic cardiomyopathy and tumours [4][5][6].
Most previous studies have focused on the contrast of white adipose between the obesity and the lean [7][8][9]. However, not all obesity contributes to insulin resistance [10,11]. Hardy and his colleagues demonstrated not only that five genes including CCL2, CCL3, CCL4, CCL18 and IL8/CXCL8 were most highly expressed independent of body mass index (BMI) in the human omental adipose tissue of insulin-resistant obese (IRO) subject, as compared with insulin-sensitive obese (ISO) subjects, but that increased macrophage infiltration in the omental adipose tissue was correlated to insulin resistance. It was of great significance for their demonstration, however, the study only focused on that BMIindependent inflammation in omental adipose tissue associated with insulin resistance in morbid obesity [12].
In the present study, for the better understanding of the molecular biomarkers, the potential mechanisms and potential therapeutic agents for white adipocyte insulin resistance, diabetes, and other metabolic diseases, we downloaded two mRNA expression profiles from Gene Expression Omnibus (GEO, http://www. ncbi.nlm.nih.gov/geo/), which is an international public repository providing freely high-throughput microarray and relevant functional genomic data sets [13]. The total 30 samples of white adipocytes in two databases were collected from the human omental adipose tissue of IRO subjects and ISO subjects, respectively. With the performance of limma package in R software, 86 differentially expressed genes (DEG) which would be the novel diagnostic biomarkers, were screened between the IRO and ISO subjects. The potential mechanisms of obesity-induced insulin resistance such as kidney mesenchyme development, cellular response to interferon-gamma and cytokine-cytokine receptor interaction and so on were explored with the performance of Gene Set Enrichment Analysis (GSEA), the Gene Ontology (GO) analysis, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and Protein-Protein Interaction (PPI) network. Ten hub genes (IL6, MMP9, CXCL8, CCL4, CXCL10, PTGS2, CCL2, SELE, CCL2, BCL2A1) were identified with the Cytohubba, including three genes (CCL2, CCL4, CXCL8) which had been identified in the previous study. CMap analysis was performed to discover several small molecular compounds to reverse the altered DEGs, including dropropizine, aceclofenac, melatonin, and so on. Our output could empower the novel and more comprehensive diagnostic and therapeutic targets for omental white adipocyte insulin resistance, and white adipocyte insulin resistance-induced diabetes and other chronic metabolic diseases.

Microarray data archives
The expression profiles by an array of GSE15773 and GSE20950 were retrieved from GEO database. The samples in two databases were the human omental (for visceral) white adipocytes collected from insulinresistant obese (IRO) and insulin-sensitive obese (ISO) subjects undergoing gastric bypass surgery between 2005 and 2009 at the University of Massachusetts Medical School [12]. GSE20950 collected 10 omental samples from IRO subjects and 10 omental samples from ISO subjects, and GSE15773 contained five IRO samples and five ISO samples. Totally, 15 omental samples from IRO subjects and 15 omental samples from ISO subjects. The statistical analyses for age, gender, height, weight, BMI, total cholesterol, high-density lipoprotein (HDL) cholesterol, low-density lipoprotein (LDL) cholesterol, triglycerides, and the number of lipids lowering therapy, between ISO and IRO subjects had no statistical significance. However, fasting glucose, fasting insulin and homeostatic model of assessment for insulin resistance (HOMA2-IR) between two groups had statistical significance [12]. The expression profilings of both databases were based on GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array) platform. Series matrix files and data table header descriptions of two databases were downloaded from the GEO database to screen and verify hub genes involved in the IRO subjects.

Microarray data and degs identification
Following two databases annotated and consolidated by the performance of Perl script, sva package in R software (version 3.5.3) (University of California, Berkeley, CA) was applied for background expression value correction and data normalization [14]. DEGs with the threshold criterion of adjusted p < 0.05 and | log FC|; (fold change) >1 between the IRO and ISO subjects were screened in limma package in R software [15]. Pheatmap package in R software was subsequently performed to plot the heatmap of DEGs [16].

GO and pathway enrichment analyses
GO is a commonly used bioinformatic tool that provides comprehensive information on gene function of individual genomic products based on defined features. GO analysis of all detected genes was conducted by GSEA software (version 3.0) [17]. GSEA is a promising, widely used software package, which derives gene sets to determine different biological functions between two groups.
GO and KEGG pathway analyses of DEGs were performed via The Database for Annotation, Visualization, and Integrated Discovery (DAVID 6.8, http://david. ncifcrf.gov) [18]. The GO analysis consists of biological processes (BP), and cellular components (CC), molecular functions (MF). KEGG is a database resource for understanding high-level biological functions and utilities. Gene count >2 and p < 0.05 were set as the threshold.

PPI network creation and hub gene identification
PPI network of DEGs was constructed by Search Tool for the Retrieval of Interacting Genes (STRING10.5; https://string-db.org/) with a combined score >0.4 as the cut-off point [19]. Hub genes were identified using Cytohubba, a plug-in of Cytoscape software (Cytoscape, 3.7.1) and significant modules in the PPI network were identified by molecular complex detection (MCODE 1.5.1), another plug-in of Cytoscape software [20,21]. The parameters of DEGs clustering and scoring were set as follows: MCODE score ≥4, degree cut-off = 2, node score cut-off = 0.2, max depth = 100, and k-score = 2.

Correlation between hub genes and diabetes
Correlation between hub genes and diabetes was performed with the Attie Lab Diabetes database (http:// diabetes.wisc.edu). The Attie Lab Diabetes database is a searchable resource of the gene expression data that is used to display the gene expression profiles of different experimental groups (lean and obese BTBR mice at 4 and 10 weeks of age) in any of six tissues, including adipose [22].

CMap analysis
The Connectivity Map (CMap) (https://portals.broadin stitute.org/cmap) is an open resource that links disease, genes, and drugs by similar or opposite gene expression profiles [23]. CMap analysis is used to predict potential small molecular compounds that can reverse altered expression of DEGs in cell lines. Mean < −0.4 and p < 0.05 were set as the screening criteria.

Statistical analysis
The statistical analyses of DEGs were done in R software. The p-values in GSEA analysis were analyzed with GSEA software (version 3.0). The p-value in the correlation between hub genes and diabetes were obtained from Attie Lab Diabetes database (http://dia betes.wisc.edu). The p-values in CMap analysis were analyzed in the CMap (https://portals.broadinstitute. org/cmap). Whenever asterisks are used to indicate statistical significance, *p < 0.05, **p < 0.01, and ***p < 0.001.

Identification of DEGs related to insulin-resistant obese
To identify DEGs in the omental white adipocytes between ISO and IRO subjects, we retrieved relevant microarray expression profiles of GSE15773 and GSE20950 from GEO database. After consolidation and normalization of the microarray data, 86 DEGs between ISO and IRO subjects were screened by limma package (|logFC| >1, adjusted p < 0.05) as shown in the heatmap ( Figure 1). Among them, 14 genes were upregulated and 72 genes were downregulated ( Figure 2, Table 1).

GO enrichment analysis of all detected genes
To identify gene sets with a statistically significant difference in the omental white adipocytes between ISO and IRO subjects, GSEA was performed, which showed most enriched gene sets of all detected genes in the IRO subjects. The top-three most significantenriched gene sets negatively correlated with the IRO subjects were kidney mesenchyme development, sex determination, positive regulation of synapse assembly (Figure 3a-c), meanwhile, the top-three most significant-enriched gene sets positively correlated with the IRO subjects were leukocyte chemotaxis, chemokinemediated signalling pathway, positive regulation of inflammatory response (Figure 3d-f).

GO enrichment analysis of DEGs
To determine the biological features of DEGs, GO analysis was accomplished by DAVID online tools. The BP analysis revealed that the DEGs were major enriched in cellular response to interferon-gamma, chemokinemediated signalling pathway, cellular response to interleukin-1, non-canonical Wnt signalling pathway via JNK cascade ( Figure 4). The CC analysis showed that DEGs were enriched in extracellular space, extracellular region, extracellular exosome and proteinaceous extracellular matrix ( Figure 4). Changes in MF of DEGs were significantly enriched in chemokine activity, heparin binding, protein binding, and peptidase activity ( Figure 4).

KEGG enrichment analysis of DEGs
To explore the potential mechanism of these DEGs, KEGG pathway analysis was performed using DAVID online tools. The results of KEGG analysis revealed that DEGs were mainly involved in cytokine-cytokine receptor interaction, TNF signalling pathway, pathways in cancer, NF-kappa B signalling pathway ( Figure 5).

PPI network analysis
To identify the most significant clusters of the DEGs, PPI network of DEGs was constituted by STRING. As shown in Figure 6(a), there were 47 nodes and 102 edged in the PPI network. The most significant modules (score = 8.5) were recognized by MCODE, a plugin of Cytoscape. (Figure 6(b)).

Hub genes recognition
To identify the hub gene in the DEGs, Cytohubba, a plug-in Cytoscape was performed. All the gene code and edge were calculated. The top 10 genes were identified as hub genes (Table 2). To find the correlation between hub genes and diabetes, the Attie Lab Diabetes database was performed. BTBR mice become severely diabetic with obesity at 10 weeks of age. We checked the hub genes using the Attie Lab Diabetes database to identify the correlation between the hub genes and diabetes. We could find that the expression of CCL2, IL6, CCL4 were significantly upregulated in the 10weeks BTBR obese diabetic mice (Figure 7).

CMap analysis
To search for potential small molecular compounds to reverse altered expression of DEGs, CMap analysis was performed. The most three significant small molecular compounds were dropropizine, aceclofenac, melatonin (Table 3).

Discussion
Insulin resistance is defined as the metabolic disordered situation that even higher concentration of insulin is insufficient to control the value of glycemia. During the recent decades, white adipose tissue has been emerged as an important regulator in the metabolism. Increasing studies have discovered that white adipose insulin resistance is strongly associated with the diabetes, cardiovascular diseases, and tumorigenesis [24][25][26]. Traditionally, white adipose includes subcutaneous adipose and visceral adipose. However, metabolic disorders are associated more strongly with visceral adiposity, rather than with subcutaneous adiposity [27]. The great concern is thus given to the diagnosis and therapeutic targets of visceral insulin resistance [28]. In the present study, bioinformatic methods are promising methods to analyze the critical genes and pathways which were associated with omental white adipose insulin resistance.
In the present study, a total of 21,755 genes were included. GSEA provided evidence that the most significant-enriched gene sets negatively correlated with the IRO subjects was kidney mesenchyme development. It has been discovered that BMP7, one of the gene ontology annotations in GO kidney mesenchyme development, could augment insulin sensitivity in mice with type 2 diabetes by potentiating PI3K/AKT pathway [29]. It will provide a new perspective on the therapeutic strategy on the insulin resistance and type 2 diabetes. Otherwise, GSEA provided further evidence that inflammation played a critical role in adipocyte insulin resistance, for the gene sets in GO that positively correlated with the IRO subjects were enriched in leukocyte chemotaxis and chemokine-mediated signalling pathway.
Based on the mRNA expression data, the 86 DEGs were identified between ISO and IRO groups. The analysis of BP in GO annotation indicated the DEGs were significantly enriched in cellular response to interferon-gamma, which was consistent with the previous demonstration that interferon-gamma released from omental adipose tissue of insulin-resistant humans  impaired the response to insulin [30]. The most enriched gene set of DEGs in the BP of GO was inflammatory response, which was well consistent with the demonstration by Hardy and his colleagues [12]. The most gene set of DEGs in the CC of GO was enriched in extracellular exosome, which included 25 DEGs. Exosomes are extracellular microvesicles (30 to 150 nm in diameter) derived from various cells, transferring different proteins, non-coding RNA and coding RNA, which have been looked as diseases biomarkers or cellcell communication factors [31,32]. Increasing studies have unveiled that exosomes derived from the insulinresistant adipocyte were implicated in the skeletal muscle insulin resistance, obesity-related liver disease, atherosclerosis, and lung cancer [33][34][35][36]. Given the broad spectrum of the discoveries of the function of these exosomes, it is not surprising that exosomes derived from insulin-resistant adipocytes, functioned as independent metabolic units, which might provide a promising therapeutic target on the insulin resistance, diabetes, and related metabolic disorders [37]. The MF analysis of GO suggested that the DEGs were the most    significant enriched in protein binding, suggesting that the interaction of two or more proteins played an important role in the adipocyte insulin resistance. Additionally, KEGG enrichment analysis of DEGs showed that these DEGs were mapped in cytokinecytokine receptor interaction, TNF signalling pathway, toll-like receptor signalling pathway, pathways in cancer, all which were consistent with the previous demonstration that white adipocyte insulin resistance had cross-talking with inflammation and tumorigenesis [38,39].
In the present study, we found 10 hub genes including MMP9, IL6, CXCL8 (IL8), CCL4, CXCL10, PTGS2 (COX-2), CCL2 (MCP-1), SELE, CCL21 and BCL2A1. MMP9 has been reported to be positively correlated with omental adipocyte insulin resistance and MMP9 was decreased in response to pioglitazone [40,41]. Hoene et al. demonstrated that IL6 could induce insulin resistance and IL6 played a pivotal role in the metabolic process [42]. Kobashi and his colleagues suggested that IL-8 could induce insulin resistance via the inhibition of insulin-induced Akt phosphorylation in adipocytes [43]. Po-Shiuan et al. reported that COX-2 activation in visceral fat inflammation might crucially contribute to the development of insulin resistance and fatty liver in high-fat induced obese rats [44]. Kanda et al. demonstrated that abundance of MCP-1 mRNA in adipose tissue was increased in genetically obese diabetic (db/db) mice. Their research also revealed that insulin resistance induced by a high-fat diet was improved extensively in MCP-1 homozygous KO mice compared with WT animals and that acute expression of a dominant-negative mutant of MCP-1 ameliorated insulin resistance in db/db mice, which made it confirmed that MCP-1 played a critical role in adipocyte insulin resistance [45]. To the date, there is still no reports on the correlation between the genes of CCL4, CXCL10, SELE, BCL1A1 and CCL21 with adipocyte insulin resistance.
In the present study, we found several potential small molecular compounds to reverse the altered expression of the DEGs, which might improve white adipocyte insulin resistance. It was reported that the replacement therapy of melatonin might contribute to restore insulin resistance of cardiomyocytes and skeletal muscle [46][47][48]. Withaferin A also has been demonstrated that it played an important role on improving high-fat diet-induced obesity and palmitic acid-induced endothelial insulin resistance through attenuation of oxidative stress and inflammation [49,50]. Levodopa is known as a precursor to dopamine. The previous studies revealed that the altered dopamine turnover contributed to the behavioural disorders in brain insulin resistance, implying that dopamine might be a protective molecule [51]. However, whether dopamine and its precursor levodopa could improve    [52]. However, the other small molecular compounds have not been reported to have the function to reverse insulin resistance or diabetes. All these small molecular compounds could be explored as the novel therapeutic targets to treat insulin resistance, diabetes, and related metabolic diseases.
In the present study, though we identified 10 hub genes of adipocyte insulin resistance and potential mechanism of white adipose insulin resistance with the bioinformatic analysis, further studies are urgently demanded to validate the hub genes, and further mechanisms would be uncovered. All the output will pave way to the potential therapeutic strategy to treat insulin resistance, diabetes and related metabolic disease.

Disclosure statement
No potential conflict of interest was reported by the authors.