In silico prediction and validation of potential therapeutic genes in pancreatic β-cells associated with type 2 diabetes.

Diabetes mellitus is becoming a major health burden worldwide. Pancreatic β-cell death is a characteristic of type 2 diabetes (T2D), but the underlying mechanisms of pancreatic β-cell death remain unknown. Therefore, the aim of the present study was to identify potential targets in the pancreatic islet of T2D. The GSE20966 dataset was obtained from the Gene Expression Omnibus (GEO) database, and differentially expressed genes (DEGs) were identified by using the GEO2R tool. The Gene Ontology terms and Kyoto Encyclopedia of Genes and Genomes Pathway enrichment analysis of DEGs were further assessed using the Database for Annotation, Visualization and Integrated Discovery. Furthermore, protein-protein interaction (PPI) networks were constructed for the up- and downregulated genes using STRING databases and were then visualized with Cytoscape. The body weight, fasting blood glucose (FBG), pancreatic index and biochemistry parameters were measured in db/db mice. Moreover, the morphology of the pancreas was detected by hematoxylin and eosin staining, and hub genes were assessed using reverse transcription-quantitative PCR (RT-qPCR) and western blot analysis. In total, 570 DEGs were screened, including 376 upregulated and 194 downregulated genes, which were associated with 'complement activation, classical pathway', 'proteolysis', 'complement activation' and 'pancreatic secretion pathway'. It was found that the body weight, FBG, alanine aminotransferase, aspartate aminotransferase, total cholesterol, triglycerides, blood urea nitrogen, creatinine, fasting serum insulin, glucagon and low-density lipoprotein cholesterol levels were significantly higher in db/db mice, while high-density lipoprotein cholesterol levels and the pancreatic index were significantly decreased. Furthermore, albumin, interleukin-8, CD44, C-C motif chemokine ligand 2, hepatocyte growth factor, cystic fibrosis transmembrane conductance regulator, histone cluster 1 H2B family member n, mitogen-activated protein kinase 11 and neurotrophic receptor tyrosine kinase 2 were identified as hub genes in PPI network. RT-qPCR and western blotting results demonstrated the same expression trend in hub genes as found by the bioinformatics analysis. Therefore, the present study identified a series of hub genes involved in the progression of pancreatic β-cell, which may help to develop effective therapeutic strategy for T2D. Copyright: © Zhou et al.

Introduction

Diabetes mellitus (DM), a common syndrome of disordered metabolism, is a growing major public health issue that involves the deficiency of insulin production by β-cells in type 1 diabetes (T1D), but approximately 95% of patients have T2D (1). T2D is characterized by peripheral insulin resistance (IR) and decompensation of the pancreatic β-cells (1,2). With rapidly increasing incidence and prevalence rates, 415 million individuals have been diagnosed with T2D in 2015 globally, and this figure is expected to rise to 642 million individuals by 2040 (2-4). Moreover, previous studies have shown that the association between several genetic and environmental factors contribute to the risk of developing T2D, causing IR and the dysfunction of pancreatic β-cells (4,5).

The pancreas is an important endocrine organ, and pancreatic β-cells play a vital role in regulating the stability and metabolism of blood glucose, and are the only class of cells characterized by synthesizing and secreting insulin in the body (6-8). Furthermore, progressive pancreatic β-cells dysfunction and apoptosis are recognized as fundamental pathologies of T2D (9,10). Previous studies have revealed that the abnormal transcription of lncRNAs growth arrest specific 5 (Gas5), MAF bZIP transcription factor A (MAFA), MAF bZIP transcription factor B (MAFB), NK6 homeobox 1 (NKX6.1), and pancreatic and duodenal homeobox 1 (PDX1), may be associated with islet functional injury (7,11). Additionally, Robertson et al (12) reported that deficiency in PDX1, a central administrator of insulin promoter activity, is an essential mechanism resulting in β-cell dysfunction. Moreover, Okauchi et al (13) showed that suppressing sodium-glucose co-transporter 2 using luseogliflozin has protective effects on pancreatic β-cell mass and function to improve glycemic control, due to the increase in β-cell proliferation and reduction of β-cell apoptosis. In addition, reactive oxygen species are known to enhance NF-κB activity, which potentially induces β-cell apoptosis (14). Therefore, these results indicate that transcription factors are involved in maintaining the development processes of adult β-cells. Significant advances have been made in understanding glucose homeostasis and the pathophysiology of T2D over the past few decades (4,15); however, the identification of islet-specific molecular targets and their mechanisms of pathology in T2D remains unknown.

