Identification of potential markers for type 2 diabetes mellitus via bioinformatics analysis.

Type 2 diabetes mellitus (T2DM) is a multifactorial and multigenetic disease, and its pathogenesis is complex and largely unknown. In the present study, microarray data (GSE201966) of β‑cell enriched tissue obtained by laser capture microdissection were downloaded, including 10 control and 10 type 2 diabetic subjects. A comprehensive bioinformatics analysis of microarray data in the context of protein‑protein interaction (PPI) networks was employed, combined with subcellular location information to mine the potential candidate genes for T2DM and provide further insight on the possible mechanisms involved. First, differential analysis screened 108 differentially expressed genes. Then, 83 candidate genes were identified in the layered network in the context of PPI via network analysis, which were either directly or indirectly linked to T2DM. Of those genes obtained through literature retrieval analysis, 27 of 83 were involved with the development of T2DM; however, the rest of the 56 genes need to be verified by experiments. The functional analysis of candidate genes involved in a number of biological activities, demonstrated that 46 upregulated candidate genes were involved in 'inflammatory response' and 'lipid metabolic process', and 37 downregulated candidate genes were involved in 'positive regulation of cell death' and 'positive regulation of cell proliferation'. These candidate genes were also involved in different signaling pathways associated with 'PI3K/Akt signaling pathway', 'Rap1 signaling pathway', 'Ras signaling pathway' and 'MAPK signaling pathway', which are highly associated with the development of T2DM. Furthermore, a microRNA (miR)‑target gene regulatory network and a transcription factor‑target gene regulatory network were constructed based on miRNet and NetworkAnalyst databases, respectively. Notably, hsa‑miR‑192‑5p, hsa‑miR‑124‑5p and hsa‑miR‑335‑5p appeared to be involved in T2DM by potentially regulating the expression of various candidate genes, including procollagen C‑endopeptidase enhancer 2, connective tissue growth factor and family with sequence similarity 105, member A, protein phosphatase 1 regulatory inhibitor subunit 1 A and C‑C motif chemokine receptor 4. Smad5 and Bcl6, as transcription factors, are regulated by ankyrin repeat domain 23 and transmembrane protein 37, respectively, which might also be used in the molecular diagnosis and targeted therapy of T2DM. Taken together, the results of the present study may offer insight for future genomic‑based individualized treatment of T2DM and help determine the underlying molecular mechanisms that lead to T2DM.

Introduction

Type 2 diabetes mellitus (T2DM) has become the third main chronic non-infectious disease following tumors and cardiovascular disease, and threatens human health worldwide (1). In total, ~425 million adults are currently living with diabetes in the world, with the majority of cases being T2DM (2). The International Diabetes Federation reported that in 2045, ~629 million individuals globally will suffer from diabetes, of which ~90% will be T2DM (2). It is well known that insulin resistance and pancreatic β-cell dysfunction are major pathophysiological characteristics of T2DM. Pancreatic β-cells are needed to yield more insulin to meet mounting requirements when insulin resistance occurs (3). Previous studies of pancreatic β-cells provide a basis for improved insight into the pathogenesis and pathophysiology of T2DM, as pancreatic β-cells help in the regulation of the blood glucose level (4,5). T2DM is a complex, polygenic disease that results from the interplay of environmental and genetic factors. Candidate gene association high-throughput methods have been carried out to uncover the genetic aspects of the pathogenesis of T2DM (6–8).

In recent years, single gene research and genome-wide association studies have determined genetic susceptibility genes for the increased risk of T2DM (9–11). Previous studies of gene expression in T2DM demonstrated that decreased expression of insulin (12,13), and a reduced expression of syntaxin 1A and transcription factor 7 like 2 contributed to impaired insulin secretion (12,14,15). The downregulation of FXYD domain containing ion transport regulator 2-stimulated β-cell proliferation (13) and the upregulation of genes δ like non-canonical Notch ligand 1, diacylglycerol kinase β and zinc finger MIZ-type containing 1 were implicated in T2DM (11,16). Previous genetic studies identified several dozen genes leading to monogenic diabetes due to impaired insulin secretion (17,18). These genes play a key role in pancreatic β-cell lineage, phenotype and function. How genetic and epigenetic factors are involved in β-cell development, proliferation, differentiation and function requires further investigation. Understanding of the underlying mechanisms is vital to the development of new therapeutic methods to prevent β-cell dysfunction and failure in the development of T2DM. The identification of T2DM candidate genes has been challenging in biomedical research and the majority of the genes have yet to be discovered. The aim of the present study was to contribute to research efforts to identify the biological markers and signaling pathways associated with T2DM. These molecular mechanisms may provide insight for aspects of T2DM pathogenesis or pathophysiology.

