Bioinformatics analysis of hepatic gene expression profiles in type 2 diabetes mellitus.

Type 2 diabetes mellitus (T2DM) is characterized by hyperglycemia. The liver has a critical role in regulating glucose homeostasis. The present study aimed to analyze hepatic gene expression profiles and to identify the key genes and pathways involved in T2DM. Gene expression profiles of 10 patients with T2DM and 7 subjects with normal glucose tolerance were downloaded from the Gene Expression Omnibus database. Subsequently, differentially expressed genes (DEGs) were identified and functional enrichment analysis was performed. In addition, a protein-protein interaction network was built and hub genes were identified. In total, 1,320 DEGs were identified, including 698 up- and 622 downregulated genes, and these were mainly enriched in positive regulation of transcription from RNA polymerase II promoter, cell adhesion, inflammatory response, positive regulation of apoptotic process, signal transduction and the Tolllike receptor signaling pathway. A total of 8 hub genes (G-protein subunit gamma transducin 2, ubiquitinconjugating enzyme E2 D1, glutamate metabotropic receptor 1, G-protein signaling modulator 1, C-X-C motif chemokine ligand 9, neurotensin, purinergic receptor P2Y1 and ring finger protein 41) were screened from the network. The present study may contribute to the elucidation of the hepatic pathology of T2DM. Copyright: © Chen et al.

Introduction

Type 2 diabetes mellitus (T2DM), resulting from insulin resistance and impaired β-cell function, constitutes a major health problem throughout the world (1). Exploration of the underlying pathological mechanisms and potential therapeutic targets for T2DM is becoming increasingly important (2).

The liver is involved in glucose metabolism, including gluconeogenesis, glycogenolysis, glycogenesis and insulin extraction (3). Dysregulation of glucose metabolism in the liver contributes to the development of T2DM (4). Disruption in the process of hepatic glucose release gives rise to insulin resistance or diabetes and liver diseases may exacerbate insulin resistance by disturbing the physiological effects of insulin on liver cells (5). A previous study reported that targeted inactivation of the hepatic insulin receptor gene resulted in diabetes-like symptoms, demonstrating a direct involvement of insulin regulation in liver metabolism (6). A further study also revealed that selective inactivation of insulin to disrupt hepatic glucose release and fatty acid synthesis led to insulin resistance in the liver, further corroborating that the liver is a significant target for the effect of insulin (7). Impaired fatty acid metabolism in the liver also causes the development of T2DM (8–10). In addition, a clinical study revealed an elevated incidence of newonset diabetes when patients received liver grafts with steatosis, which is strongly linked to hepatic insulin resistance (11).

Genomic data relevant to various diseases are archived in public repositories that are easily accessed to obtain meaningful information and to make novel discoveries (12). Searching in public repositories has been widely applied to investigate the pathology of T2DM, including the identification of underlying pathways and coexpression networks in islets of patients with T2DM (13–15). The gene expression in the liver of a T2DM mouse model has also been analyzed (16). However, to the best of our knowledge, differentially expressed genes (DEGs) in the liver of T2DM patients vs. subjects with normal glucose tolerance (NGT) have remained to be identified. Therefore, the mechanisms underlying the putative hepatic pathology of T2DM remain to be explored.

In the present study, hepatic DEGs in subjects with T2DM vs. NGT were identified, and subsequently, functional enrichment analysis was performed. A protein-protein interaction (PPI) network was also built to identify hub genes. The results of the present study may contribute towards the elucidation of the hepatic pathology of T2DM.

Materials and methods

Microarray data

A gene expression profile (accession no. GSE23343) was obtained from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/). The GEO database stores abundant highthroughput data, particularly those generated by DNA microarray technology (17). A total of 10 patients (6 males and 4 females) with T2DM and 7 subjects (4 males and 3 females) with NGT were included in this GEO dataset, and their clinical characteristics are available from the supplementary information online (18). The array data were acquired from the Affymetrix Human Genome U133 Plus 2.0 array [GPL570; transcript (gene) version].

DEG analysis

