Comprehensive analysis of endoplasmic reticulum stress-related mechanisms in type 2 diabetes mellitus

Abstract

BACKGROUND

The endoplasmic reticulum (ER) is closely related to a wide range of cellular functions and is a key component to maintain and restore metabolic health. Type 2 diabetes mellitus (T2DM) is a serious threat to human health, but the ER stress (ERS)-related mechanisms in T2DM have not been fully elucidated.

AIM

To identify potential ERS-related mechanisms and crucial biomarkers in T2DM.

METHODS

We conducted gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA) in myoblast and myotube form GSE166502, and obtained the differentially expressed genes (DEGs). After intersecting with ERS-related genes, we obtained ERS-related DEGs. Finally, functional analyses, immune infiltration, and several networks were established.

RESULTS

Through GSEA and GSVA, we identified several metabolic and immune-related pathways. We obtained 227 ERS-related DEGs and constructed several important networks that help to understand the mechanisms and treatment of T2DM. Finally, memory CD4⁺ T cells accounted for the largest proportion of immune cells.

CONCLUSION

This study revealed ERS-related mechanisms in T2DM, which might contribute to new ideas and insights into the mechanisms and treatment of T2DM.

Key Words: Endoplasmic reticulum stress; Type 2 diabetes mellitus; Biomarkers; Memory CD4⁺ T cells

Core Tip: This study revealed endoplasmic reticulum stress-related mechanisms in type 2 diabetes mellitus (T2DM), which might contribute to new ideas and insights for the mechanisms and treatment of T2DM.

INTRODUCTION

Diabetes is a chronic disease that occurs either when the pancreas does not produce enough insulin or when the body cannot effectively use the insulin it produces[1,2]. Hyperglycemia is a common effect of uncontrolled diabetes and over time leads to serious damage to many of the body’s systems, especially the heart, blood vessels, eyes, kidneys, and nerves[3,4]. Recently, the estimated prevalence of diabetes among children, adolescents, and adults has increased[5,6]. The majority of people with diabetes have type 2 diabetes mellitus (T2DM)[7]. Simple lifestyle measures have been shown to be effective in preventing or delaying the onset of T2DM[8]. Recently, with the in-depth understanding of the mechanisms of T2DM, many new drugs, such as sodium-glucose cotransporter-2 inhibitors, glucagon-like peptide-1 analogs, and dipeptidyl peptidase-4 inhibitors, have been gradually applied to clinical practice and achieved good results[9-11]. However, the residual risk of these populations remains high, especially when combined with other diseases[12].

The endoplasmic reticulum (ER) is closely related to a wide range of cellular functions and is a key component to maintain and restore metabolic health[13]. Protein handling, modification, and folding in the ER are tightly regulated processes that determine cell function, fate, and survival[14]. Many genetic and environmental damages hinder the ability of cells to correctly fold and post-translationally modify secreted and transmembrane proteins in the ER, resulting in the accumulation of misfolded proteins in this organelle, which is called ER stress (ERS)[15]. Chronic ERS is becoming a key factor in more human diseases, including T2DM[16,17]. Recently, the biological mechanisms of ERS in T2DM have been gradually explored. YIPF5 mutations can disrupt the ER-to-Golgi trafficking, thereby resulting in T2DM[16]. Inositol-requiring enzyme 1alpha upregulates miR-200a degradation and stimulates TXINP/NLRP3-pathway-mediated pyroptosis and renal damage in T2DM[18]. Mfn2 plays an important role in ERS, and Mfn2 silencing prevents mitochondrial Ca2⁺ overload-mediated mitochondrial dysfunction[19]. ATF5 is a regulator of ERS and β-cell apoptosis in different models of diabetes mellitus[20]. Lactogens modulate the ERS pathway, causing enhanced β-cell survival and reduced T2DM incidence[21]. The development of ERS for the treatment of T2DM has also emerged in clinical trials. A randomized placebo-controlled crossover trial indicated that decreased ERS may lead to improvement of insulin sensitivity mediated by hyperbaric oxygen[22]. Nevertheless, the role of ERS in T2DM, especially the related markers and mechanisms, is still lacking.