High-throughput sequencing is becoming an important tool, extensively applied in life sciences, including in cancer detection (19–21) and for identifying global gene expression changes in T2DM (22). Knowledge of the subcellular localization of proteins provides new insight into protein function and the complex pathways that modulate biological processes on a sub-cellular level, contributing to the current understanding of the proteins that interact with each other and with other molecules in the cellular environment (23). Accordingly, subcellular proteomics, as an important step to functional proteomics, has been the focus of the prediction of subcellular protein location, which is associated with molecular cell biology, proteomics, systems biology and drug discovery (24–26); it is used to better understand complex diseases (27), such as breast cancer (28), ovarian carcinoma (29), ischemic dilated cardiomyopathy (30), esophageal squamous cell carcinoma (31) and asthma (32). It was previously demonstrated that an integrative analysis of gene expression and a protein-protein interaction (PPI) network could offer insight of the molecular mechanisms of a variety of diseases (33–35). Consequently, the present study proposed a comprehensive bioinformatics analysis of gene expression data combining protein subcellular localization information and the construction of a layered PPI network (as opposed to a traditional PPI network) to identify candidate genes. Functional enrichment analyses were performed for candidate genes. A microRNA (miRNA/miR)-target gene regulatory network and transcription factor (TF)-target gene regulatory network were also constructed to identify miRNAs and TFs, which could be involved in T2DM development. The findings of the present study may help in the discovery of potentially novel predictive and prognostic markers for T2DM, and provide insight into the underlying molecular mechanisms of T2DM.

Materials and methods

Data acquisition, preprocessing and differentially expressed genes (DEGs) analysis

