Downregulation of ALDH6A1 is a New Marker of Muscle Insulin Resistance in Type 2 Diabetes Mellitus.

PURPOSE: Skeletal muscle insulin resistance (IR) is an important etiology of type 2 diabetes mellitus (T2DM); however, its molecular mechanism is yet to be fully defined. This study attempted to identify the gene expression patterns and molecular disorders in T2DM patients' skeletal muscle samples.
METHODS: First, the difference in genetic expression among GSE25462 data was analyzed. Next, PPI network analysis of differential genes was carried out, after which the maladjustment module was identified. Then, an enrichment analysis and gene set enrichment analysis (GSEA) were carried out. Finally, the transcription factors that regulate the modular genes by raid were predicted.
RESULTS: Most differentially expressed genes were found to be able to form an interaction network and cluster into 9 modules. These modular genes were shown to possess a significant correlation with immune inflammation and metabolic response. Importantly, the top 15 genes of area under receiver operating characteristic curve (AUC) were identified, and the expression of 10 genes by GSE12643, GSE18732 and GSE29221 was confirmed. The expression and AUC value of ALDH6A1 were then verified according to three sets of data, where ALDH6A1 was found to be negatively correlated with follicular helper T cells. However, among the predicted transcription regulators, HDAC was shown to have a better regulatory effect.
CONCLUSION: The findings highlight that the dysregulation of ALDH6A1 expression in IR of T2DM may serve as a potential therapeutic target. ALDH6A1 is involved in the immune inflammation and metabolic pathways.

Introduction

Diabetes has been recognized as one of the four major non-infectious diseases, which has resulted in about 1.6 million deaths worldwide. It is regarded as the third risk factor of premature death due to complications, including cardiovascular diseases.1 It has been estimated that the global prevalence of diabetes is 8.8% (as of 2017), which is expected to rise to 9.9% by 2045.2 In China, the prevalence of diabetes reached 11.6% in 2010, affecting about 113.9 million adults.3 Type 2 diabetes mellitus (T2DM) is the most common type of diabetes, accounting for 95% of all global cases.4 Skeletal muscle insulin resistance (IR) has been noted to be an early feature in the progression of T2DM.5

Fasting blood glucose, 2 hours postprandial blood glucose, glycosylated hemoglobin, insulin resistance and other clinical indicators are widely used in the diagnosis and classification of T2DM.6 Compared with the general population, people with type 2 diabetes have an increased risk of cardiovascular disease, heart failure and death.7 Moreover, studies have shown that the excess risk associated with diabetes is mainly caused by hyperglycemia.8,9 The pathogenesis of T2DM includes two major abnormalities: insulin resistance and insulin secretion dysfunction, which lead to an inability in regulating blood glucose levels.10,11 Previous animal experiments have found that IR is related to tyrosine metabolism and affects the development of T2DM.12 Interestingly, studies have also shown that low-grade inflammation is related to T2DM risk, subclinical inflammation can lead to insulin resistance, and is related to features of metabolic syndrome, including hyperglycemia.13,14 Increased tumor necrosis factor (TNF)-α secretion is associated with obesity-related insulin resistance. Obesity is a risk factor for the development of T2DM,15,16 and blood pressure has also been shown to be independently associated with the development of T2DM.17

Recently, researchers have conducted a comprehensive and systematic study on the changes of gene expression in T2DM patients using genomics technology, in which its biomarkers and possible metabolic pathways were studied, providing a theoretical basis for the pathogenesis of T2DM. In order to fully understand the pathological mechanism of T2DM and provide ideas and references for future research on the pathogenesis of T2DM, the genetic expression of skeletal muscle samples is analyzed in the present study. Furthermore, this study attempts to provide a theoretical basis for the disease mechanism and molecular markers related to IR in the skeletal muscle of T2DM patients.

Materials and Methods

Identification of Differentially Expressed Genes (DEGs) from Microarray Data

Transcriptome sequencing data of skeletal muscle samples in the GSE25462, GSE12643, GSE18732 and GSE29221 datasets were collected from the GEO (gene expression omnibus) database.18 GSE25462 included 10 patients with T2DM and 15 healthy controls; GSE12643 included 10 obese T2DM patients and 10 healthy controls; GSE18732 included 45 patients with T2DM and 47 healthy controls; and GSE29221 included 3 patients with T2DM and 3 healthy controls. The pre-processing of the original data of each microarray included background correction, normalization and log 2 conversion. The R language software package, DESeq, was used to determine the DEGs between the T2DM samples and control samples. P < 0.05 was taken as the cut-off criterion.