Previous studies using high-throughput microarray platforms have facilitated knowledge discovery and are increasingly valued for identifying the initiation and progression of complex diseases from multiple scales (16,17). With the continuous progression of the systems biology research field, network pharmacology is becoming a novel subject for scientific researchers to investigate and identify interconnecting disease features, bioactive agents and drug targets in complex biological systems (16-18), as well as bioinformatics, an interdisciplinary field that has key roles in the analysis of high-throughput data (19). Microarray analysis of T2D, combining network pharmacology with bioinformatics methods, provides a promising strategy to identify molecular biological information, and are also powerful tools for the development of anti-T2D targets (19). A large amount of core data on complex diseases, such as T2D, has been produced from gene combining chromatin immunoprecipitation, which are uploaded on the Gene Expression Omnibus (GEO) public database of the National Center of Biotechnology Information (20). Furthermore, extracting and re-analyzing molecular data can provide a valuable perspective on understanding the molecular mechanisms of complex diseases (20).

The aims of the present study were to identify and assess the transcriptional expression levels of differently expressed genes (DEGs) in pancreatic β-cells associated with T2D. Network pharmacology and bioinformatics re-analysis were performed on original data downloaded from GEO (GSE20966) of β-cells separated from the pancreases of ten healthy controls and ten T2D samples, from a study by Marselli et al (21). The validation of key candidate DEGs expressions in pancreatic β-cells, which were isolated from C57BLK/J-db/db and db/m mice, was performed using reverse transcription-quantitative PCR (RT-qPCR). The present results provided further understanding of the pathogenesis of T2D at the molecular level and identified potential targets for individualized treatment of DM.

Materials and methods

Data collection

Gene expression profiling of GSE20966, using the Agilent-014850 Whole Human Genome Microarray 4x44 K G4112F platform (Agilent Technologies, Inc.), deposited by Marselli et al (21) was downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). A total of 20 samples (GSM524151-GSM524170) including 10 controls and 10 T2D subjects, along with a laser capture microdissection technique (21), were used to extract β-cells from pancreatic tissue sections of the samples.

Data processing and identification of DEGs