Here, we conducted gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA) in both proliferating myoblasts and differentiated myotubes, which are important in T2DM. Then, the differentially expressed genes (DEGs) and ERS-related DEGs between T2DM patients and healthy populations were investigated, sequentially. Furthermore, functional enrichment analysis [Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG], immune infiltration analysis, and three networks [transcription factor (TF)–mRNA, miRNA–mRNA, and drug–mRNA] were detected to explore the mechanisms and potential therapeutic agents of ERS in T2DM. The flow chart is shown in Figure 1.

Figure 1 Study flow chart. GSEA: Gene set enrichment analysis; GSVA: Gene set variation analysis; ER: Endoplasmic reticulum; DEG: Differentially expressed genes; PPI: Protein–protein interaction; TF: Transcription factor.

MATERIALS AND METHODS

Acquisition and processing of raw data

The raw data of the microarray expression dataset GSE166502[23] and its annotation file GPL10558 (Illumina HumanHT-12 V4.0 Expression BeadChip) were obtained from Gene Expression Omnibus[24]. GSE166502 holds the mRNA expression in proliferating myoblasts and differentiated myotubes in patients with T2DM (n = 13) or controls (n = 13).

GSEA and GSVA

We selected and downloaded c2.cp.v7.2.symbols.gmt gene set data through the GSEA database[25], and conducted GSEA on the proliferating myoblasts and differentiated myotubes through the clusterProfiler package (version 3.14.3)[26]. The statistical process of GSEA was to calculate the enrichment score, estimate the significance, and correct the multiple hypothesis tests. We also selected the same data from GSEA and conducted GSVA. The different pathways were obtained through the limma package (version 3.42.2)[27].

Identification of DEGs

After the processing of raw data, we analyzed the data using the limma package with a fold change and P for DEGs. The threshold of DEGs was |log₂fold change| > 0.263 and P < 0.05 as described previously, and the results were visualized as a heat map and volcano map using the pheatmap package (version 1.0.12).

Acquisition of ERS-related DEGs

GeneCards provides annotated and predicted human gene information, which integrates gene data from about 150 network sources, including genomics, transcriptomics, proteomics, genetics, and clinical and functional information[28]. In this study, ERS-related genes were downloaded through GeneCards with “endoplasmic reticulum stress” as the search keyword. Taking the intersection of DEGs and ERS-elated genes, we got the ERS-related DEGs and the Venn diagram was drawn through the Venndiagram package (version 1.6.20).

Functional enrichment analysis

GO and KEGG pathway analysis can contribute to the interpretation of system-level data and enable discoveries[29,30]. In this work, GO terms and KEGG analysis of ERS-related DEGs and potential molecular complex were carried out using the clusterProfiler package with P < 0.05, and then visualized by the ggplot2 package (version 3.3.3), as described previously[31].

Protein–protein interaction analysis

Protein–protein interaction (PPI) is one of the cores of cellular processing. The analysis of PPI makes the relationships among proteins clear and helps the function explanation of potential protein complexes or functional modules. In this work, PPI information was surveyed using the String database (version 11.0)[32]. The PPI network of ERS-related DEGs was uploaded to Cytoscape (version 3.8.2)[33] and the NetworkAnalyzer plugin was used to further processing and analysis. The cytoHubba plugin was used to select the top 20 key genes[34].

Network analysis

Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining (TRRUST, version 2) manually curated database of human and mouse transcriptional regulatory networks[35]. Current TRRUST contains 8444 and 6552 TF-target regulatory relationships of 800 human and 828 mouse TFs. TRRUST database also provides information on the mode of regulation (activation or repression). miRWalk (version 3.0) stores predicted data obtained with a machine-learning algorithm including experimentally verified miRNA–target interactions[36]. The drug–gene interaction database (DGIdb, version 4.2.0) builds drug–gene interactions mined from DrugBank, PharmGKB, Chembl, Drug Target Commons, Therapeutic Target Database, and others[37]. DGIdb contains > 40000 genes and > 10000 drugs involved in > 100000 drug–gene interactions or belonging to one of 42 potentially druggable gene categories. We obtained the TFs, miRNAs, and drugs of ERS-related DGEs, respectively, and then constructed the regulation relationship networks through Cytoscape.