Construction of the Protein–Protein Interactions (PPI) Network

Using the online search tool to search the interaction gene (String) website (), the comprehensive score was >0.4, and the PPI network19 of DEGs was then constructed. Then, the molecular complexity detection (MCODE) plug-in20 in Cytoscape was used to filter the modules in the PPI network according to the default settings.

Enrichment Analysis and GSEA of Module Genes

The R-package clusterProfiler21 was used to enrich and analyze the module genes in the gene ontology (GO) and Kyoto genome Encyclopedia (KEGG). P < 0.05 was defined as the cut-off for the significant function and pathway analysis. Additionally, the mgeneSim function of the R package in GOSemSim was used to calculate the functional similarity of GO.

In order to identify significantly dysregulated biological pathways in T2DM, GSEA22 was performed by the clusterProfiler R package under the functional annotations of the KEGG database. According to the default weighted enrichment statistical method, each analysis was repeated 1000 times. Activated and suppressed pathways with P < 0.05 in each dataset were then merged, after which those with a higher frequency (found in ≥3 datasets) were identified as dramatically changed KEGG pathways in T2DM.

Prediction of Transcription Factors

The transcription factor (TF) with a regulatory effect on the module gene was then predicted. The human transcription factor target data in TRRUST version 2 was downloaded so as to select the TF related to the module gene. Then, the rentrez R package was utilized to acquire the interactions in the NCBI database and filter the TF. Finally, the TFs that regulated more than 20 module genes were screened, which were tested using the hypergeometric method.

Results and Discussion

Identification of Differentially Expressed Genes in the Skeletal Muscle of Diabetic Patients

In order to identify the DEGs related to IR in skeletal muscle of T2DM, the gene expression of GSE25462 samples in the GEO database was analyzed. Accordingly, 1772 DEGs () were identified, of which 1080 were upregulated and 692 were downregulated (Figure 1A and B). The top 10 upregulated and downregulated DEGs are listed in Table 1. In order to better understand the action network of these DEGs, PPI network identification was carried out. Using the MCODE plug-in, 732 genes with high interaction were identified, which were then clustered into 9 modules (Figure 1C). These nine modules may represent different molecular mechanisms related to IR of T2DM in skeletal muscle.

Figure 1

PPI network formed by differentially expressed genes in the skeletal muscle of diabetic patients. (A) Volcanic map of differentially expressed genes in the skeletal muscle of diabetic patients. (B) Heatmap of differentially expressed genes in the skeletal muscle of diabetic patients. Red node represents an upregulated gene, while green represents a downregulated gene. (C) Nine subnetworks of PPI network for differentially expressed genes.

Table 1

Top 10 Upregulated and Downregulated Differentially Expressed Genes

DEG	logFC	AveExpr	P value
CALML6	1.331931	7.067413	0.003534
PDK4	1.327601	8.387643	0.012429
ABRA	1.277571	9.459666	2.95E-06
PVALB	1.217981	5.047922	0.006289
IL32	0.936305	8.338253	0.015874
MYH1	0.852828	12.69676	0.041078
OTUD1	0.741183	6.076541	0.004083
APOD	0.70477	9.926132	0.004646
LMOD1	0.663079	7.294198	0.002226
CYR61	0.629836	5.686033	0.000434
ALDH6A1	−0.6227	6.718785	0.000105
CIART	−0.63551	8.051595	0.003097
PTGR2	−0.65085	6.332798	0.005591
LRP1B	−0.71353	5.595935	0.0357
LOC100134445	−0.71682	6.749796	0.041365
GCLM	−0.72606	7.03357	0.001565
ENPP5	−0.79438	6.245947	0.003334
TECRL	−0.80555	8.394449	0.004057
RCAN2	−0.86372	9.346666	0.000415
LDHB	−0.87802	10.1337	0.00303

Transcriptional Regulators of Regulatory Module Genes

