Identification of Hub Genes in Type 2 Diabetes Mellitus Using Bioinformatics Analysis.

BACKGROUND: Type 2 diabetes mellitus (T2DM) is one of the most common chronic diseases in the world with complicated pathogenesis. This study aimed to identify differentially expressed genes (DEGs) and molecular pathways in T2DM using bioinformatics analysis.
MATERIALS AND METHODS: To explore potential therapeutic targets for T2DM, we analyzed three microarray datasets (GSE50397, GSE38642, and GSE44035) acquired from the Gene Expression Omnibus (GEO) database. DEGs between T2DM islet and normal islet were picked out by the GEO2R tool and Venn diagram software. Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) to identify the pathways and functional annotation of DEGs. Then, protein-protein interaction (PPI) of these DEGs was visualized by Cytoscape with Search Tool for the Retrieval of Interacting Genes/Proteins (STRING).
RESULTS: In total, we identified 36 DEGs in the three datasets, including 32 up-regulated genes and four down-regulated genes. The improved functions and pathways of the DEGs enriched in cytokine-cytokine receptor interaction, pathways in cancer, PI3K-Akt signaling pathway, and Rheumatoid arthritis. Among them, ten hub genes with a high degree of connectivity were selected. Furthermore, via the re-analysis of DAVID, four genes (IL6, MMP3, MMP1, and IL11) were significantly enriched in the Rheumatoid arthritis pathway.
CONCLUSION: Our study, based on the GEO database, identified four significant up-regulated DEGs and provided novel targets for diagnosis and treatment of T2DM.

Introduction

Diabetes mellitus (DM) is a metabolic disease characterized by hyperglycemia.1 According to the data from the International Diabetes Federation (IDF), there were 463 million diabetic patients worldwide in 2019 and this number was expected to reach 700 million in 2045. Besides, 352 million individuals with impaired glucose tolerance were at increased risk of developing diabetes.2 In clinical practice, diabetes can be divided into different types like type 1, type 2 and pregnancy type among which over 90% of DM cases are type 2 DM (T2DM).3,4 T2DM can cause not only hyperglycemia but also serious complications such as nephropathy, neuropathy, retinopathy, stroke, and coronary heart diseases, making significant contributions to mortality.5 The leading causes of T2DM are insulin resistance and defective insulin secretion, but the pathogenesis of T2DM has not been fully elucidated.6 Increasing pieces of evidence have shown that genetic factors and environmental effects may greatly induce the onset of T2DM.7,8 For better understand the mechanism of T2DM and reduce the complications of T2DM, it is especially important to clarify the pathophysiological mechanism at the gene transcript level. Currently, high-throughput microarray technology and bioinformatics analysis have been widely used to screen changes of gene expression at the transcriptome level, which could help us identify differentially expressed genes (DEGs) and functional pathways involved in T2DM patients. Although some studies have explored the mRNA expression profiles in pancreatic islets of T2DM in Gene Expression Omnibus (GEO), few studies were performed to integrate all islet data in GEO for bioinformatics analysis.9–11

In this study, we chose three microarray datasets GSE50397,9 GSE38642,10 and GSE4403511 from GEO. The GEO2R online tool () and Venn diagram software were utilized to obtain the DEGs in the three datasets. Gene Ontology (GO) functional annotation analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed for the screened DEGs. Moreover, we established a protein–protein interaction (PPI) network of the DEGs, and hub genes were selected. The present study may be useful to explore the molecular mechanism of the pathogenesis in the diabetic islet.

Materials and Methods

Microarray Data Information

GEO () is a public functional genomics data repository of high throughput gene expression data, chips, and microarrays.12 Three gene expression datasets GSE50397, GSE38642, and GSE44035 were downloaded from GEO (Affymetrix GPL6244 platform, Affymetrix Human Gene 1.0 ST Array). The GSE50397 dataset contained 66 non-diabetic pancreatic islets samples and 11 diabetic pancreatic islets samples. GSE38642 included 54 non-diabetic pancreatic islets samples and nine diabetic pancreatic islets samples. GSE44035 contained 9 non-diabetic pancreatic islet samples and one diabetic pancreatic islet samples.

Data Processing of DEGs