Correlation analysis of immune infiltration

CIBERSORT (version 1.03) calculates the proportion of different types of cells according to LM22[38]. The proportion of different cell types can be calculated after the nonnegative matrix decomposition of the expression matrix. In this study, the immune infiltration of GSE166502 was analyzed by CIBERSORT, and the infiltration of 22 kinds of immune cells in the sample was analyzed. Finally, we analyzed the correlation between the expression of the top 20 key genes in the PPI network and the immune infiltration.

RESULTS

GSEA and GSVA

Through GSEA, we found that neuroactive ligand–receptor interaction, hypertrophic cardiomyopathy, DNA replication, cell cycle, and cardiac muscle contraction were the top five pathways in proliferating myoblasts (Figure 2A–F). DNA replication, cell cycle, cardiac muscle contraction, neuroactive ligand–receptor interaction, and hypertrophic cardiomyopathy were activated, whereas glycosaminoglycan biosynthesis heparan sulfate, glycosaminoglycan biosynthesis chondroitin sulfate, glycosaminoglycan degradation, other glycan degradation, and lysosome were suppressed (Figure 2G). Other pathways, such as arachidonic acid metabolism, mismatch repair, P53 signaling pathway, metabolism of xenobiotics by cytochrome P450, and prion diseases, were also enriched (Figure 2H). Similarly, we found that viral myocarditis, steroid hormone biosynthesis, hematopoietic cell lineage, focal adhesion, and extracellular matrix (ECM) receptor interaction were the top five pathways in differentiated myotubes (Figure 2I–N). Neuroactive ligand–receptor interaction, gap junction, pathways in cancer, focal adhesion, and ECM receptor interaction were activated, whereas steroid hormone biosynthesis, cardiac muscle contraction, viral myocarditis, hematopoietic cell lineage, and steroid biosynthesis were suppressed (Figure 2O). Vascular endothelial growth factor (VEGF) signaling pathway, cell adhesion molecules cams, mitogen-activated protein kinase (MAPK) signaling pathway, and apoptosis were also enriched (Figure 2P).

Figure 2 Gene set enrichment analysis. A: Top five gene set enrichment analysis in proliferating myoblasts; B: Neuroactive ligand–receptor interaction; C: Hypertrophic cardiomyopathy; D: DNA replication; E: Cell cycle; F: Cardiac muscle contraction; G: Bubble plot in proliferating myoblasts; H: Ridgeline plot in proliferating myoblasts; I: Top 5 gene set enrichment analysis in differentiated myotubes; J: Viral myocarditis; K: Steroid hormone biosynthesis; L: Hematopoietic cell lineage; M: Focal adhesion; N: Extracellular matrix–receptor interaction; O: Bubble plot in differentiated myotubes; P: Ridgeline plot in differentiated myotubes.

Through GSVA, eight pathways were enriched in proliferating myoblasts (Figure 3A and B). RNA degradation, DNA replication, and mismatch repair were upregulated, and glycosaminoglycan biosynthesis chondroitin sulfate, other glycan degradation, lysosome, glycosaminoglycan biosynthesis heparan sulfate, and steroid biosynthesis were downregulated. Two pathways (steroid hormone biosynthesis and steroid biosynthesis) were enriched in differentiated myotubes (Figure 3C and D).

Figure 3 Gene set variation analysis. A: Volcano plot of proliferating myoblasts; B: Volcano plot of differentiated myotubes; C: Heat map of proliferating myoblasts; D: Heat map of differentiated myotubes.

Identification of ERS-related DEGs

We performed DEG analysis on proliferating myoblasts and differentiated myotubes. We obtained 426 DEGs (188 upregulated and 238 downregulated, Figure 4A and B) and 281 DEGs (135 upregulated and 146 downregulated, Figure 4C and D) from proliferating myoblasts and differentiated myotubes, respectively. Through intersecting with 6893 ERS-related genes, we obtained 227 ERS-related DEGs (Figure 4E).