Transcriptional regulation plays an important role in gene function. Therefore, the transcription factors that regulate the module genes were predicted. Among the 333 transcription factors predicted and screened through the interactions search in the NCBI gene database, 16 transcription factors that interacted with more than 20 target genes were identified. Hypergeometric analysis of the 16 TFs showed a high correlation between transcription factors, especially HDAC, which was itself coupled with transcription factors in order to play a regulatory role (Figure 2A). Numerous types of target regulatory relationships exist between key regulatory factors and target genes (Figure 2B). Therefore, it is believed that HDAC is the core gene in the regulatory network.

Figure 2

Transcription factors that regulate modular genes. (A) Sixteen transcription factors with more than 20 regulatory module genes were identified as key transcription factors. (B) Targeting relationship between regulatory factors and target genes. The red node is a transcription factor, while green is a gene regulated by more than two regulators and yellow is a gene regulated by only one regulator.

Biological Function and Signal Pathway of Modular Genes

In order to study the biological function of the identified module genes, the data was analyzed by an enrichment analysis and GSEA. The enrichment analysis showed that the module genes participated in the 4577 BP, 448 CC, 779 MF and 229 KEGG signaling pathways (). GO functions mainly included the regulation of leukocyte activation, positive regulation of neurogenesis, regulation of cytoskeleton organization and other biological functions (Figure 3A). Among them, module 1 was mainly found to be related to the regulation of the muscle system process. Module 2 was shown to be mainly related to cell morphogenesis involved in neuron differentiation. Meanwhile, module 3 was mainly related to the regulation of leukocyte activation; module 4 was mainly related to the regulation of body fluid levels; module 5 was mainly related to the small molecule catabolic process; module 6 was mainly related to cell amino acid metallic process; module 7 was mainly related to muscle organ development; module 8 was mainly related to hormone secretion; and module 9 was mainly related to the steroid metabolic process.

Figure 3

The biological function and signal pathway of module genes. (A) Module genes participate in biological processes. (B) The KEGG pathway of modular genes. (C) Fifteen GSEA pathways were the same as the KEGG pathway.

However, the KEGG pathways involved in the module genes were observed to be mainly neuroactive live receiver interaction, MAPK signaling pathway and chemokine signaling pathway (Figure 3B). GSEA showed that 15 pathways (Figure 3C) had the same signal pathway as that of KEGG enrichment, including the IL-17 signaling pathway, osteoclast differentiation, and estrogen signaling pathway. These results indicated that the differentially expressed genes in skeletal muscle of diabetic patients are mainly related to immune and neural responses.

Identification of the Key Dysregulation Gene

In order to identify the key genes that play an important role in skeletal muscle IR in diabetic patients, the area under receiver operating characteristic curve (AUC) of the module gene was calculated, and genes with the top 15 AUC values were then screened (Figure 4A). Importantly, the AUC values of 15 genes were verified by GSE12643, GSE18732 and GSE29221, of which 10 genes were verified (Figure 4B), which were upregulated or downregulated in GSE25462 (Figure 4C). Therefore, these 10 genes were believed to be potential key genes for insulin resistance in the skeletal muscle of diabetic patients. It was found that the functional similarity of SRF, JUND and NR4A1 was high (Figure 4D). Among them, the expression of the ALDH6A1 gene and AUC was verified by three sets of data (Figure 4E). Moreover, the enrichment analysis showed that ALDH6A1 was significantly related to metabolism-related functions and signaling pathways (Figure 4F). In addition, the correlation between ALDH6A1 and immune cells was calculated (Figure 4G). ALDH6A1 was found to have a negative correlation with most immune cells in the diabetic samples. In particular, the negative correlation with follicular helper T (TFH) cells was strongest (Figure 4H). Therefore, ALDH6A1 can serve as a key candidate gene for insulin resistance in diabetic skeletal muscle.

Figure 4

Key dysfunctional genes in skeletal muscle of type 2 diabetes mellitus. (A) The module genes with top 15 AUC value. (B) GSE12643, GSE18732 and GSE29221 dataset verify the AUC values of 10 genes. (C) Expression of 10 genes in GSE25462. (D) The functional similarity of GO function among the 10 genes. (E) Three sets of data were used to verify the expression and AUC of ALDH6A1. (F) ALDH6A1 is involved in the signaling pathway, as well as other genes contained in the signaling pathway. (G) Correlation between ALDH6A1 and 24 immune cells. (H) ALDH6A1 has the strongest negative correlation with TFH.