The DEGs between DM and non-diabetic samples were screened using GEO2R (). GEO2R is an online GEO analysis tool, which can perform gene differential expression analysis on GEO sample data. For a gene symbol corresponding to multiple probe IDs, the average value of these probes was calculated as the representative expression level of this gene. Genes that met the criteria |log2 fold change (FC)| >1 and adj. P-value <0.05 was considered statistically significant. Then, the raw data in TXT format were checked in Venn software online to detect the commonly DEGs among the three datasets.

Gene Ontology (GO) Enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Enrichment Analysis

GO analysis is a commonly used approach for grouping genes and their RNA or protein products into GO categories to identify unique biological properties of high-throughput transcriptome or genome data.13 KEGG is a collection of databases dealing with genomes, diseases, biological pathways, drugs, and chemical materials.14 The two analyses were available in the Database for Annotation, Visualization, and Integrated Discovery (DAVID; ) (version 6.8),15 which is an online biological information database composed of an integrated biology knowledge data and analysis tools to extract vital information. In the present study, we performed a GO term analysis and KEGG pathway analysis using the DAVID online tool. Normal P value<0.05 was considered as the threshold for a significant enrichment result.

PPI Network and CytoHubba Analysis

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) 10.0 () is an online software of interactions of genes and proteins, including direct (physical) as well as indirect (functional) associations.16 In the present study, the PPI network of DEGs was constructed using the STRING database, and interaction with a combined score >0.4 as the cut-off criterion. Besides, we have drawn the PPI networks using Cytoscape,17 and the hub genes in the PPI networks were identified using cytoHubba.18 This Cytoscape plugin CytoHubba can find hub genes more accurately by using 11 algorithms. In our study, we selected the top ten genes as hub genes.

Results

Identification of DEGs in DM

After standardizing the microarray results, we identified 467 DEGs in GSE50397. Among these DEGs, 255 genes were up-regulated, and 212 genes were down-regulated. In GSE38642, we obtained 358 DEGs, among which 174 were up-regulated, and 184 were down-regulated. 645 DEGs were identified from GSE44035, among them, 383 up-regulated genes and 262 down-regulated genes. Subsequently, we used the Venn diagram to analyze the DEGs in the three datasets. (Figure 1) Finally, 36 DEGs were detected, consisting of 32 up-regulated genes and four down-regulated genes (Table 1).

Figure 1

Venn diagram of DEGs common to all three GEO datasets (GSE50397, GSE38642 and GSE44035).

Table 1

36 Differentially Expressed Genes (DEGs) Were Detected from Three Datasets, Including 32 Up-Regulated Genes and 4 Down-Regulated Genes

DEGs	Genes Name
Up-regulated	GJA1 CCL22 GFPT2 RASD2 IL24 BDKRB1 LTBP1 ABCA8 SPON1 NT5E ABI3BPSFRP4 SRGN MMP1 IL6 MFAP4 MMP3 FGF7 CEMIP THBS2 CYTIP TGFA IL7RVCAN IL33 CXCL5 PDGFRA IL11 SLIT2 MMP10 CYR61 INHBA
Down-regulated	CACNA2D1 WNT4 MCF2L2 INPP5F

GO and KEGG Enrichment Analyses

All 36 DEGs were analyzed using DAVID. The enriched GO terms were divided into the biological process (BP), molecular function (MF) and cellular component (CC). For BP, DEGs were enriched in cell-cell signaling, positive regulation of osteoblast differentiation, positive regulation of cell proliferation, wound healing, negative regulation of cell growth, and positive regulation of gene expression. MF analysis showed that the DEGs were significantly enriched in heparin-binding, growth factor activity, cytokine activity, calcium ion binding, metalloendopeptidase activity, and hyaluronic acid-binding. CC analysis showed that the DEGs were significantly enriched in the extracellular region, extracellular space, proteinaceous extracellular matrix, extracellular matrix, cell surface, and microfibril. KEGG pathway analysis showed that the DEGs were mainly enriched in cytokine–cytokine receptor interaction, Rheumatoid arthritis, pathways in cancer, Jak-STAT signaling pathway, PI3K-Akt signaling pathway, and Hematopoietic cell lineage (Tables 2 and 3).

Table 2

Gene Ontology Analysis of Differentially Expressed Genes in T2DM