The gene expression profiles of liver samples from subjects with T2DM and NGT in the dataset GSE23343 were compared to identify DEGs. This analysis was performed using GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/) through Rbased analysis of the microarray data (12). |log (fold change)|≥1 and P<0.05 were the cut-off criteria. A heatmap of these DEGs was drawn using MeV 4.9.0 (https://sourceforge.net/projects/mevtm4/).

Enrichment analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the DEGs were performed using the Database of Annotation Visualization and Integrated Discovery (DAVID 6.8; http://david.ncifcrf.gov/) (19). The GO categories were biological process (BP), molecular function (MF) and cellular component (CC). P<0.05 was considered to indicate a statistically significant difference. The results of the enrichment analysis were visualized in a bubble chart using the OmicShare tools 3.0, a free online platform for data analysis and visualization (http://www.omicshare.com/tools).

PPI network analysis

The Search Tool for the Retrieval of Interacting Genes (STRING 10.5; http://stringdb.org/) was used to construct a PPI network. This website offers predicted and verified interactions among numerous proteins (20). A combined score >0.7 was selected as the cutoff criterion. Subsequently, the screened PPI network was imported into Cytoscape3.2.1 (http://www.cytoscape.org/) to identify critical gene modules and hub genes. Nodes with a high degree (≥2 fold the median number of connections with other nodes) were considered as significant nodes and nodes with a higher degree (≥5 fold the median number of connections with other nodes) were considered as hub nodes. Submodules of the network were screened using Molecular Complex Detection (MCODE 1.4.2) (21), with the criteria of node number >10 and MCODE score >10. Finally, enrichment analysis of the submodules was performed using DAVID.

Results

A total of 1,320 DEGs in liver samples of patients with T2DM vs. NGT samples were identified, including 698 up- and 622 downregulated genes. The heat-map of the top 50 up- and top 50 downregulated genes is presented in Fig. 1.

Figure 1.

Heat-map of the top 50 up- and top 50 downregulated DEGs in T2DM (P<0.05). The red color represents a higher expression value, whereas the blue color represents a lower expression value. DEG, differentially expressed gene; T2DM, type 2 diabetes mellitus.

GO analysis

In the GO category BP, upregulated genes were mainly enriched in positive regulation of transcription from RNA polymerase II (RNAP II) promoter, cell adhesion, inflammatory response, positive regulation of apoptotic process and extracellular matrix organization (Table I), whereas downregulated genes were mainly associated with signal transduction, multicellular organism development, positive regulation of GTPase activity, visual perception and axon guidance (Table II). In the GO category MF, upregulated genes were mainly involved in calcium ion binding, extracellular matrix structural constituent, SMAD binding, Rho guanylnucleotide exchange factor activity and 3′,5′-cyclic AMP phosphodiesterase activity (Table I), whereas downregulated genes were mainly involved in actin binding, receptor activity, RNAP II transcription factor activity, sequence-specific DNA binding, calmodulin binding and protein tyrosine kinase activity (Table II). Finally, concerning the GO category CC, upregulated genes were mainly involved in the plasma membrane, integral component of plasma membrane, extracellular region, cell junction and cytoskeleton (Table I), whereas downregulated genes were mainly involved in nuclear envelope, myosin complex, microvillus growth cone membrane and actomyosin (Table II).

Table I.

GO analysis of up- and downregulated genes in type 2 diabetes mellitus (P<0.05).

A, Upregulation
Category/term	N (%)	P-value
BP
GO:0045944-Positive regulation of transcription from RNA polymerase II promoter	39 (5.972)	0.049
GO:0007155-Cell adhesion	23 (3.522)	0.017
GO:0006954-Inflammatory response	20 (3.063)	0.017
GO:0043065-Positive regulation of apoptotic process	17 (2.603)	0.016
GO:0030198-Extracellular matrix organization	16 (2.450)	0.001
CC
GO:0005886-Plasma membrane	151 (23.124)	0.001
GO:0005887-Integral component of plasma membrane	64 (9.801)	0.000
GO:0005576-Extracellular region	59 (9.035)	0.049
GO:0030054-Cell junction	24 (3.675)	0.008
GO:0005856-Cytoskeleton	19 (2.910)	0.024
MF
GO:0005509-Calcium ion binding	31 (4.747)	0.033
GO:0005201-Extracellular matrix structural constituent	8 (1.225)	0.003
GO:0046332-SMAD binding	7 (1.072)	0.001
GO:0005089-Rho guanyl-nucleotide exchange factor activity	7 (1.072)	0.025
GO:0004115-3′,5′-cyclic-AMP phosphodiesterase activity	4 (0.613)	0.009
B, Downregulation
Category/term	N (%)	P-value
BP
GO:0007165-Signal transduction	46 (7.931)	0.002
GO:0007275-Multicellular organism development	25 (4.310)	0.003
GO:0043547-Positive regulation of GTPase activity	25 (4.310)	0.007
GO:0007601-Visual perception	12 (2.069)	0.011
GO:0007411-Axon guidance	9 (1.552)	0.043
CC
GO:0005635-Nuclear envelope	9 (1.552)	0.039
GO:0016459-Myosin complex	5 (0.862)	0.032
GO:0005902-Microvillus	5 (0.862)	0.049
GO:0032584-Growth cone membrane	3 (0.517)	0.011
GO:0042641-Actomyosin	3 (0.517)	0.033
MF
GO:0003779-Actin binding	14 (2.414)	0.022
GO:0004872-Receptor activity	12 (2.069)	0.020
GO:0000981-RNA polymerase II transcription factor activity, sequence-specific DNA binding	11 (1.897)	0.010
GO:0005516-Calmodulin binding	11 (1.897)	0.020
GO:0004713-Protein tyrosine kinase activity	10 (1.724)	0.006

BP, biological process; CC, cellular component; MF, molecular function; GO, gene ontology.

Table II.

Kyoto Encyclopedia of Genes and Genomes pathway analysis of up- and downregulated genes in type 2 diabetes mellitus (P<0.05).

A, Upregulated genes
Term	N (%)	P-value	Genes
hsa05202: Transcriptional misregulation in cancer	11 (1.685)	0.039	MAX, CD86, FLI1, SP1, CCND2, PML, ETV1, MDM2, JMJD1C, ETV5, MYCN
hsa04620: Toll-like receptor signaling pathway	10 (1.531)	0.006	IFNA2, CD86, IFNA7, MAPK14, CXCL9, MAPK10, CXCL11, TLR6, TLR8, SPP1
hsa04750: Inflammatory mediator regulation of	9 (1.378)	0.012	PRKCQ, PLA2G4A, IL1R1, PTGER4, MAPK14,
TRP channels			F2RL1, MAPK10, HTR2B, PRKCB
hsa04724: Glutamatergic synapse	9 (1.378)	0.027	SLC17A8, PLA2G4A, GNGT2, GRIK1, GRIN1, SLC38A1, GRM1, SHANK2, PRKCB
hsa04974: Protein digestion and absorption	8 (1.225)	0.021	SLC8A1, COL14A1, COL13A1, PRCP, COL1A2, COL12A1, ATP1A1, COL5A2
B, Downregulated genes
Term	Count	P-value	Genes
hsa04144: Endocytosis	12 (2.069)	0.031	ARFGAP1, IGF1R, CBLC, RET, PIP5KL1, FOLR1, SNX5, RAB35, KIF5C, CYTH4, GRK4, DNM1
hsa04530: Tight junction	8 (1.379)	0.035	SHROOM4, MYH2, EXOC4, MYH14, MYH8, CLDN23, MYH7B, AKT2
hsa05218: Melanoma	6 (1.034)	0.022	FGF6, IGF1R, FGF18, CDKN2A, AKT2, FGF4

Hsa, Homo sapiens; TRP, transient receptor potential; MAX, myc-associated factor X; CD86, T-lymphocyte activation antigen CD86; FLI1, friend leukemia integration 1 transcription factor; SP1, transcription factor sp1; CCND2, cyclin-D2; PML, promyelocytic leukemia; ETV1, ETS translocation variant 1; MDM2, E3 ubiquitin-protein ligase mdm2; JMJD1C, jumonji domain-containing protein 1C; ETV5, ETS variant transcription factor 5; MYCN, N-myc proto-oncogene; IFNA2, interferon α-2; IFNA7, interferon α-7; MAPK14, mitogen-activated protein kinase 14; CXCL9, C-X-C motif chemokine 9; MAPK10, mitogen-activated protein kinase 10; CXCL11, C-X-C motif chemokine 11; TLR6, Toll-like receptor 6; TLR8, Toll-like receptor 8; SPP1, secreted phosphoprotein 1; PRKCQ, protein kinase C θ; PLA2G4A, phospholipase A2 group IVA; IL1R1, interleukin-1 receptor type 1; PTGER4, prostaglandin E receptor 4; F2RL1, coagulation factor II receptor-like 1; HTR2B, 5-hydroxytryptamine receptor 2B; PRKCB, protein kinase C β; SLC17A8, solute carrier family 17 member 8; GNGT2, G protein subunit γ transducin 2; GRIK1, glutamate receptor ionotropic kainate 1; GRIN1, glutamate receptor ionotropic NMDA 1; SLC38A1, solute carrier family 38 member 1; GRM1, glutamate metabotropic receptor 1; SHANK2, SH3 and multiple ankyrin repeat domains 2; SLC8A1, solute carrier family 8 member A1; COL14A1, collagen type XIV α 1 chain; COL13A1, collagen type XIII α 1 chain; PRCP, prolylcarboxypeptidase; COL1A2, collagen type I α 2 chain; COL12A1, collagen type XII α 1 chain; ATP1A1, ATPase Na+/K+ transporting subunit α 1; COL5A2, collagen type V α 2 chain; ARFGAP1, ADP-ribosylation factor GTPase-activating protein 1; IGF1R, insulin-like growth factor 1 receptor; CBLC, E3 ubiquitin-protein ligase CBL-C; RET, ret proto-oncogene; PIP5KL1, phosphatidylinositol-4-phosphate 5-kinase like 1; FOLR1, folate receptor α; SNX5, Sorting nexin-5; RAB35, Ras-related protein Rab-35; KIF5C, kinesin family member 5C; CYTH4, Cytohesin-4; GRK4, G protein-coupled receptor kinase 4; DNM1, dynamin 1; SHROOM4, shroom family member 4; MYH2, myosin heavy chain 2; EXOC4, exocyst complex component 4; MYH14, myosin heavy chain 14; MYH8, myosin heavy chain 8; CLDN23, Claudin-23; MYH7B, myosin-7B; AKT2, RAC-β serine/threonine-protein kinase; FGF6, fibroblast growth factor 6; FGF18, fibroblast growth factor 18; CDKN2A, cyclin dependent kinase inhibitor 2A; FGF4 fibroblast growth factor 4.

KEGG pathway analysis

Upregulated genes were mainly enriched in transcriptional misregulation in cancer, Toll-like receptor (TLR) signaling pathway, inflammatory mediator regulation of transient receptor potential (TRP) channels, glutamatergic synapse, and protein digestion and absorption, whereas downregulated genes were mainly associated with endocytosis, tight junction and melanoma (Table II). The results of the enrichment analysis were visualized in Figs. 2 and 3, respectively.

Figure 2.

GO and KEGG pathway analysis of the top 20 upregulated genes in type 2 diabetes mellitus (P<0.05). The Y-axis represents GO categories, including BP, MF and CC, whereas the Xaxis represents the enrichment factor. The enrichment factor is the ratio of the number of DEGs annotated to the term to the number of all genes annotated to it. In addition, the dot size represents the number of DEGs annotated to the term, whereas the dot color indicates the significance of gene enrichment. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEG, differentially expressed gene; TRP, transient receptor potential; BP, biological process; CC, cellular component; MF, molecular function.

Figure 3.

GO and KEGG pathway analysis of the top 18 downregulated genes in type 2 diabetes mellitus (P<0.05). The Y-axis represents GO categories, including BP, MF and CC, whereas the X-axis represents the enrichment factor. The enrichment factor is the ratio of the number of DEGs annotated to the term to all the genes annotated to it. In addition, the dot size represents the number of DEGs annotated to the term, whereas the dot color indicates the significance of gene enrichment. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEG, differentially expressed gene; BP, biological process; CC, cellular component; MF, molecular function.

As presented in Fig. 4, the PPI network of DEGs consisted of 443 nodes and 996 edges. A total of 11 genes were selected as candidates for hub genes. In addition, two submodules were selected, one of which had 28 nodes and 197 edges, while the other module had 14 nodes and 91 edges (Fig. 5). Finally, as presented in Table III, eight hub genes involved in these two submodules were identified, including Gprotein subunit gamma transducin 2 (GNGT2), ubiquitin-conjugating enzyme E2 D1 (UBE2D1), glutamate metabotropic receptor 1 (GRM1), G-protein signaling modulator 1 (GPSM1), C-X-C motif chemokine ligand 9 (CXCL9), neurotensin (NTS), purinergic receptor P2Y1 (P2RY1) and ring finger protein 41 (RNF41). However, no enrichment was observed in these two submodules.

Figure 4.

Protein-protein interaction network of the DEGs in type 2 diabetes mellitus. The red nodes represent the upregulated genes, whereas green nodes represent the downregulated genes. Specifically, light blue nodes indicate the significant genes with a high degree (≥2-fold the median number of connections with other nodes), whereas dark blue nodes indicate the hub genes with a higher degree (≥5-fold the median number of connections with other nodes). In addition, purple nodes indicate other genes associated with the DEGs that were identified. The node size indicates the node counts. DEG, differentially expressed gene.

Figure 5.

(A) One of the two significant submodules. (B) Core protein-protein interaction network of differentially expressed genes in type 2 diabetes mellitus. (C) The other significant submodule.

Table III.

Hub nodes in the network of differentially expressed genes in type 2 diabetes mellitus.

Hub node	Description	Degree	MCODE score	Count	Up/downregulation
GNGT2	G protein subunit gamma transducin 2	39	14	41	Up
UBE2D1	Ubiquitin-conjugating enzyme E2 D1	19	13	29	Up
GRM1	Glutamate metabotropic receptor 1	17	13	20	Up
GPSM1	G-protein signaling modulator 1	16	14	18	Down
CXCL9	C-X-C motif chemokine ligand 9	16	14	24	Up
NTS	Neurotensin	16	13	38	Up
P2RY1	Purinergic receptor P2Y1	16	13	17	Up
RNF41	Ring finger protein 41	16	13	17	Up

Discussion

In the present study, 698 up- and 622 downregulated DEGs were screened from the hepatic genes of patients with T2DM and normal subjects. GO term analysis revealed that the upregulated DEGs were mainly associated with positive regulation of transcription from RNAP II promoter, cell adhesion, inflammatory response, positive regulation of apoptotic process and extracellular matrix organization. Hepatocyte nuclear factor 4 (HNF4) regulates numerous pivotal metabolic pathways and exert significant effects on recruiting RNAP II to synthesize gene promoters. Abnormalities in the hepatic HNF4 transcription network are accountable for diabetes and fatty liver (22). Alterations in cell adhesion may disturb significant cellular processes, leading to the causation of various diseases. Targeted inactivation of carcinoembryonic antigen-related cell adhesion molecule 1 (CEACAM1) in the liver was reported to cause insulin resistance and promote hepatic adipogenesis, suggesting a critical role of CEACAM1 in regulating insulin clearance in the liver (23). Hyperglycemiainduced oxidative stress induces liver tissue injury and the ensuing derangement of protein, carbohydrate and lipid metabolism leads to increased oxidative stress, further triggering the inflammatory cascade (24). Hepatocyte inflammation significantly downregulates insulin signaling components, including insulin receptor substrate (IRS)-1, IRS-2, PI3K, Akt and mTOR (25). Inflammatory regulators induced by hepatocyte apoptosis-associated damage are able to regulate the insulin signaling pathway, and these insulin resistanceassociated regulators may, in turn, affect hepatocyte apoptosis (5). Endoplasmic reticulum stressinduced apoptosis of hepatocytes and adipocytes is also important in the development of diabetes, characterized by increased insulin resistance (26). The downregulated DEGs were mainly involved in signal transduction, multicellular organism development and positive regulation of GTPase activity. In the process of metabolic alterations, cellular responses to extracellular stimulation require signal transduction, contributing to physiological events including increased uptake of blood glucose (27). GTPases are also important in signal transduction at the intracellular domain of transmembrane receptors (28).

The KEGG pathway enrichment analysis indicated that the upregulated DEGs were accumulated in the TLR signaling pathway, inflammatory mediator regulation of TRP channels and protein digestion and absorption, and that the downregulated DEGs were enriched in endocytosis and tight junction. Diabetes frequently occurs in combination with other metabolic diseases, including hyperlipidemia, hypertension and non-alcoholic fatty liver disease (29). Deposition of fatty acids in the liver, particularly saturated fatty acids, activates the TLR pathway, which is associated with the inflammatory response (30). Hepatic inflammation is closely correlated with insulin resistance (25). A previous study reported that TRP cation channel subfamily V member 4 effectively regulates the expression of various pro-inflammatory genes in adipose tissue, and that these pro-inflammatory genes are closely associated with insulin resistance (31). Tight-junction proteins, besides their function as integral proteins of tight junctions that form barriers in the gut and the liver, may also be expressed outside the tight junction to regulate signaling, trafficking and gene expression. A hallmark is their regulation of epithelial-to-mesenchymal transition (32). A previous study demonstrated that the endocytosis impairment of specific ligands or other macromolecules may represent an important pathology mechanism in diabetes (33). The biological processes and pathways identified and discussed above may indicate an important role of the liver in the pathology of T2DM.

In the present study, the following eight hub genes were also selected: GNGT2, UBE2D1, GRM1, GPSM1, CXCL9, NTS, P2RY1 and RNF41. GNGT2 was reported to be involved in β-arrestin-1induced Akt phosphorylation and NF-κB activation (34). Activation of NFκB in the liver may result in hepatic insulin resistance (5). The low-density lipoprotein (LDL) receptor is indispensable for the uptake of LDL cholesterol and for regulating the levels of plasma lipoprotein (35). The E3 ubiquitin ligase inducible degrader of LDL receptor/ubiquitin-conjugating enzyme E2D complex is effectively responsible for determining LDL receptor activity (36). Sirtuin 1, a type of nicotinamide adenine dinucleotide-dependent deacetylase, also regulates the pathogenesis of metabolic disease, aging and tumorigenesis (37). Sirtuin 1-mediated epigenetic regulation of the expression of the metabotropic glutamate receptor 1/5 (encoded by the GRM1/5 gene) was reported to be involved in the development of neuropathic pain in a rat model of T2DM (38,39). The GPSM1 locus has been demonstrated to be associated with the insulinogenic index and with the fasting glucose level (40,41), and GPSM1 has been identified as one of the three novel T2DM loci in East Asian populations (42). Previous studies have also suggested an important role of CXCL9 in diabetic neuropathy, diabetic retinopathy and diabetic nephropathy (43–45). Advanced glycation end products were reported to promote apoptosis and inflammation in mouse podocytes via CXCL9-regulated activation of the JAK2/STAT3 pathway (46). The fasting plasma levels of pro-NTS produced in equimolar amounts with NTS were indicated to be positively associated with the risk of diabetes, cardiovascular disease and mortality (47). Obese and insulin-resistant patients had higher plasma concentrations of proNTS, and NTS-deficient mice on a highfat diet exhibited lower levels of fasting plasma glucose and insulin compared with wildtype mice (48). Furthermore, P2Y1 receptorknockout mice exhibited high levels of plasma insulin, plasma glucose and increased body weight, indicating an important regulatory role of the P2Y1 receptors in glucose homeostasis (49). Finally, RNF41, an E3 ubiquitin ligase, was identified to be essential for activation of the type 1 interferon pathway to sustain insulin sensitivity in the muscle tissue of obese patients (50). Pancreatic islet β-cells require normal mitochondrial function in terms of their high metabolic activity. Stabilization of the C-type lectin domain-containing 16A/RFP41/ubiquitin-specific peptidase 8 mitochondrial autophagy complex is essential for cellular respiration and insulin secretion. However, a study reported that elevated levels of glucose and fatty acids destabilized the complex, causing β-cell apoptosis (51).

In conclusion, a comprehensive analysis of hepatic DEGs in T2DM was performed in the present study, revealing an important role of the liver in the pathological mechanisms of T2DM. However, considering the absence of clinical validation in the present study, further investigation of these mechanisms underlying the hepatic pathology of T2DM is required to confirm these results.