It is both clinically and scientifically important to study the main molecular defects present in the pathophysiology of human T2DM. This study showed that 1772 DEGs were present in the skeletal muscle of T2DM patients compared with those of healthy controls. Gene expression analysis can be used to identify markers or early defects related to disease phenotypes. Previous studies have shown that the expression of mitochondrial oxidative phosphorylation related nuclear coding genes in T2DM patients was reduced.23 In addition, in T2DM patients, insulin signaling and glucose uptake related genes and signaling pathways were also shown to be changed significantly.24 These DEGs can be identified by a PPI network with significant interaction. Accordingly, 9 subnetworks of gene interaction were obtained in this study, each of which may represent a mechanism of action.25

The conducted enrichment analysis demonstrated that the module genes were mainly involved in the regulation of immune and metabolic processes. First, genes related to immune inflammation were shown to play an important role in the progression of T2DM.10 Due to the low degree of inflammation in T2DM, these genes can affect insulin sensitivity and insulin secretion.26 Various studies have shown that low-level inflammatory markers may be helpful in distinguishing classic adult type 1 diabetes and T2DM.27 The present data showed that most modular genes were involved in the MAPK and IL-17 signaling pathways. P38 mitogen activated protein kinase (MAPK) signaling pathway was also found to be involved in the apoptosis of a variety of inflammatory cells, which may serve as a potential therapeutic target of T2DM characterized by low-grade inflammation.28 In addition, p38 γ MAPK activation has been shown to inhibit insulin signal transduction at the level of the insulin receptor substrate.29 However, IL-17 has been shown to play an important role in inflammation, insulin resistance and type 2 diabetes.30 Compared with the healthy controls, the expression of IL-17 in Foxp3-CD4 T cells and Foxp3 + Treg cells in T2DM patients was demonstrated to be significantly higher.31 Moreover, the current analysis also found that metabolic function was found to be significantly enriched in the T2DM skeletal muscle samples. Past studies have analyzed the metabolites of 1940 patients with T2DM, in which a rise in isoleucine, leucine and proline led to an increase in the relative risk of T2DM.32 In addition, glycine was observed to be positively correlated with insulin sensitivity, and β - hydroxypyruvate, an intermediate metabolite of glycine, was also considered to be a predictor of diabetes.33,34 The metabolism of bile acids, short chain fatty acids and metabolites cooperate in the metabolic regulation system, leading to IR.35 Therefore, it is believed that the imbalance of skeletal muscle gene expression in T2DM patients may serve as a potential factor leading to IR.

Importantly, ALDH6A1 was shown to play a vital role in the diagnosis and treatment of IR in the skeletal muscle of T2DM patients. ALDH6A1 has been identified as a new adipose tissue marker related to T2DM obese people.36 In the IR atherosclerotic mouse model, ALDH6A1 is the highest ranked genes for the oxidative phosphorylation pathway.37 ALDH6A1 encoding methylmalonate-semialdehyde dehydrogenase, consistent with impaired branched-chain amino acids (BCAAs) metabolism.38 BCAAs functional attenuation has previously been shown to be involved in T2DM or IR associated mitochondrial dysfunction and impairment of oxidative metabolism.39 In this study, the mRNA level of ALDH6A1 in T2DM skeletal muscle was found to be relatively downregulated. However, whether the decrease in ALDH6A1 transcription activity in skeletal muscle is the cause or result of IR remains unanswered.

Histone deacetylase (HDAC) was found to be the most interactive gene among transcription factors by predicting the transcription regulators of the modular genes. Histone modifying functions of HDACs can regulate gene expression and affect pancreatic development.40 HDAC is the therapeutic target of T2DM, and its inhibitor has been demonstrated to be a new therapeutic drug of T2DM.41,42 HDAC inhibitors, which constitute a new class of drugs to treat diabetes and its complications, play an important role in preventing β - cell resistance to insulin, destroying β - cells and preventing cytokine mediated attacks on pancreatic cells.43

Conclusions

This study’s findings demonstrate that abnormal immunoinflammation and metabolic response may play a role in the pathogenesis of skeletal muscle insulin resistance in patients with T2DM. In addition, ALDH6A1 may serve as a potential marker, and HDAC is a transcriptional regulator of T2DM maladjusted genes rather than only a class of enzyme.