Category	Term	Count	P-value
GOTERM_CC_DIRECT	GO:0005576~extracellular region	22	1.97E-14
GOTERM_CC_DIRECT	GO:0005615~extracellular space	16	4.82E-09
GOTERM_CC_DIRECT	GO:0005578~proteinaceous extracellular matrix	9	2.57E-08
GOTERM_CC_DIRECT	GO:0031012~extracellular matrix	9	5.54E-08
GOTERM_CC_DIRECT	GO:0009986~cell surface	5	1.77E-02
GOTERM_CC_DIRECT	GO:0001527~microfibril	2	1.85E-02
GOTERM_MF_DIRECT	GO:0008201~heparin binding	5	3.24E-04
GOTERM_MF_DIRECT	GO:0008083~growth factor activity	5	3.39E-04
GOTERM_MF_DIRECT	GO:0005125~cytokine activity	5	4.65E-04
GOTERM_MF_DIRECT	GO:0005509~calcium ion binding	7	3.25E-03
GOTERM_MF_DIRECT	GO:0004222~metalloendopeptidase activity	3	2.29E-02
GOTERM_MF_DIRECT	GO:0005540~hyaluronic acid binding	2	4.66E-02
GOTERM_BP_DIRECT	GO:0007267~cell-cell signaling	6	1.71E-04
GOTERM_BP_DIRECT	GO:0045669~positive regulation of osteoblast differentiation	4	2.62E-04
GOTERM_BP_DIRECT	GO:0008284~positive regulation of cell proliferation	7	3.62E-04
GOTERM_BP_DIRECT	GO:0042060~wound healing	4	6.11E-04
GOTERM_BP_DIRECT	GO:0030308~negative regulation of cell growth	4	2.02E-03
GOTERM_BP_DIRECT	GO:0010628~positive regulation of gene expression	5	2.07E-03

Table 3

The Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analysis in T2DM

Category	Term	Count	P-value	Genes
KEGG_PATHWAY	hsa04060: Cytokine–cytokine receptor interaction	7	1.44E-04	INHBA, IL6, CCL22, CXCL5, IL24, IL7R, IL11
KEGG_PATHWAY	hsa05323: Rheumatoid arthritis	5	2.19E-04	IL6, CXCL5, MMP3, MMP1, IL11
KEGG_PATHWAY	hsa05200: Pathways in cancer	7	1.87E-03	WNT4, IL6, FGF7, PDGFRA, TGFA, BDKRB1, MMP1
KEGG_PATHWAY	hsa04630: Jak-STAT signaling pathway	4	1.34E-02	IL6, IL24, IL7R, IL11
KEGG_PATHWAY	hsa04151: PI3K-Akt signaling pathway	5	2.98E-02	IL6, FGF7, PDGFRA, IL7R, THBS2
KEGG_PATHWAY	hsa04640: Hematopoietic cell lineage	3	3.64E-02	IL6, IL7R, IL11

A total of 35 nodes and 204 edges were imported in the PPI network, including 32 up-regulated genes and three down-regulated genes (Figure 2). INPP5F has no interaction with any other DEGs, so it is not included in the PPI network. Then we applied CytoHubba to calculate the degree of each protein node. The results demonstrated that the top 10 genes ranked by the node degree were IL6, MMP3, INHBA, MMP1, VCAN, FGF7, PDGFRA, THBS2, IL11, and CYR61, respectively (Figure 3) (Table 4).

Figure 2

Protein–protein interaction network of the differentially expressed genes (red circles meant up-regulated DEGs and green circles meant down-regulated DEGs).

Figure 3

Top ten hub genes with higher degree of connectivity from PPI network.

Table 4

Function and Connectivity of the Ten Hub Genes

Gene Symbol	Gene Description	Function	Degree
IL6	Interleukin-6	IL6 is implicated in a wide variety of inflammation associated disease states	25
MMP3	Matrix metallopeptidase 3	MMP3 can involve in the breakdown of extracellular matrix	23
INHBA	Inhibin subunit beta A	INHBA may involve in regulating insulin secretion	20
MMP1	Matrix metallopeptidase 1	MMP1 may related to tumor-induced angiogenesis.	20
VCAN	Versican core protein	VCAN is involved in cell adhesion, proliferation, proliferation migration and angiogenesis	19
FGF7	Fibroblast growth factor 7	FGF7 Plays an important role in the regulation of embryonic development, cell proliferation and cell differentiation	19
PDGFRA	Platelet derived growth factor receptor alpha	PDGFRA plays an essential role in the regulation of embryonic development, cell proliferation, survival and chemo taxis	19
THBS2	Thrombospondin 2	THBS2 can mediates cell-to-cell and cell-to-matrix interactions, inhibit angiogenesis.	18
IL11	Interleukin-11	IL11 can support the proliferation of hematopoietic stem cells and megakaryocyte progenitor cells.	18
CYR61	Cysteine rich angiogenic inducer 61	CYR61 may involve in angiogenesis, inflammation and matrix remodeling	17