Figure 4 Endoplasmic reticulum stress-related differentially expressed genes. A: Heat map of proliferating myoblasts; B: Volcano plot of proliferating myoblasts; C: Heat map of differentiated myotubes; D: Volcano plot of differentiated myotubes; E: Endoplasmic reticulum stress-related differentially expressed genes. ER: Endoplasmic reticulum.

Function enrichment analysis

GO terms include biological processes, molecular functions, and cellular components. There were 227 ERS-related DEGs enriched in 875 biological process terms, 103 molecular function terms, and 81 cellular component terms. The results indicated that numerous biological processes were involved in extracellular structure organization, collagen fibril organization, ECM organization, cellular response to external stimulus, response to ketone, cellular response to fatty acid, cellular response to prostaglandin stimulus, response to fatty acid, response to mechanical stimulus, respiratory tube development, cellular response to extracellular stimulus, response to alcohol, and cell–substrate adhesion (Figure 5A). The results indicated that numerous cellular components were involved in the collagen-containing ECM, ECM component, collagen trimer, ER lumen, ER–Golgi intermediate compartment, Golgi-associated vesicle membrane, the complex of collagen trimers, membrane raft, membrane microdomain, membrane region, phagocytic vesicle, neuronal cell body, Golgi-associated vesicle, focal adhesion, cell–substrate adherens junction, cell–substrate junction, postsynaptic endosome, transport vesicle, COPII-coated ER to Golgi transport vesicle, and ER to Golgi transport vesicle membrane (Figure 5B). The results indicated that numerous molecular functions were involved in prostaglandin receptor activity, ECM structural constituent, prostanoid receptor activity, icosanoid receptor activity, growth factor binding, ECM structural constituent conferring tensile strength, platelet-derived growth factor binding, transmembrane receptor protein kinase activity, heat shock protein binding, sulfur compound transmembrane transporter activity, transmembrane receptor protein tyrosine kinase activity, transmembrane-ephrin receptor activity, oxidoreductase activity, ephrin receptor activity, virus receptor activity, and hijacked molecular function (Figure 5C). Through KEGG function enrichment analysis, 26 pathways were significant, such as axon guidance, protein digestion and absorption, focal adhesion, protein processing in the ER, cortisol synthesis and secretion, Fc gamma R-mediated phagocytosis, renin secretion, glutathione metabolism, AMPK signaling pathway, ECM–receptor interaction, DNA replication, calcium signaling pathway, thyroid cancer, aldosterone synthesis, and secretion, lipid and atherosclerosis, P53 signaling pathway, other glycan degradation, biosynthesis of amino acids, ABC transporters, phospholipase D signaling pathway, and steroid biosynthesis (Figure 5D).

Figure 5 Functional enrichment analysis. A: Gene Ontogeny biological processes; B: Gene Ontogeny cellular components; C: Gene Ontogeny molecular function; D: Kyoto Encyclopedia of Genes and Genomes.

PPI analysis

The network consisted of 227 nodes and 416 edges (Figure 6A). We used the NetworkAnalyzer plugin to calculate the degree and combine the score (Figure 6B). We obtained 20 key genes (2 modules with 67 interactions) via the cytoHubba plugin (Table 1 and Figure 6C).

Figure 6 Protein–protein interaction network. A: Protein–protein interaction (PPI) network; B: PPI network by NetworkAnalyzer; C: PPI network of 20 key genes.

Table 1 20 key genes.