GEO2R (Version: 2.40.0; http://www.ncbi.nlm.nih.gov/geo/geo2r/), an interactive web-based tool, was used to identify DEGs between T2D and healthy pancreatic β-cells. Furthermore, the GEO2R tool was used to view the distribution of the expression values in the samples, which were graphically demonstrated as a box plot. Subsequently, the linear models for microarray data package (limma; Version: 3.26.8; https://www.bioconductor.org/) in R/Bioconductor (Version: R 3.2.3; https://www.r-project.org/) were used to screen the DEGs between T2D and healthy pancreatic β-cells, with P<0.05 and log2|(FC)|≥1 used as the cut-off criteria. Statistically significant DEGs were visualized on a heatmap.

Functional enrichment analysis of DEGs

Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/), an online biological information database, was used to annotate specific biological attributes for large-scale molecular datasets using high-throughput microarrays (22,23). The KEGG pathway database (https://www.kegg.jp/) was used to perform GO functional and KEGG pathway enrichment analyses for DEGs, via DAVID. P<0.05 and Benjamini-value of <0.05 were considered as statistically significant.

Establishment of PPI network and hub DEGs screening

The STRING database (https://string-db.org/) is an online tool, which provides both experimental and predicted interaction information for evaluating PPI networks (24). Therefore, the interaction of the upregulated and downregulated DEGs were analyzed by STRING, and a confidence with interaction score of ≥0.4 was used as the cut-off criteria (25,26). Subsequently, complex networks of the upregulated and downregulated DEGs were visualized using Cytoscape (version 3.2.1; http://cytoscape.org/). The hub genes in the up- and down-regulated DEGs PPI networks were identified and the network topology analysis was performed using the Network Analyzer plug-in in the Cytoscape (v3.6.1; https://cytoscape.org/). The hub genes, which may play a critical role in the network, were identified by using the front value of node degree, node betweenness and closeness of nodes (25).

Animals

The db/db male mice (weight, 31-36 g; age, 7 weeks; n=8) and male non-diabetic control (db/m) mice (weight, 17-26 g; age, 7 weeks; n=8) were purchased from the Nanjing Biological Medicine Research Institute Affiliated with Nanjing University. All the mice were housed in cages at 25˚C with 50-70% humidity and maintained on a 12-h light/darkness cycle, with free access to water and food in the Animal Center of Zhejiang Chinese Medicine University. All animal experimental procedures were based on the guidelines from, and were approved by, the Animal Care and Welfare Committee of Zhejiang Chinese Medical University (permit no. ZSLL-2018-196). db/db mice and db/m mice were fed with basic feed for 7 weeks. At the 14th week, the body weight, fasting blood glucose (FBG) and biochemical parameters were measured. At the defined time points, mice were euthanized using sodium pentobarbital (150 mg/kg, intraperitoneal) after an 8-h fast. The eyeballs were removed, and then 1.5 ml blood samples were collected. When the heart of the mouse stopped completely, pancreas tissues were removed, weighed and the index of the pancreas to body weight was calculated. A part of pancreas was frozen immediately in liquid nitrogen and stored at -80˚C until further analysis, and the remaining tissue sections were stored in formalin at room temperature, and the islet were isolated. Serum samples were acquired by centrifugation at 1,000 x g for 10 min at 4˚C for the measurement of alanine aminotransferase (ALT), aspartate aminotransferase (AST), total cholesterol (TC), triglycerides (TG), blood urea nitrogen (BUN), creatinine (CREA), low-density lipoprotein cholesterol (LDL-C) and high-density lipoprotein cholesterol (HDL-C) levels using a 7020 fully automatic biochemical analyzer (Hitachi Ltd.). In addition, the fasting serum insulin (FINS) and glucagon levels were measured by radioimmunoassay (RIA) method using a commercial kit (cat. no. 2013000920; Shanghai Ruiqi Biotechnology Co., Ltd.).

Morphology

Samples from the pancreas were fixed in 10% buffered formalin for 24 h at room temperature, embedded in paraffin and cut into 6-µm sections. The sections were subsequently stained with hematoxylin and eosin (H&E staining kit; cat. no. G1120; Beijing Solarbio Science & Technology Co., Ltd.) at 37˚C for 20 min. Pancreatic lesions were examined under a light microscope (magnification, x200; DMI 3000B; Leica Microsystems GmbH).

Isolation of pancreatic islets

Pancreatic islets were isolated as previously described (27). The pancreas was excised and digested by 0.5 mg/l collagenase P (Roche Diagnostics, Basel, Switzerland) for 13 min at 37˚C and acquired using a density gradient. Islet cells at a density of 5x104/ml were cultured in six-well plates with RPMI-1640 medium (Gibco; Thermo Fisher Scientific, Inc.) containing 11.1 mM glucose and 10 mM HEPES, supplemented with 10% FBS (Gibco; Thermo Fisher Scientific, Inc.), 100 mg/ml streptomycin and 100 U/ml penicillin for 12 h in a 95% O2/5% CO2 atmosphere to recover from the isolation. Islet cells were then harvested to investigate the transcriptional expression levels of hub genes.

RT-qPCR

Total RNAs were extracted from isolated mouse islets using TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. The quality and quantity of RNA samples were measured using spectrophotometry at 260/280 nm (Thermo Fisher Scientific, Inc.), and cDNA was obtained using a PrimeScript RT Master Mix kit (cat. no. RR036A; Takara Biotechnology Co., Ltd.) according to the manufacturer's instructions. qPCR was carried out using a CFX96 Touch Deep Well detection system (Bio-Rad Laboratories, Inc.). Each reaction (10 µl) contained 1 µl cDNA template, 0.4 µl forward and reverse primers, 3.2 µl RNase-free ddH2O and 5 µl SYBR-Green Master Mix (Takara Biotechnology Co., Ltd.). The following thermocycling conditions were used for the qPCR: Initial denaturation at 95˚C for 30 sec; 40 cycles of denaturation at 95˚C for 5 sec and annealing at 60˚C for 30 sec and dissolution curve conditions at 65˚C for 0.05 sec and 95˚C for 0.5 sec. The forward and reverse primer sequences used for PCR are shown in Table I. The transcriptional expression levels of the hub genes were normalized to GAPDH by comparing the cycle threshold (Cq) values, calculated by the 2-ΔΔCq method (28).

Table I

Primers used for reverse transcription-quantitative PCR.

Gene	Forward primer (5'→3')	Reverse primer (5'→3')
ALB	TGGTGAAATGGCTGACTGCT	CTCTGGTCTCACCAATCGGG
IL-8	GGACAACAGAGAGGTGTGCT	ACAGCGGTGCATCAGAATTG
CD44	CCCATTCGACAACAGGGACA	TGGGGTGTGAGATTGGGTTG
CCL2	GATGCAATCAATGCCCCAGTC	TTTGGGACACTTGCTGCTGG
HGF	ACCCTGGTGTTTCACAAGCA	GCAAGAATTTGTGCCGGTGT
CFTR	CCCAGCCATTTTTGGCCTTC	GACGCCTGTAACAACTCCCA
HIST1H2BN	CACCGGTATCTCGTCCAAGG	ATGGTCGAGCGCTTGTTGTA
MAPK11	CGACGAGCACGTTCAATTCC	TCACAGTCCTCGTTCACAGC
NTRK2	ACTGTGAAAGGCAACCCCAA	CAGCCCATGAAGTGAGCAGA
GAPDH	ATTCCATGGCACCGTCAAGG	TCGCCCCACTTGATTTTGGA

ALB, albumin; IL, interleukin; CCL2, C-C motif chemokine ligand 2; HGF, hepatocyte growth factor; CFTR, cystic fibrosis transmembrane conductance regulator; HIST1H2BN, histone cluster 1 H2B family member n; MAPK11, mitogen-activated protein kinase 11; NTRK2, neurotrophic receptor tyrosine kinase 2.

Western blot analysis

Isolated mouse islets were lysed by RIPA buffer (Beyotime Institute of Biotechnology) to obtain the total protein. A bicinchoninic acid protein assay kit (cat. no. P0011; Beyotime Institute of Biotechnology) was used to measure the concentration of the protein samples. The samples were boiled for 5 min, and 30 ng protein from each cell lysate were resolved by SDS-PAGE on 8 or 12% gels and transferred to a PVDF membrane (EMD Millipore). Then, the membranes were blocked with 5% skim milk at 37˚C for 1 h and incubated at 4˚C overnight with the primary antibodies against: Albumin (ALB; 1:2,000; cat. no. ab207327; Abcam), interleukin (IL)-8 (1:1,000; cat. no. sc-376750; Santa Cruz Biotechnology, Inc.), CD44 (1:200; cat. no. sc-7297; Santa Biotechnology, Inc.), C-C motif chemokine ligand 2 (CCL-2; 1:1,000; cat. no. ab9669; Abcam), hepatocyte growth factor (HGF; 1:2,000; cat. no. ER1910-54; HuaBio, Inc.), cystic fibrosis transmembrane conductance regulator (CFTR; 1:500; cat. no. ab2784; Abcam), histone cluster 1 H2B family member n (HIST1H2BN; 1:100; cat. no. sc-515808; Santa Biotechnology, Inc.), mitogen-activated protein kinase 11 (MAPK11; 1:1,000; cat. no. ab183208; Abcam), neurotrophic receptor tyrosine kinase 2 (NTRK2; 1:1,000; cat. no. ET1701-16; HuaBio, Inc.) and GAPDH (1:5,000; cat. no. 60004-1-1 g; HuaBio, Inc.), which served as a loading control. Subsequently, the membranes were washed with TBS containing 0.1% Tween-20, and incubated with horseradish peroxidase-conjugated secondary antibodies (cat. nos. A9294 and A6715; 1:5,000; Sigma-Aldrich; Merck KGaA) at room temperature for 1 h. To visualize the specific protein bands, an ECL advanced western blot detection kit (Thermo Fisher Scientific, Inc.) was used, and the intensity of bands was measured using ImageJ software (version 1.46; National Institutes of Health).

Statistical analysis

Statistical analysis was performed using SPSS version 21.0 software (IBM Corp.). Data are presented as the mean ± SD. Student's t-test was used to analyze differences between db/db mice and control groups. P<0.05 was considered to indicate a statistically significant difference.

Results

Data processing and DEG screening

As shown in Fig. 1A, raw expression data were normalized after preprocessing; median-centered values demonstrated that the data were normalized and thus it was possible to cross-compare between control and T2D samples. Using the threshold of P<0.05 and log2|FC|≥1, a total of 570 genes were screened to be differentially expressed in GSE20966 (Fig. 1B). Among them, 376 were upregulated and 194 were downregulated DEGs. Moreover, the top 40 upregulated and 40 downregulated DEGs were analyzed by heatmap clustering based on Euclidean distance (Fig. 1C).

Figure 1

Distribution of differentially expressed genes in the control and type 2 diabetic samples. (A) Box plots of the GSE20966 datasets following normalization. The x-axis represents the samples, and the y-axis represents the distribution of expression levels. (B) Volcano plot of the up- and downregulated DEGs. The horizontal axis is logFC and the vertical axis is the P-value after the negative logarithm conversion. Up- and downregulated genes are marked in red and green, respectively. (C) Heat map of the upregulated and downregulated genes. The horizontal axis represents the sample and the vertical axis represents the genes. The red and green blocks represent the expression pattern of the gene; green represents low expression, and red represents high expression.

Enrichment analysis of DEGs

DAVID was used to investigate GO function and KEGG pathway enrichment analysis of upregulated and downregulated DEGs. Following GO function analysis, the DEGs were divided into biological process, molecular function and cellular component terms. As shown in Table II, it was found that the majority of enriched GO terms for biological process were ‘complement activation, classical pathway’, ‘proteolysis’ and ‘complement activation’. Furthermore, the most enriched GO terms of cellular component were the ‘extracellular region’, ‘extracellular space’, ‘blood microparticle’ and ‘extracellular exosome’. It was also identified that the enrichment for molecular function was for ‘serine-type endopeptidase activity’. Following the KEGG pathway analysis, it was found that the 570 DEGs were primarily involved in ‘pancreatic secretion’, ‘complement and coagulation cascades’, ‘protein digestion and absorption’ and ‘TNF signaling pathway’, with a P-value and Benjamini-value of <0.05.

Table II

Gene Ontology and pathway analysis of the differentially expressed genes in GSE20966.

GO	GO ID	GO term	P-value	Benjamini	Count
GO_BP	GO:0006958	Complement activation, classical pathway	4.49x10-7	9.57x10-4	15
GO_BP	GO:0006508	Proteolysis	2.57x10-6	2.74x10-3	34
GO_BP	GO:0006956	Complement activation	2.29x10-5	1.62x10-2	12
GO_CC	GO:0005615	Extracellular space	9.22x10-14	3.16x10-11	86
GO_CC	GO:0005576	Extracellular region	4.44x10-10	7.61x10-8	87
GO_CC	GO:0072562	Blood microparticle	1.46x10-7	1.67x10-5	19
GO_CC	GO:0070062	Extracellular exosome	4.42x10-6	3.79x10-4	114
GO_MF	GO:0004252	Serine-type endopeptidase activity	4.37x10-7	2.87x10-4	24
PATHWAY	hsa04972	Pancreatic secretion	8.40x10-8	1.97x10-5	16
PATHWAY	hsa04610	Complement and coagulation cascades	5.31x10-6	6.24x10-4	12
PATHWAY	hsa04974	Protein digestion and absorption	2.75x10-4	0.021326216	11
PATHWAY	hsa04668	TNF signaling pathway	3.07x10-4	0.017869834	12
PATHWAY	hsa04978	Mineral absorption	3.88x10-4	0.018059595	8

GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function; TNF, tumor necrosis factor.

PPI network analysis of upregulated and downregulated DEGs

The PPI network of upregulated genes is shown in Fig. 2, and was constructed with 198 nodes and 451 edges. Moreover, three topological properties, ‘degree’, ‘betweenness’ and ‘closeness’ were calculated to screen for the putative hub genes, and six central node genes were identified as putative hub genes, by filtering of degree >20 criteria in the network (red in Fig. 2): ALB (degree=61, node betweenness=0.50283976 and closeness=0.5327381), IL-8 (degree=33, node betweenness=0.14382652 and closeness=0.45316456), CD44 (degree=27, node betweenness=0.13017011 and closeness=0.42316785), CCL2 (degree=23, node betweenness=0.05333124 and closeness=0.42417062), HGF (degree=21, node betweenness=0.0329843 and closeness=0.4124424) and CFTR (degree=21, node betweenness=0.09656952 and closeness=0.40406321).

Figure 2

Protein-protein interaction network of upregulated DEGs. In total, six upregulated DEGs in red were identified as putative candidate genes. DEGs of low degree values are presented as small-sized nodes of bright colors. DEGs, differentially expressed genes.

The downregulated PPI network consisted of 58 nodes and 47 edges (Fig. 3). In addition, three down-regulated genes (in green) were selected as putative hub genes in the downregulated PPI networks: HIST1H2BN (degree=4, node betweenness=0.38562092 and closeness=0.31034483), MAPK11 (degree=4, node betweenness=0.60130719 and closeness=0.36) and NTRK2 (degree=4, node betweenness=0.21568627 and closeness=0.29508197.

Figure 3

Protein-protein interaction network of downregulated DEGs. In total, three downregulated DEGs in green were identified as putative candidate genes. DEGs of low degree values are presented as small-sized nodes of bright colors. DEGs, differentially expressed genes.

Measuring the differences between db/db and control group mice

A number of characteristic differences were measured between db/db and control mice at the 14th week (Fig. S1). The results indicated that the body weights of db/db mice were significantly increased compared with the control group. Moreover, the FBG value was significantly upregulated in the db/db mice group. The biochemical parameters of ALT, AST, TG, TC, BUN, CREA, LDL-C, HDL-C, FINS and glucagon were measured from sera collected at the experimental end point. The results indicated that the db/db mice group had higher levels of ALT, AST, TG, TC, BUN, CREA, FINS, glucagon and LDL-C compared with the control group (P<0.01). However, the levels of HDL-C and pancreatic/body weight coefficient were significantly decreased in db/db mice.

Islet morphology changes in db/db mice

Light microscopy examination of HE stained sections identified the appearance of pancreatic acinar structures in the control mice group (Fig. 4A). Moreover, pancreatic islets of healthy pancreatic tissue were round or oval and had a clear boundary, with the cells arranged in a neat, medium-sized uniform distribution. However, in the db/db mice group, pancreatic morphologies were abnormal and the boundaries were not clear, and the cell number was decreased (Fig. 4B).

Figure 4

Islet morphology changes in db/db mice. (A) Control mice group. (B) db/db mice group. Magnification, x200.

Hub genes expression

The islets of pancreas were removed to detect the mRNA expression levels of hub genes in the upregulated and downregulated PPI network, including ALB, IL-8, CD44, CCL2, HGF, CFTR, HIST1H2BN, MAPK11 and NTRK2. It was found that the mRNA expression levels of ALB, IL8, CD44, CCL2, HGF and CFTR were significantly increased in the db/db mice group compared with the control group (P<0.05; Fig. 5A-F). Furthermore, the results indicate that the mRNA expression levels of HIST1H2BN, MAPK11 and NTRK2 were significantly decreased compared with the control group (P<0.05; Fig. 5G-I). Therefore, the RT-qPCR results of the DEGs were consistent with the microarray analysis.

Figure 5

Validation of microarray data with RT-qPCR. Microarray and RT-qPCR for (A) ALB, (B) IL-8, (C) CD44, (D) CCL2, (E) HGF, (F) CFTR, (G) HIST1H2BN, (H) MAPK11 and (I) NTRK2. RT-qPCR results were consistent with the microarray data, and the P-values were calculated by Student's t-test. ▲P<0.05, ▲▲P<0.01. RT-qPCR, reverse transcription-quantitative PCR; ALB, albumin; IL, interleukin; CCL2, C-C motif chemokine ligand 2; HGF, hepatocyte growth factor; CFTR, cystic fibrosis transmembrane conductance regulator; HIST1H2BN, histone cluster 1 H2B family member n; MAPK11, mitogen-activated protein kinase 11; NTRK2, neurotrophic receptor tyrosine kinase 2.

In addition, it was found that the protein expression levels of ALB and IL-8 were significantly higher in db/db mice compared with the control group (P<0.01; Fig. 6A-C). However, there was no significant difference in CD44 and CCL2 protein expression levels (Fig. 6D and E). Moreover, the results suggested that db/db mice had significantly higher HGF and CFTR protein expression levels compared with control mice (P<0.01 and P<0.05, respectively; Fig. 6F and G). It was also found that there were no significant differences in the expression of HIST1H2BN between db/db and control mice groups (Fig. 6H). In addition, the results indicated that the protein expression levels of MAPK11 and NTRK2 were significantly lower in db/db mice compared with the control group (P<0.05 and P<0.01, respectively; Fig. 6I and J).

Figure 6

Western blotting. (A) Detection of protein expression levels of (B) ALB, (C) IL-8, (D) CD44, (E) CCL-2, (F) HGF, (G) CFTR, (H) HIST1H2BN, (I) MAPK11 and (J) NTRK2. *P<0.05, **P<0.01 vs. control mice. ALB, albumin; IL, interleukin; CCL2, C-C motif chemokine ligand 2; HGF, hepatocyte growth factor; CFTR, cystic fibrosis transmembrane conductance regulator; HIST1H2BN, histone cluster 1 H2B family member n; MAPK11, mitogen-activated protein kinase 11; NTRK2, neurotrophic receptor tyrosine kinase 2; Con, control.

Discussion

T2D, the most common form of DM worldwide, is a complex metabolic disease and is characterized by varied degrees of dysfunction in pancreatic β-cells, resulting in IR (29,30). However, the pathophysiology of T2D has not been as extensively studied, with regards to defective insulin secretion and impaired insulin action in peripheral tissues. In the present study, putative hub genes and enriched pathways of pancreatic β-cells were identified using bioinformatics analysis, and assessed using RT-qPCR and western blot analysis. Gene expression data from GSE20966 were re-analyzed, and it was found that there were 376 upregulated and 194 downregulated DEGs in patients with diabetes compared with the control subjects. The functional annotation analysis results suggested that these DEGs participated in ‘complement activation and classical pathways’. Furthermore, abnormal insulin secretion was one of the prominent molecular events of pancreas. It was also demonstrated that the ‘pancreatic secretion pathway’ was the most significantly enriched KEGG pathway.

In total, nine hub genes were identified, which were screened from the calculation of network topology in the upregulated and downregulated PPI networks. In the present study, the body weight, FBG and biochemical parameters for the db/db mice were significantly increased compared with the control mice. The leptin receptor-deficient db/db mouse model is a generally-used research model, in which the development of progressive diabetes is characterized by reduced proliferation and increased apoptosis in pancreatic β-cells (31). Similarly, the pancreatic/body weight coefficient and islet morphology changes identified in the present study were consistent with previous studies (31,32). Moreover, using the GSE20966 dataset, Ding et al (33) revealed DEGs that are enriched in the pancreatic secretion pathway, and that serpin family G member 1 and alanyl aminopeptidase membrane are associated with the development of T2D. In addition, Zhong et al (34) reported that ATP-citrate lyase is hub gene involved in the molecular mechanism of T2D. Based on previous studies, the characteristic of the PPI network were further analyzed in the present study, which identified additional hub genes, including ALB, IL-8, CD44, CCL2, HGF, CFTR, HIST1H2BN, MAPK11 and NTRK2. It was found that these hub genes exhibited the same expression levels as the microarray data of GSE20966, and the protein expression levels of ALB, IL-8, CD44, HGF, CFTR, MAPK11 and NTRK2 were consistent with the RT-qPCR results.

ALB is a useful marker of nutrition and general health status, and it has been shown that ALB concentration is associated with metabolic syndromes, diabetes, coronary heart disease, all-cause mortality and IR (35-37). T2D is also characterized by systemic low-grade inflammation, which promotes the development of IR (38). IL-8, a pro-inflammatory polypeptide, is one of the most central members of the CXC chemokine superfamily, and is involved in macrophage infiltration, systemic immunity and activation of adipose tissue (38). Furthermore, it has been reported that IL-8 may play a major role in the pathogenesis of T2D (38). Previous studies revealed that inflammatory stimulation induces IL-8 production in human adipocytes via the p38 MAPK pathway and/or ERK pathway, and can induce IR via the inhibition of Akt phosphorylation in human adipocytes (38). Moreover, Zozuliñska et al (39) showed that increased IL-8 expression in a small cohort of both patients with T1D or T2D results in poor glycemic control compared with healthy control subjects (38,39).

CD44 is a multifunctional adhesion molecule receptor expressed at the cell surface of extracellular matrix (ECM) components, which contains 20 exons and generates multiple protein isoforms, resulting in standard (CD44s) and variant (CD44v) isoforms, via alternative transcript splicing (40,41). Previous studies have shown CD44 expression to be positively associated with adipose inflammation and the index of IR in adipose tissue, and that serum CD44s levels are positively associated with IR and glycemic control in patients with diabetes (42,43). Kobayashi et al (41) reported that CD44v is highly expressed in pancreatic islets of diabetic mice and that it inhibits insulin secretion by attenuating amino acid transport mediated by the L-type amino acid transporter LAT1. Similar to these previous studies, the present study also found upregulation of CD44 in db/db mice, indicating that CD44 may be a promising therapeutic target for DM.

The pro-inflammatory chemokine CCL2 is a member of the CC chemokine family, which contains 76 amino acids, and plays a key role in the inflammatory process (44). Studies on the intervention of CCL2 or CCL2-deficient diabetic mouse kidneys revealed that a lack of CCL2 results in reduced accumulation of interstitial myofibroblasts, renal production of transforming growth factor-β1 and deposition of ECM (45,46). As reported by Hopfgarten et al (47), inflammation of the islet promotes the loss of β-cells and T2D, as shown by increased expression CCL2 in the islets of donors with T1D, whereas it is undetected in healthy donors. The results of the present study is consistent with the findings of Hopfgarten et al (47), indicating that CCL2 is a crucial potential therapeutic target for the treatment of inflammation in T2D.

HGF is a mesenchymal-derived pleiotropic factor, which has multiple roles, including cellular proliferation, morphogenesis, motogenic, anti-apoptotic, angiogenic and angioprotective, as well as tissue regeneration and repair (48-50). Mellado-Gil et al (51) showed that the disruption of the endogenous pancreatic HGF/c-Met signaling pathway enhances pancreatic β-cell death and accelerates the onset of DM in diabetogenic conditions. Moreover, it has been revealed that there is increased HGF protein expression in adipose and liver tissues with IR (51). Araújo et al (52) also reported that HGF treatment on isolated β-cells from rats results in increased levels of insulin receptor substrate 2 tyrosine phosphorylation, as well as AKTser473 and ERK phosphorylation, thus suggesting that the correlation between HGF and β-cell mass increase may involve these proteins from insulin signaling pathways. Consistent with previous studies, the present results suggested that the expression of HGF in the islet of db/db mice was higher compared with the control mice.

Cystic fibrosis-related (CF-related) diabetes (CFRD) is identified as a form of type 3c (pancreatogenic) DM, as CFRD shows characteristics of both T1D and T2D, with regards to the occurrence and development of chronic infection, pancreatic insufficiency and varying degrees of IR (53). Furthermore, CFRD is caused by mutations in CFTR (53,54), which is a Cl- channel that belongs to the ATP-binding cassette-transporter proteins family and is activated by cAMP. It has been shown that mutations in the CFTR impair the function of ion channels, causing cystic fibrosis (55). The present study found that CFTR levels were significantly increased in db/db mice. These data suggested that abnormal expression and mutations of CFTR may be a potential molecule for therapeutic response of DM.

HIST1H2BN encodes for members of the histone H2B family in humans located on human chromosome 6 (6p21-p22), regarded as the histone H2B-type 1-N protein, and also functions in the compaction of chromatin into higher-order structures (56). In the present study, it was demonstrated that the transcriptional expression level of HIST1H2BN was downregulated in GSE20966; however, additional prospective experiments are required to investigate whether HIST1H2BN acts as a therapeutic target in T2D.

P38 MAPK belongs to the MAPK superfamily, which is composed of four isoforms: MAPK14, MAPK11, MAPK12 and MAPK13. Among them, MAPK14 and MAPK11 are the predominant forms, and play important roles in the occurrence and development of T2D, due to increased apoptosis of pancreatic β-cells (57,58). Huang et al (58) reported that diazoxide prevents the onset and development of DM in OLETF rats by inhibiting β-cell apoptosis via increasing MAPK11, and elevating the Bcl-2/Bax ratio to ameliorate insulin secretory capacity. In addition, the mRNA and protein expression of MAPK11 were significantly lower in the db/db mice of the present study. The result indicated that MAPK11 might be a novel therapeutic strategy to prevent and manage T2D. To the best of our knowledge, this is the first report in elucidating the abnormality of MAPK11 in db/db mice, and further research of MAPK 11 on β-cell function is required.

NTRK2, a possible candidate gene on 9q21, belongs to the family of NTRK proteins and is the membrane-bound receptor for brain-derived neurotrophic factor (BDNF) (59). In addition, BDNF binding triggers the signaling cascades of phospholipase Cγ, MAPK and PI3K pathways with tyrosine phosphorylation in the intracellular domains of NTRK2(58). Furthermore, BDNF/NTRK2 deficiency and mutations of NTRK2 are correlated with increased weight and obesity in humans and mice (59-61). In the present study, it was identified that NTRK2 had significantly suppressed expression in db/db pancreatic β-cells.

Therefore, the aforementioned results suggested that these hub genes may be potential therapeutic targets and may play a critical role in pancreatic β-cell death of T2D; however, this process is currently not fully understood. Moreover, it has been revealed that body mass index (BMI) is the most common predictor to identify IR (62); however, the relationship between BMI and IR was not investigated in the present study, and only body weight was measured. In addition, the sample sizes were relatively few in the present study, and thus should be performed with clinical cases. Therefore, future comprehensive studies are required with larger clinical samples to identify the mechanisms related to the potential targets to increase understand their roles in T2D.

In conclusion, the present results suggested ALB, IL-8, CD44, CCL2, HGF, CFTR, HIST1H2BN, MAPK11 and NTRK2 may be potential core genes of T2D, which may facilitate the development of treatment options for T2D. However, further studies are required to investigate the significant pancreatic β-cell-associated genes that have potential for targeted therapy of T2D.