Re-Analysis of 10 Hub Genes via KEGG Pathway Enrichment

To investigate the possible pathways of these ten potential hub genes, we re-performed KEGG pathway enrichment using the DAVID (Table 5). The results showed that four genes (IL6, MMP3, MMP1, and IL11) were significantly enriched in the Rheumatoid arthritis pathway (Figure 4).

Table 5

The Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analysis of 10 Hub Genes

Category	Term	Count	P-value	Genes
KEGG_PATHWAY	hsa05323:Rheumatoid arthritis	4	1.08E-04	IL6, MMP3, MMP1, IL11
KEGG_PATHWAY	hsa04151:PI3K-Akt signaling pathway	4	5.80E-03	IL6, FGF7, PDGFRA, THBS2
KEGG_PATHWAY	hsa05200:Pathways in cancer	4	8.35E-03	IL6, FGF7, PDGFRA, MMP1
KEGG_PATHWAY	hsa04060:Cytokine–cytokine receptor interaction	3	3.02E-02	INHBA, IL6, IL11

Figure 4

The result of re-analysis of 10 hub genes using the KEGG pathway enrichment. The red indicator represents four genes (IL6, MMP3, MMP1, and IL11) were significantly enriched in the rheumatoid arthritis pathway.

Discussion

T2DM has become a significant public health concern.19 DM not only brings about physical and mental injury to patients who suffer from it but also imposes heavy financial burdens on individuals and countries. In the current medical situation, DM is still a permanent disease.20,21 Thus, it urgently demanded potential, highly effective diagnostic, and therapeutic markers. Microarray technology enables us to explore the genetic alterations in T2DM and has proven to be a useful approach to identify new biomarkers in other diseases. Zhou et al identified the potential predictive biomarkers of diabetic peripheral neuropathy (DPN) from two GEO datasets by functional enrichment analysis, PPI network, and Receiver-Operator Characteristic (ROC) curve analysis.22 Li et al analyzed three mRNA microarray datasets using bioinformatics analysis methods, and sixteen hub genes were identified as potential diagnostic biomarkers for hepatocellular carcinoma (HCC).23 Xiong et al used the gene expression profiles from GEO datasets to identify the molecular mechanisms and the key hub genes that may contribute to synovial inflammation in patients with rheumatoid arthritis.24

To identify more useful prognostic biomarkers in T2DM, we used bioinformatical methods based on three gene expression profile datasets (GSE50397, GSE38642, and GSE44035). Via GEO2R and Venn software, we revealed a total of 36 commonly changed DEGs (|log2FC| > 1 and adj. P-value < 0.05), most [ie, 32 (88.89%)] of which were up-regulated. Then, GO and pathway enrichment analysis using DAVID methods showed that 1) for BP, DEGs were mainly enriched in cell-cell signaling, positive regulation of osteoblast differentiation, positive regulation of cell proliferation, wound healing, negative regulation of cell growth, and positive regulation of gene expression; 2) for MF, DEGs were enriched in heparin-binding, growth factor activity, cytokine activity, calcium ion binding, metalloendopeptidase activity, and hyaluronic acid-binding; 3) for CC, up-regulated DEGs were significantly enriched in extracellular region, extracellular space, proteinaceous extracellular matrix, extracellular matrix, cell surface, and microfibril. For the KEGG pathway enrichment analysis, up-regulated DEGs were significantly enriched in cytokine–cytokine receptor interaction, Rheumatoid arthritis, pathways in cancer, Jak-STAT signaling pathway, PI3K-Akt signaling pathway, and Hematopoietic cell lineage while down-regulated DEGs in no noteworthy signaling pathways (P< 0.05). Next, the Cytoscape tool visually constructed a PPI network with 35 nodes and 204 edges, which was predicted using STRING (a combined score of > 0.4). Furthermore, the top 10 hub genes were screened from the PPI network complex by CytoHubba analysis. Finally, we re-analyzed the ten genes via DAVID for KEGG pathway enrichment. Four genes (IL6, MMP3, MMP1, and IL11) were found enriched in Rheumatoid arthritis, which reaches statistical significance (P<0.05), which could be regarded as diagnostic biomarkers for T2DM.