Gene	Description	MCC-score	Degree	Closeness	Betweenness
COL3A1	Collagen Type III Alpha 1 Chain	11179	16	64.5	1189.37723
COL6A3	Collagen Type VI Alpha 3 Chain	11168	12	62.03333	567.89507
COL4A1	Collagen Type IV Alpha 1 Chain	11167	12	59.91667	480.86034
COL4A2	Collagen Type IV Alpha 2 Chain	10806	10	56.23333	84.68666
COL11A1	Collagen Type XI Alpha 1 Chain	10800	9	54.5	18.34225
PLOD2	Procollagen-Lysine,2-Oxoglutarate 5-Dioxygenase 2	10326	10	56.88333	112.4573
COL22A1	Collagen Type XXII Alpha 1 Chain	10080	8	51.80952	4.88333
SERPINH1	Serpin Family H Member 1	6001	10	60.3	1308.02191
P4HA1	Prolyl 4-Hydroxylase Subunit Alpha 1	5166	9	55.07857	68.74063
MCM7	Minichromosome Maintenance Complex Component 7	1804	14	68.24286	1418.45764
MCM2	Minichromosome Maintenance Complex Component 2	1801	13	63.94524	784.74991
MSH2	MutS Homolog 2	1645	12	62.22619	646.66407
TIPIN	TIMELESS Interacting Protein	1596	10	54.1631	37.22401
SMC2	Structural Maintenance of Chromosomes 2	1567	10	53.77976	230.10198
POLA2	DNA Polymerase Alpha 2, Accessory Subunit	1560	8	52.6131	4.22778
ESPL1	Extra Spindle Pole Bodies Like 1, Separase	962	10	56.52976	455.86734
POSTN	Periostin	860	15	66.51667	1995.77049
SKP2	S-Phase Kinase Associated Protein 2	771	13	64.5619	1861.813
TRIM32	Tripartite Motif Containing 32	722	8	54.52976	282.93712
SIAH2	Siah E3 Ubiquitin Protein Ligase 2	722	8	57.06905	982.96735

MCC: McCormick.

Network analysis

We obtained 27 TFs and 49 target genes from TRRUST to build the TF–mRNA network (Figure 7A). We obtained 51 miRNAs and 25 target genes from miRWalk to build the miRNA–mRNA network (Figure 7B). We also obtained 59 drugs and 22 target genes from DGIdb to build the drug–mRNA network (Figure 7C).

Figure 7 Networks. A: Transcription factor–mRNA network; B: miRNA–mRNA network; C: Drug–mRNA network.

Correlation analysis of immune infiltration

We demonstrated that memory CD4⁺ T cells accounted for the largest proportion of 22 immune cell types (Figure 8A). Figure 8B showed the distribution of different immune cells in each sample. Moreover, we evaluated the correlation between immune infiltration and each sample (Figure 8C). In 20 key genes, the enrichment degree of each immune cell was different (Figure 8D).

Figure 8 Immune infiltration. A: Histogram of immune infiltration distribution; B: Histogram of immune infiltration sample distribution; C: Heat map of immune infiltration correlation; D: Correlation diagram of 20 key genes.

DISCUSSION

T2DM is a complex metabolic disease driven by interactions among diverse environmental and genetic susceptibilities[39]. Although environmental and epigenetic factors clearly play a contributory role in the pathogenesis of T2DM, genetic factors appear to be the primary contributors to the recent rise in T2DM prevalence[40]. More studies have shown that ERS is involved in T2DM[41]. In the present study, we first explored the potential pathways in proliferating myoblasts and differentiated myotubes, and obtained 227 ERS-related DEGs in T2DM, which may contribute to the occurrence and development of T2DM. Later enrichment analysis, immune infiltration, TF–mRNA network, and miRNA–mRNA network revealed the mechanisms of T2DM, which provided a way for clinical treatment of T2DM. In particular, the drug–mRNA network provided new insights and perspectives into the therapeutic reagents.

In GSEA and GSVA, we confirmed that DNA replication, cell cycle, neuroactive ligand–receptor interaction, glycosaminoglycan biosynthesis heparan sulfate, glycosaminoglycan biosynthesis chondroitin sulfate, glycosaminoglycan degradation, other glycan degradation, lysosome, arachidonic acid metabolism, mismatch repair, metabolism of xenobiotics by cytochrome P450, steroid hormone biosynthesis, focal adhesion, and ECM–receptor interaction, neuroactive ligand–receptor interaction, gap junction, steroid biosynthesis, and cell adhesion molecules were enriched. Moreover, the P53 signaling pathway, VEGF signaling pathway, MAPK signaling pathway, and apoptosis may contribute to T2DM. Previous studies have indicated that these biological processes are related to T2DM[42-44]. SRT2104 enhanced renal SIRT1 expression and activity, deacetylated P53, and activated NRF2 antioxidant signaling, providing remarkable protection against T2DM[45]. The p-ERK/p-JNK/VEGF/PKC signaling pathway may play an important role in pathological T2DM conditions[46]. TREM-2 negatively regulates p38 MAPK-mediated inflammatory response in T2DM[47]. These previous findings are consistent with our findings in this study.