GSE20966 (36), the gene chip datasets of β-cells acquired from cadaver pancreases of non-diabetic subjects (control group, n=10) and T2D subjects (T2D group, n=10), was assessed and obtained from the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/). The annotation information of GeneChip was acquired from the GPL1352 Affymetrix Human X3P Array platform (Affymetrix; Thermo Fisher Scientific, Inc.). The probes were mapped to gene names based on the GPL1352 platform and the average expression value for the probes was calculated when there was more than one gene corresponding to the same probe. In the original gene expression profiles, after normalization, MATLAB 2018a (https://www.ilovematlab.cn/forum.php?mod=home) was used to identify the DEGs by value of a fold change >1.5 and a false discovery rate_Benjamini & Hochberg (fdr_BH) <0.1. The differences in gene expression between the control and T2DM subjects were assessed using hierarchical clustering and principal component analysis (PCA).

PPI network, layering and network analysis

The PPI data were retrieved from the Human Protein Reference Database v9.0 (HPRD) (37), BioGRID v3.5 (38), IntACT v4.2 (39) and STRING v10.5 (40) databases. First, single nodes, self-loops and duplicates were removed from the PPI data. Second, the total DEGs were mapped to the PPI data. To improve the reliability, only the direct interaction proteins of these DEGs were matched. Third, the integrated PPI network was visualized and analyzed using Cytoscape v3.6.1 (41). Then, the subcellular localization information for each protein in the integrated PPI network obtained from the HPRD, the UniProt database (http://www.uniprot.org/help/uniprotkb) and the Human Protein Atlas database (http://www.proteinatlas.org/) (42) was input as a node attribute. The Cerebral plug-in in Cytoscape was applied to redistribute nodes on the basis of subcellular localization without changing their interactions (43). The layered PPI network was split into five layers: Extracellular, plasma membrane, cytoplasm, nucleus and mitochondria. Hub protein nodes that encoded DEGs in the layered network with a connectivity degree >8 were screened as candidate genes.

Functional interpretations for the candidate genes

To investigate the functions of the candidate genes, functional enrichment analysis was performed using the ClueGO and the CluePedia plug-ins (44) for Cytoscape v3.6.1 software (45). ClueGO was used to decipher functionally grouped Gene Ontology (GO) (46) and pathway annotation networks to understand their implication in three different classifications; biological process (BP), molecular function (MF) and cell component (CC), in addition to the Kyoto Encyclopedia of Genes and Genomes (KEGG) (47) signaling pathway. The relationship between the terms was calculated using κ statistics and the ClueGO network was built based on the similarity of their related genes. The CluePedia plug-in is a search tool for new markers potentially associated to pathways, and can provide a broad viewpoint of a pathway using integrated experimental and in silico data. In the present study, the enrichment analysis of gene-BP and gene-pathway was statistically validated using the Cytoscape plug-ins ClueGO and CluePedia. BPs/signaling pathways were functionally split into several groups with κ score ≥0.4 and a network was constructed, where the node represents a BP/pathway and the edge between two nodes indicates that the two BPs/pathways share common genes.

Prediction of target miRNAs and TFs for the candidate genes

Genes need to interact to react to the environment of an organism, as they cannot alone to regulate the organism. Gene expression is modulated by TFs and miRNAs at the transcriptional and post-transcriptional levels. Information on TFs, miRNAs and their corresponding target genes could provide insight into the processes of T2DM. miRNet v2.0 (http://www.mirnet.ca/) (48) was used to predict the miRNAs associated with candidate genes noted in miRTarBase v7.6 (49) and miRecords (50). The 8 most captivating groups (top 15) and a minimum of two genes for each of the groups were picked as the threshold. Then, the TFs encoded by candidate genes were used for prediction coupled with human TF information (NetworkAnalyst v3.0; http://www.networkanalyst.ca) (51) noted in Binding and Expression Target Analysis v1.0.7 (BETA) (http://cistrome.org/BETA/) (52). The miRNA-target gene regulatory network and TF-target gene regulatory network were visualized using Cytoscape.

Results

Screening for DEGs

Following data preprocessing, 108 DEGs were identified to be differentially expressed in 10 control subjects and 10 T2DM subjects, with 66 upregulated and 42 downregulated genes, as presented in the heat map of the cluster analysis of DEGs, according to the cut-off criteria of fold change >1.5 and fdr_BH <0.1 (Fig. 1A; Table SI). The PCA plot demonstrated that the DEGs could roughly divide the majority of T2DM subjects from the non-diabetic controls (Fig. 1B).

Figure 1.

Hierarchical clustering and principal component analysis of DEGs between control and T2DM subjects. (A) Hierarchical clustering analysis of DEGs was performed using MATLAB software, which split samples into groups with similar models in gene expression data. Red represents upregulated genes and green represents downregulated genes. (B) Principal component analysis of control and T2DM subjects based on DEGs. Yellow dots and blue dots refer to controls and T2DM subjects, respectively. DEGs, differentially expressed genes; T2DM; type 2 diabetes mellitus.

PPI network, layering network construction and network analysis

The identification of proteins that interact directly with proteins encoded by the 108 DEGs could help understand the molecular mechanism underlying T2DM pathophysiology. In the present study, a PPI network was built from the 108 DEGs with Cytoscape and was composed of 1,546 nodes and 1,842 edges, including 52 proteins that were encoded by upregulated genes, 36 that were encoded by downregulated genes and 1,458 nodes marking proteins that were not encoded by DEGs (Fig. 2; Table I).

Figure 2.

Integrated protein-protein interaction network. Different color nodes represent the proteins that were encoded by differentially expressed genes. Red nodes are proteins that encoded upregulated genes and the green nodes are proteins that encoded downregulated genes. Pink nodes are proteins that were not encoded by significant differentially expressed genes.

Table I.

Distribution of nodes.

Localization	Upregulated	Downregulated	Unchanged	Total
Extracellular	13	5	132	150
Plasma membrane	7	9	244	260
Cytoplasm	21	8	427	456
Nucleus	8	14	516	538
Mitochondrion	3	0	139	142
Total	52	36	1,458	1,546

Subcellular protein localization is a crucial process in numerous cells. Following synthesis, proteins are transported to distinct compartments depending on their molecular function within the cell. Certain proteins are even transported to distant sites. Protein localization data can contribute to the elucidation of protein functions. The subcellular localization information for each protein in the integrated PPI network obtained from the HPRD, the UniProt database (http://www.uniprot.org/help/uniprotkb) and the Human Protein Atlas database (http://www.proteinatlas.org/) (42), was input as a node attribute. Then, the layered network was created from the PPI network using the Cerebral plug-in (43) of Cytoscape, according to the subcellular localization information of 1,546 proteins, which was split into five layers as follows: Extracellular, plasma membrane, cytoplasm, nucleus and mitochondrion (Fig. 3).

Figure 3.

Layered protein-protein interaction network. Red nodes are upregulated genes, green nodes are downregulated genes and pink nodes are unchanged genes.

The degree distribution of a network is a standard feature of scale-free networks. The degree distributions of the layered network closely followed the power law distribution, with an R2=0.865. This suggested that the integrated PPI network is a true cellular complex biological network and scale-free. The other topological parameters of the network are presented in Table II. The results also suggested that a small number of nodes are hubs with a number of links to nodes. A total of 83 DEGs were identified as hub genes with an interaction degree ≥8 and were selected as candidate genes (Table SII). The top 20 candidate genes are presented in Table III. Of these candidate genes, ISG15 ubiquitin like modifier (ISG15) had the highest degree (185), followed by phosphoenolpyruvate carboxykinase 1 (PCK1) (85) and neural precursor cell expressed, developmentally downregulated 9 (NEDD9) (50). The present study identified 83 candidate genes for T2DM (Table IV). The identified candidate genes may serve as biological markers for future T2DM treatment research.

Table II.

Topological parameters of network.

Parameters	Value
y=βxα	y=73.313×−1.347
R2	0.865
Correlation	0.946
Clustering coefficient	0.2
Network centralization	0.118
Network density	0.002

Table III.

Top 20 hub genes.

No.	Gene name	Degree	Gene expression
1	ISG15	185	Upregulated
2	PCK1	85	Upregulated
3	NEDD9	50	Downregulated
4	TMSB4X	44	Upregulated
5	SYT1	44	Upregulated
6	IGFBP3	41	Downregulated
7	NEFL	41	Upregulated
8	TENC1	37	Upregulated
9	TIMM44	36	Upregulated
10	HMMR	35	Upregulated
11	ASB9	34	Downregulated
12	TRAF3IP2	34	Downregulated
13	KLF5	32	Upregulated
14	SERPINA3	31	Upregulated
15	NEFM	30	Upregulated
16	PHC1	29	Downregulated
17	RASGRF1	29	Downregulated
18	SERPING1	28	Upregulated
19	MYCN	27	Downregulated
20	SERPINA5	26	Upregulated

Table IV.

Candidate gene identification for type 2 diabetes mellitus.

A, Upregulated genes, n=46
Gene name	PIMID
CCR4	PMID: 17244154, PMID: 12464673
CEL	PMID: 19760265
CKS2
CPA2
CTRB1
DKK3
DPYSL3
EFHD2
EFNA3
FHIT
FXYD3	PMID: 25058609
GAD1
GIPC2
HMMR
ISG15	PMID: 25031023
KANK4
KLF5
MDFIC
NCF4
NEFL
NEFM
NEU3	PMID: 17292733
NHLH2
PAIP2B
PCK1	PMID: 24089092, PMID: 25997216
PLA2G1B	PMID: 16567514
PLXNA1
PRIM2
PRSS2
REG1A
REG1B
RENBP
RPS19BP1
SERPINA3	PMID: 28150914
SERPINA5
SFTPD
SIX6	PMID: 23478426
SLC22A3
SMOC1	PMID: 28163738
SPINK1
SYT1
TCTEX1D1	PMID: 15144884
TENC1
TIMM44	PMID: 25749183
TMPRSS2	PMID: 25749183, PMID: 9419343
TMSB4X
B, Downregulated, n=37
Gene names	PIMID
SERPING1	PMID: 23277452
ANKRD23	PMID: 26398569
ASB9
CDK2AP2
CHD5
CHL1	PMID: 22768844
CTGF	PMID: 22045431
EDN3
ESM1	PMID: 27756187
FAM105A	PMID: 20644627
GALNT14
GLIS3	PMID: 23737756
IGFBP3	PMID: 22554827, PMID: 26880678
KCNG3
LPAR3
MBNL3
MCOLN3
MYCN
NEDD9
NXPH1
ODC1
PCOLCE2
PHC1
PIGA
PLA1A
PPM1E	PMID: 20801214
PPP1R1A	PMID: 25489054
RASGRF1
RCAN2
SLC2A2	PMID: 28052964
TMEM27	PMID: 24905913
TMEM37	PMID: 29185012
TRAF3IP2	PMID: 23085260
UBL3
UNC13B
ZNF610
ZNF697

PIMID, PubMed number.

Functional enrichment analysis

To clarify possible biological functions of candidate genes, and examine the relationship between the functional groups and their underlying annotations in the networks, BP enrichment analyses were performed for the 46 upregulated and 37 downregulated candidate genes using ClueGO and CluePedia. A κ score >0.4 was set as the criterion. The results are presented in Fig. 4. Specifically, for the upregulated groups, the results yielded BPs related to the activation of ‘cellular response to cytokine stimulus’, ‘chemotaxis’, ‘inflammatory response’, ‘lipid metabolic process’, ‘macromolecule catabolic process’, ‘positive regulation of biosynthetic process’, ‘neurogenesis’, ‘regulation of cell differentiation’, ‘regulation of peptidase activity’, ‘regulation of transport’, ‘response to external biotic stimulus’ and ‘response to wounding’ (Fig. 4A). For the downregulated groups, the results yielded BPs related to the activation of ‘ion transport’, ‘neuron differentiation’, ‘positive regulation of cell death’, ‘positive regulation of cell proliferation’, ‘positive regulation of cellular component biogenesis’, ‘positive regulation of signaling’, ‘regulation of cell cycle process’ and ‘regulation of ion transport’ (Fig. 4B).

Figure 4.

Enriched GO network groups. (A) BP-enrichment analysis using ClueGO and CluePedia for upregulated genes (red). (B) BP-enrichment analysis using ClueGO and CluePedia for downregulated genes (green). Each node is a BP. Edges are links between the nodes and the length of edge indicates the degree of relatedness of two processes. The most significant parent or child term per group is displayed in the ClueGO grouped network as a group title. The size of the nodes indicates enrichment significance of the GO terms. Node color indicates the class. Mixed colors indicate that the particular node is owned by multiple classes. GO, Gene Ontology; BP, biological process.

To obtain an improved understanding of the functional involvement of these candidate genes, pathway-based functional enrichment analyses was performed using ClueGO and CluePedia. A κ score >0.4 was set as the criterion. These genes were involved in pathways associated with ‘amyotrophic lateral sclerosis (ALS)’ [neurofilament light (NEFL) neurofilament medium], ‘axon guidance’ [ephrin A3 (EFNA3) and plexin A1 (PLXNA1)], ‘cellular senescence’ [insulin like growth factor binding protein 3 (IGFBP3) and TRAF3 interacting protein 2 (TRAF3IP2)], ‘complement and coagulation cascades’ [serpin family A member 5 (SERPINA5) and serpin family G member 1 (SERPING1)], ‘fat digestion and absorption’ [carboxyl ester lipase (CEL) and phospholipase A2 group IB (PLA2G1B)], ‘glucagon signaling pathway’ (PCK1 and solute carrier family 2 member 2), ‘influenza A’ [serine protease 2 (PRSS2) and transmembrane serine protease 2 (TMPRSS2)], ‘MAPK signaling pathway’ [EFNA3 and Ras protein specific guanine nucleotide releasing factor 1 (RASGRF1)], ‘neuroactive ligand-receptor interaction’ [lysophosphatidic acid receptor 3 (LPAR3) and PRSS2], ‘PI3K-Akt signaling pathway’ (EFNA3, LPAR3 and PCK1), ‘pancreatic secretion’ [CEL, carboxypeptidase A2 (CPA2), chymotrypsinogen B1 (CTRB1), PLA2G1B and PRSS2], ‘pathways in cancer’ [CDC28 protein kinase regulatory subunit 2 (CKS2) and LPAR3), ‘phagosome’ [neutrophil cytosolic factor 4 (NCF4) and surfactant protein D (SFTPD)], ‘protein digestion and absorption’ (CPA2, CTRB1 and PRSS2), ‘purine metabolism’ [fragile histidine triad diadenosine triphosphatase (FHIT) and DNA primase subunit 2), ‘Rap1 signaling pathway’ (EFNA3 and LPAR3), ‘Ras signaling pathway’ (EFNA3, phospholipase A1 member A, PLA2G1B and RASGRF1), ‘small cell lung cancer’ (CKS2 and FHIT), ‘synaptic vesicle cycle’ [synaptotagmin (SYT1) and unc-13 homolog B (UNC13B)] and ‘transcriptional misregulation in cancer’ (IGFBP3, MYCN proto-oncogene, bHLH transcription factor and TMPRSS2) (Fig. 5).

Figure 5.

Group of significant Kyoto Encyclopedia of Genes and Genomes pathways of differentially expressed genes. Each node is a main pathway and their relation to genes (red is upregulated and green is downregulated). The node size indicates the significance of the pathway and the edge between nodes indicates shared or common genes. Dissimilar colors of node indicate dissimilar functional groups. The most significant pathway of each group is highlighted in different colors.

miRNA-target gene regulatory network

The miRNAs for DEGs were predicted using the two microRNA prediction tools through miRNet. The miRNA-gene regulatory network was built, which included 22 upregulated target genes, 19 downregulated target genes and 12 miRNAs (Fig. 6). A total of 12 miRNAs were selected, including hsa-mir-335-5p (degree=15), hsa-mir-8485 (degree=4), hsa-mir-1277-5p (degree=5), hsa-mir-190a-3p (degree=5), hsa-mir-5011-5p (degree=5), hsa-mir-124-3p (degree=6), hsa-mir-7106-5p (degree=5), hsa-let-7a-5p (degree=5), hsa-mir-192-5p (degree=5), hsa-mir-26b-5p (degree=6), hsa-let-7b-5p (degree=5) and hsa-mir-98-5p (degree=5).

Figure 6.

miRNA-target gene regulatory networks. The red diamonds indicate upregulated genes, the green diamonds indicate downregulated genes and the blue squares indicate miRNAs. miRNA/miR, microRNA.

TF-target gene regulatory network

In order to understand the topology and dynamics of the transcriptional regulatory network, TFs with a P<0.05 in BETA with its target genes via network analysis were built into a TF-target gene regulatory network using Cytoscape. The network consisted of 127 edges and 66 nodes (Fig. 7). Based on the degree, the top 8 TFs were selected to be enhancers of SUZ12 polycomb repressive complex 2 subunit (SUZ12; degree=10), enhancer of zeste 2 polycomb repressive complex 2 subunit (EZH2; degree=15), BCL6 transcription repressor (BCL6; degree=9), zinc finger protein 580 (ZNF580; degree=10), Kruppel like factor 9 (KLF9; degree=8), MYC associated zinc finger protein (MAZ; degree=15), activating transcription factor 1 (ATF1; degree=12), structure specific recognition protein 1 (SSRP1; degree=10), WRN helicase interacting protein 1 (WRNIP1; degree=10), chromodomain helicase DNA binding protein 1 (CHD1; degree=10) and SMAD5 (degree=11).

Figure 7.

TF-target gene regulatory network. Red indicates upregulated genes, green indicates downregulated genes and blue indicates TFs. TF, transcription factor.

Discussion

T2DM is a multifactorial and multigenetic disease, and its pathogenesis is complex and largely unknown. A PPI network and a layered network for DEGs were constructed and it was observed that the majority of the proteins were localized in the cytoplasm, followed by the nucleus. The modules were mined from the PPI network and ISG15, PCK1, NEDD9, thymosin β 4 X-linked (TMSB4X), SYT1, IGFBP3, NEFL, tensin 2, translocase of inner mitochondrial membrane 44 (TIMM44), hyaluronan mediated motility receptor (HMMR), ankyrin repeat and SOCS box containing 9 and TRAF3IP2 were screened as the candidate genes with the highest degree of connectivity.

ISG15 has an anti-apoptotic capability on MIN6 cells (53). Upregulated ISG15 could be a potential therapeutic approach for type 1 diabetes (T1D) in pancreatic β-cells (53). PCK1 has been a candidate gene for T2DM susceptibility (54). SYT1 is a Ca2+ sensor that plays a central role in insulin release, which is a characteristic deterioration in the early stages of T2DM (55,57). Higher levels of IGFBP3 might raise the risk of T2DM (57,58). The TIMM44 gene could be a new target for T2DM therapy (59). The procession of vascular diseases can be delayed by targeting TRAF3IP2 during diabetes and atherosclerosis, as TRAF3IP2 reconciles high glucose-induced NF-κB and AP-1-dependent inflammatory signaling and endothelial dysfunction (60). TRAF3IP2 may play a role in the pathogenesis of T1D (61). Notably, to the best of the authors' knowledge, ribosomal protein S19 binding protein 1 (RPS19BP1) and SFTPD have not been previously reported to be dysregulated in T2DM. RPS19BP1 is a direct regulator of NAD-dependent deacetylase sirtuin-1 (SIRT1), which is a promising molecular target for the treatment of obesity. RPS19BP1 can serve as a prognostic indicator via the direct regulation of SIRT1 in obese patients with T2DM (62–64). SFTPD is an element of lung innate immunity that strengthens pathogen clearance and regulates inflammatory responses (65); its expression is decreased in obesity and in impaired glucose tolerance, both of which are related to the development of T2DM.

The functional assay for candidate genes using ClueGO and CluePedia in GO terms or KEGG pathways identified several molecular mechanisms, widely known to underlie the pathogenesis of T2DM. A number of vital processes/signaling pathways and key factors connected with the pathogenesis of T2DM were identified from the functional enrichment analyses. In the upregulated group, a number of the corresponding encoded proteins were distributed in the extracellular and cytoplasmic layers. In particular, it was identified that the majority of BPs were associated with ‘inflammatory response’, ‘cellular response to cytokine stimulus’ and ‘lipid metabolic process’. In previous years, a number of studies have suggested that T2DM may be a chronic inflammatory response modulated by inflammatory factors (66–69). Immune cell infiltration and high levels of cytokines were observed in the pancreas islets of T2DM (70,71), which caused differing degrees of impairment to pancreatic β-cell activity, resulting in β-cell failure (72,73). Lipotoxic effects lead to impaired insulin secretion and apoptosis of β-cells, which can give rise to the β-cell functional loss in the pathogenesis of T2DM (74). In the downregulated group, more corresponding encoded proteins were distributed in the plasma membrane and nucleus layers, and BPs were related to the ‘positive regulation of cell death’ and ‘positive regulation of cell proliferation’. Previous studies demonstrated that glucotoxic conditions are used to elevate β-cell proliferation and neogenesis, and the inhibition of apoptosis and death lead to insulin release defects, which is typical of diabetes (75,76); β-cell death is the major cause of T2DM. According to protein subcellular localization, the composition and biological value of proteins could change; analyzing PPIs may help identify the signaling pathways. The present study identified three interactions among 83 candidate genes, including SERPINA3 and CTRB1 in the extracellular matrix, SERPING1 and thymosin b4 X-linked (TMSB4X) in the cytoplasm, and SYT1 and UNC13B in the cytoplasm, which were shared between BPs, including protein input into the cytoplasm and cell-cell signaling, which might indicate that their expression was altered by the signaling cascades of the extracellular-plasma membrane-cytoplasm or nucleus and alteration in cell development. Takahashi et al (77) identified that SERPINA3 levels were significantly increased in T2DM. The rs7202877 locus for CTRB1 and CTRB2, a known diabetes risk locus, might be able to prevent diabetes via the incretin pathway (78). SERPINA5 inhibits activated protein C (APC). APC has a potential preventative role for islet β-cell damage and diabetes (79). A previous study observed that the plasma levels of APC were notably decreased in T2DM (80). SERPINA5 may be involved in T2DM via inhibited APC expression. TMSB4X is involved in cell proliferation, migration and differentiation, and its level increased in diabetic membranes (81).

The enriched KEGG pathways of candidate genes involved ‘MAPK signaling pathway’, ‘Ras signaling pathway’, ‘PI3K-Akt signaling pathway’, ‘Rap1 signaling pathway’ and ‘purine metabolism’. Previous studies demonstrated that p38 MAPK and ERK signaling were activated to inhibit obesity and associated T2DM (82,83). The Ras/Raf/ERK signaling pathway may control β-cell proliferation, and Ras is essential for normal β-cell development and function (84). Saltiel and Kahn (75) demonstrated that any obstacles in the PI3K/Akt insulin signaling pathway result in islet β-cells insulin resistance and lead to β-cell function reduction. Previous studies demonstrated that the Rap1 pathway may yield targets for β-cell dysfunction therapy in diabetes (85–88). The pathway enrichment results for candidate genes in the present study identified the MAPK. Ras, PI3K-Akt and Rap1 signaling pathways in diabetes. RASGRF1 is mainly involved in the MAPK and Ras signaling pathways (89). Suppressing the expression of Rasgrf1 may contribute to insufficient insulin secretion (90), which due to insulin resistance, causes T2DM. In addition, EFNA3 in the extracellular matrix was enriched in the PI3K-Akt, MAPK, Rap1 and Ras signaling pathways. EFNA3 is an upstream gene of the MAPK signaling pathway and the PI3K-Akt signaling pathway (91). Therefore, it was hypothesized that a low expression of EFNA3 may regulate β-cell proliferation by activating Ras/Raf/MEK/ERK. FHIT, a protein product involved in purine metabolism that participates in the T2DM pathway, is expressed in the pancreas (92). Its single-nucleotide polymorphism (rs3845971) was related to an intensified risk of T2DM (93). FHIT increases adenosine-diphosphate in the purine metabolism pathway (94). Therefore, FHIT may induce β-cell apoptosis in the pancreas due to T2DM.

miRNA-target gene interaction networks were constructed from 12 miRNAs. HMMR, ubiquitin like 3, ornithine decarboxylase 1, muscleblind like splicing regulator 3 and procollagen C-endopeptidase enhancer 2 (PCOLCE2) were regulated by hsa-miR-192-5p. Connective tissue growth factor (CTGF), family with sequence similarity 105, member A (FAM105A), MyoD family inhibitor domain containing, soluble carrier family 22 member 3, IGFBP3, CKS2 and EF-hand domain family member D2 were regulated by hsa-miR-124-3p. Protein phosphatase 1 regulatory inhibitor subunit 1 A (PPP1R1A), C-C motif chemokine receptor 4 (CCR4), SYT1, LMBR1 domain containing 2, C-type lectin domain containing 11A, NLR family pyrin domain containing 13, SERPINA3, protein phosphatase, Mg2+/Mn2+ dependent 1E (PPM1E), GIPC PDZ domain containing family member 2, NCF4, SERPING1, SIX homebox 6 CTRB1, CEL, cell adhesion molecular L1 like (CHL1), SFTPD and PLXNA1 were regulated by hsa-miR-335-5p. It was observed that hsa-miR-192-5p, hsa-miR-124-3p and hsa-miR-335-5p appeared to regulate the majority of the candidate genes identified in T2DM in the present study. It was previously identified that the downregulation of miR-192-5p usually occurs in the more extreme stages of diabetes (95). PCOLCE2, a collagen-related gene, is significantly reduced in T2DM (96). Zhu et al (97), identified that the expression level of hsa-miR-124-3p is decreased in patients with T2DM [9 high-body mass index (BMI) and 1 low-BMI] CTGF expression, a vital adjudicator of progressive pancreatic fibrosis, is elevated in T2DM (98). FAM105A is reported to be associated with T2DM (36). miR-335-5p expression was increased by islets in a diabetic Goto-Kakizaki-rat model (99). PPP1R1A has previously been identified as a potential participant and experimentally validated in the pathogenesis of islet dysfunction in T2DM (100). The ratio of CXCR3 to CCR4 receptor expression was positively correlated with the duration of T1D (r=0.947; P=0.0004) (101). The expression of SERPINA3 was increased significantly in T2DM and can be used for the early detection of T2DM (77). PPM1E is a potential drug target for diabetic therapies (102). The CTRB1/2 locus influences the susceptibility and treatment for diabetes via the incretin pathway (78). It was previously identified that mutations for the highly polymorphic CEL gene can be a rare cause of monogenic diabetes (103). CHL1 affects insulin secretion in INS-1 cells and has been identified as being potentially involved in T2DM (104). SFTPD expression was decreased in patients with T2DM (65). Therefore, has-miR-8485, has-miR-1277-3p, has-miR-190a-3p, has-miR-5011-3p, has-let-7a-5p, has-let-7b-5p, has-miR-98-5p, has-miR-7106-5p and has-miR-26b-5p may also be involved in T2DM by potentially regulating the expression of various candidate genes, such as CTGF, PCK1, PPP1RA, PCOLCE2, FAM105A, TRAF3PI2 and neuraminidase 3.

A TF-target gene regulatory network was constructed, from which 10 TFs were identified and Smad5 was a potential target for T2DM treatment (105). The Forkhead box class O/Bcl6/cyclin D2 pathway connects nutrient and growth factor status to cell cycle control in pancreatic β-cells, and should therefore be considered for its therapeutic potential in diabetes (106). Notably, 8 of these transcription regulatory factors, SUZ12, EZH2, ZNF580, KLF9, MAZ, ATF1, SSRP1, WRNIP1, CHD1 were shown to be involved in the development of T2DM by modulating the expression of various candidate genes such as ankyrin repeat domain 23 (ANKRD23), transmembrane protein 37 (TMEM37), PPP1R1A, PCK1, CTGF, ISG1, SSRP1, WRNIP1 and CHD1, which have not been previously reported to be dysregulated in T2DM. The present study predicted that these TFs might play key roles in the occurrence and development of T2DM. This result provided preliminary evidence that a lower expression of TMEM37 could reflect a decrease in β-cell numbers in T2DM (107). ANKRD23, a diabetes-related ankyrin repeat protein, was identified as a novel gene that is upregulated in the hearts of KKA(y) mice, a T2DM and insulin-resistant animal model (108). Solimena et al (110) also identified that ANKRD23, PPP1R1A and TMEM37 were enriched in β-cells and downregulated in T2DM. TMEM37 prohibits Ca2+-influx and insulin secretion in β-cells (109).

The present study has some limitations. The number of samples was relatively inadequate, although the combining of multiple datasets can compensate for missing or unreliable information in any single dataset. Additionally, the present results are preliminary and descriptive. Integrative analysis of gene profiling data cannot entirely exclude false positive results. Furthermore, the present study only discussed mRNA expression and did not refer to the protein expression of the factors identified. Due to post-transcription regulatory events, protein expression levels may or may not correlate with mRNA levels. However, the alteration of protein structure, function and interaction is the underlying mechanism of a number of diseases, including diabetes (110–112). Therefore, some experiments, such as reverse transcription-quantitative PCR (113,114), western blotting (115), cross-linking immunoprecipitation (116) or functional experimental validation, are needed to validate key genes, TFs and miRNAs and relevant proteins in T2DM development. Despite these limitations, the present study still provided insight for understanding the complicated underlying molecular mechanisms of T2DM.

Overall, the bioinformatics analysis of the present study identified potential markers that may play a potential role in the occurrence, development and treatment of T2DM. A total of 83 candidate genes were selected, and ISG15, PCK1, SYT1, IGFBP3, TIMM44 and TRAF3IP2 could be the core genes of T2DM. Certain key BPs such as ‘inflammatory response’, ‘cellular responses to cytokine stimulus’, ‘lipid metabolic process’, ‘positive regulation of cell death’ and ‘positive regulation of cell proliferation’, and certain signaling pathways associated with the PI3K-Akt, MAPK, Rap1 and Ras signaling pathways were identified to be involved in T2DM. The present study also identified miRNAs, including hsa-miR-192-5p, hsa-miR-124-5p and hsa-miR-335-5p, and TFs, including Smad5 and Bcl6, that might be potential targets for the diagnosis and treatment of T2DM. In addition, has-miR-8485, has-miR-1277-3p, has-miR-190a-3p, has-miR-5011-3p, has-let-7a-5p, has-let-7b-5p, has-miR-98-5p, has-miR-7106-5p and has-miR-26b-5p, and TFs SUZ12, EZH2, ZNF580, KLF9, MAZ, ATF1, SSRP1, WRNIP1 and CHD1 have not been previously identified to be related to T2DM, to the best of the authors’ knowledge, while they and their target genes may serve as diagnostic indicators for patients with T2DM. To obtain more reliable correlation results, it is necessary to validate the predicted results with a series of verification experiments. The present study identified candidate genes for T2DM development, which might be redefined as pathogenic genes for T2DM diagnosis and therapy. The experimental results could provide insight for future genomic individualized treatment of T2DM and help to identify the underlying molecular mechanisms that lead to T2DM.