Interleukin-6 (IL6) is a cytokine that functions in inflammation and the maturation of B cells.25 Also, it plays an essential role in the final differentiation of B-cells into immunoglobulin-secreting cells which involved in lymphocyte and monocyte differentiation.26 It is discharged into the bloodstream after muscle contraction and acts to increase the breakdown of fats and improve insulin resistance. The functioning of this gene is implicated in a wide variety of inflammation-associated disease states, including susceptibility to DM27,28 and systemic juvenile rheumatoid arthritis.29 It has been reported that IL6 exerts severe effects not only in inflammation and infection but also within the nervous and endocrine systems.30 Fosgerau et al indicated that matching levels of autoantibody (aAb)-IL6 cause obesity and hyperglycemia. A small subset of T2DM may, in part, evolve from an autoimmune attack against IL6.31 Besides, analysis of specific components of the metabolic syndrome revealed a higher prevalence of insulin resistance and dyslipidemia among IL6 −174 G-allele carriers.32 IL6 is a critical inflammatory player altered in monocytes from T2DM patients. Both hyperinsulinemia and hyperglycemia contribute to changing the expression of these genes.33

Matrix metallopeptidase 3 (MMP3) gene is a part of a cluster of MMP genes that are localized to chromosome 11q22.2. MMP3 can degrade fibronectin, laminin, gelatins of type I, III, IV, and V, activates procollagenase. A study tested the hypothesis that the rs3025058 of the MMP3 gene may be a genetic marker of subclinical carotid and coronary atherosclerosis in patients with T2DM.34 Van der Leeuw et al evaluated the ability of 23 novel biomarkers representing several pathophysiological pathways to improve the prediction of cardiovascular event (CVE) risk in patients with T2DM. Results presented that MMP3 improved CVE risk prediction in 2 separate cohorts of patients with T2DM beyond traditional risk factors.35

MMP1 gene is part of a cluster of MMP genes on chromosome 11. It is often up-regulated in several cancers, and this enzyme is involved in tumor-induced angiogenesis.36 Tumor-derived MMP1 can activate endothelial PAR1 to induce the destruction of the vascular endothelial barrier and increases vascular permeability.37,38 MMP1 has been demonstrated to be involved in developing insulin resistance and T2DM. Compared with the non-diabetic controls, MMP1 expression in vitreous specimens collected from patients with diabetic retinopathy was up-regulated. Hence, the interaction between MMP1, thrombin, PAR1, and VEGF may promote angiogenesis in these patients.39 Furthermore, compared with the healthy control group, the expression of MMP1 in the blood of T2DM patients was significantly enhanced.40

The protein encoded by interleukin-11 (IL11) gene is a member of the gp130 family of cytokines. It can promote the proliferation of hematopoietic stem cells and megakaryocyte progenitor cells, accelerate the maturation of megakaryocytes and thus increase the production of platelets.41 It also promotes the spread of hepatocytes in response to liver damage. Binding to its receptor formed by IL6ST and either IL11RA1 or IL11RA2 activates a signaling cascade that promotes cell proliferation.42 IL11 has been recognized for its role in tumorigenesis and is produced by cancer-associated fibroblasts, which are up-regulated in cancer cells.43 Reported that hyperglycemia can promote the secretion of TGFβ1 and IL11 and promote diabetic myocardial fibrosis.44

The re-analysis results showed that the Rheumatoid arthritis pathway displayed the most significant difference among all signaling pathways. Rheumatoid arthritis (RA) is a chronic autoimmune joint disease where persistent inflammation affects bone remodeling and leads to progressive bone destruction. Multiple studies have indicated that the inflammation mediators involved in RA are associated with T2DM.45 Previous studies predicted that TGM2, NF-кB, p38MAPK, TNF and CEBPA were the central regulators in the shared molecular pathways of RA and T2DM.46 Almost all DEGs were up-regulated both in RA and T2DM, indicating the interaction of inflammatory cytokines in the two diseases. This may be a new discovery of the T2DM pathological mechanism.

Conclusion

In conclusion, we identified DEGs by bioinformatics network analysis to find essential genes involved in the progression of T2DM. Our study predicted four useful DEGs (IL6, MMP3, MMP1, and IL11) and these predictions should be confirmed by a series of biological experiments in the future. These findings may provide some valuable approaches and direction for the precision medicine and molecular mechanisms of T2DM.