We identified 227 ERS-related DEGs and later function enrichment analysis demonstrated that the enriched biological processes and pathways are highly consistent with the previous GSEA and GSVA results. The immune infiltration analysis revealed that memory CD4⁺ T cells accounted for the largest proportion of 22 immune cell types. T2DM patients are present with self-reactive T cells with a memory phenotype[48]. The memory CD4⁺ T cells develop directly from effector cells and thereby preserve features of their effector precursors are reserved[49]. Depending on the immune context, memory CD4⁺ T cells can contribute to immune protection, pathology, or tissue remodeling[50]. The memory CD4⁺ T cells could act as immunological markers for predicting change in β-cell function in T2DM[51]. TFs recognize specific DNA sequences to control chromatin and transcription, forming a complex system that guides the expression of the genome[52]. Here we obtained 27 TFs, which may contribute to T2DM. MiRNA is a class of endogenous noncoding RNA encoding 19–25 nucleotides, which is involved in the post-transcriptional regulation of genes[53]. Most of them have high sequence conservation, expression timing, and tissue specificity[54]. Recent studies have shown that miRNA is involved in a variety of regulatory pathways, we here identified 51 miRNAs to further explain the mechanisms of T2DM. Importantly, we also established a drug–mRNA network map to provide new ideas and directions for the treatment of T2DM. Immune infiltration plays an important role in the occurrence and development of T2DM[55,56]. The memory CD4⁺ T cells play central roles in immunity in health and disease[57]. We also explored the relationship between immune infiltration and T2DM, and we found that memory CD4⁺ T cells were the most numerous types of immune cells in T2DM. Previous studies indicated that CD4⁺ T cells contribute to the destruction of insulin-producing β-cells in type 1 diabetes mellitus[58,59], which confirmed our results.

This study has some limitations. First, all the results of the analysis were derived from previous data. Despite the efforts we have made in the present, our results still need verification experimentally and clinically. Moreover, the TF–mRNA, miRNA–mRNA, and drug–mRNA networks we built in this study provided some new ideas and insights for the mechanisms and treatment of T2DM. However, this is only the beginning, and more work is still needed in the follow-up.

CONCLUSION

This study revealed ERS-related mechanisms in T2DM, which might contribute to new ideas and insights for the mechanisms and treatment of T2DM.

ARTICLE HIGHLIGHTS

Research background

The endoplasmic reticulum (ER) is closely related to a wide range of cellular functions and is a key component to maintain and restore metabolic health.

Research motivation

Type 2 diabetes mellitus (T2DM) is a serious threat to human health, but knowledge of the ER stress (ERS)-related mechanisms in T2DM is lacking.

Research objectives

Here, we conducted a bioinformatics analysis to identify potential ERS-related mechanisms and crucial biomarkers in T2DM.

Research methods

We conducted gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA) in myoblast and myotube form GSE166502, and obtained the differentially expressed genes (DEGs). After intersecting with ERS-related genes, we obtained ERS-related DEGs. Finally, functional analyses, immune infiltration, and several networks were established.

Research results

Through GSEA and GSVA, we identified several metabolic and immune-related pathways. We obtained 227 ERS-related DEGs and constructed several important networks that help to understand the mechanisms and treatment of T2DM. Finally, memory CD4⁺ T cells accounted for the largest proportion of immune cells.

Research conclusions

This study revealed ERS-related mechanisms in T2DM.

Research perspectives

Our study might contribute to new ideas and insights for the mechanisms and treatment of T2DM.