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

Introduction

Diabetes mellitus (DM) is a metabolic disease characterized by hyperglycemia.¹ 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.² 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.⁵ The leading causes of T2DM are insulin resistance and defective insulin secretion, but the pathogenesis of T2DM has not been fully elucidated.⁶ 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,⁹ GSE38642,¹⁰ and GSE44035¹¹ from GEO. The GEO2R online tool (https://www.ncbi.nlm.nih.gov/geo/geo2r) 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 (https://www.ncbi.nlm.nih.gov/geo) is a public functional genomics data repository of high throughput gene expression data, chips, and microarrays.¹² 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 (http://www.ncbi.nlm.nih.gov/geo/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.¹³ KEGG is a collection of databases dealing with genomes, diseases, biological pathways, drugs, and chemical materials.¹⁴ The two analyses were available in the Database for Annotation, Visualization, and Integrated Discovery (DAVID; http://david.ncifcrf.gov) (version 6.8),¹⁵ 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 (https://string-db.org/) is an online software of interactions of genes and proteins, including direct (physical) as well as indirect (functional) associations.¹⁶ 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,¹⁷ and the hub genes in the PPI networks were identified using cytoHubba.¹⁸ 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).

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

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

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

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

PPI Network and CytoHubba Analysis

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).

Table 4 Function and Connectivity of the Ten Hub Genes

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.

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

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.¹⁹ 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.²² 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).²³ 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.²⁴

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.²⁵ Also, it plays an essential role in the final differentiation of B-cells into immunoglobulin-secreting cells which involved in lymphocyte and monocyte differentiation.²⁶ 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 DM^27,28 and systemic juvenile rheumatoid arthritis.²⁹ It has been reported that IL6 exerts severe effects not only in inflammation and infection but also within the nervous and endocrine systems.³⁰ 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.³¹ Besides, analysis of specific components of the metabolic syndrome revealed a higher prevalence of insulin resistance and dyslipidemia among IL6 −174 G-allele carriers.³² IL6 is a critical inflammatory player altered in monocytes from T2DM patients. Both hyperinsulinemia and hyperglycemia contribute to changing the expression of these genes.³³

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.³⁴ 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.³⁵

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.³⁶ 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.³⁹ Furthermore, compared with the healthy control group, the expression of MMP1 in the blood of T2DM patients was significantly enhanced.⁴⁰

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.⁴¹ 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.⁴² IL11 has been recognized for its role in tumorigenesis and is produced by cancer-associated fibroblasts, which are up-regulated in cancer cells.⁴³ Reported that hyperglycemia can promote the secretion of TGFβ1 and IL11 and promote diabetic myocardial fibrosis.⁴⁴

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.⁴⁵ Previous studies predicted that TGM2, NF-кB, p38MAPK, TNF and CEBPA were the central regulators in the shared molecular pathways of RA and T2DM.⁴⁶ 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.