Differentially Expressed Genes Regulating Glutathione Metabolism, Protein-Folding, and Unfolded Protein Response in Pancreatic β-Cells in Type 2 Diabetes Mellitus

Impaired redox homeostasis in the endoplasmic reticulum (ER) may contribute to proinsulin misfolding and thus to activate the unfolded protein response (UPR) and apoptotic pathways, culminating in pancreatic β-cell loss and type 2 diabetes (T2D). The present study was designed to identify differentially expressed genes (DEGs) encoding enzymes for glutathione metabolism and their impact on the expression levels of genes regulating protein folding and UPR in β-cells of T2D patients. The GEO transcriptome datasets of β-cells of diabetics and non-diabetics, GSE20966 and GSE81608, were analyzed for 142 genes of interest using limma and GREIN software, respectively. Diabetic β-cells showed dataset-specific patterns of DEGs (FDR ≤ 0.05) implicated in the regulation of glutathione metabolism (ANPEP, PGD, IDH2, and CTH), protein-folding (HSP90AB1, HSP90AA1, HSPA1B, HSPA8, BAG3, NDC1, NUP160, RLN1, and RPS19BP1), and unfolded protein response (CREB3L4, ERP27, and BID). The GCLC gene, encoding the catalytic subunit of glutamate–cysteine ligase, the first rate-limiting enzyme of glutathione biosynthesis, was moderately down-regulated in diabetic β-cells from both datasets (p ≤ 0.05). Regression analysis established that genes involved in the de novo synthesis of glutathione, GCLC, GCLM, and GSS affect the expression levels of genes encoding molecular chaperones and those involved in the UPR pathway. This study showed for the first time that diabetic β-cells exhibit alterations in the expression of genes regulating glutathione metabolism, protein-folding, and UPR and provided evidence for the molecular crosstalk between impaired redox homeostasis and abnormal protein folding, underlying ER stress in type 2 diabetes.


Introduction
Diabetes mellitus is one of the fastest-growing global health challenges of the 21st century [1]. Based on the 2021 International Diabetes Federation Diabetes Atlas, 537 million adults aged 20-79 years (a prevalence of 10.5%) are estimated to be living with diabetes mellitus worldwide, and 783 million (a prevalence of 12.2%) are projected to be living with the condition by 2045 [1]. Russia is in second place among all European countries in terms of the number of diabetic patients, 90% of whom are type 2 diabetics (T2D) [1,2].
Type 2 diabetes is known as a chronic multifactorial disease characterized by peripheral insulin resistance, impaired hepatic glucose production regulation, and declining pancreatic beta-cell function attributed to cellular death with unknown etiology [3,4]. Dysfunction of pancreatic beta-cell is primary in T2D and is detected at the earliest disease stage, and the loss of beta cells is responsible for disease progression and complications [5][6][7]. It has been observed that oxidative stress in T2D, caused by a shift in redox homeostasis towards decreased antioxidant defense and excessive production of reactive oxygen and nitrogen species, is strongly associated with dysfunction and apoptosis of beta cells, which are extremely sensitive to oxidative damage [8,9]. Endogenous glutathione deficiency has been identified in T2D and is widely agreed to be the primary cause of oxidative stress and mitochondrial dysfunction linked to disease development and progression [10][11][12][13][14][15][16].
In type 2 diabetes, pancreatic beta-cell dysfunction and apoptosis are key processes attributed to the activation of endoplasmic reticulum stress affecting beta cells through PERK-eIF2α-mediated damage [17,18]. We recently hypothesized [16] that glutathione deficiency may be a major cause of impaired folding of proinsulin, which has been identified in type 2 diabetes mellitus [19][20][21]. The hypothesis was based on the assumption that glutathione is essential for regulating the formation of native disulfide bonds within proteins entering the secretory pathway, and numerous studies have demonstrated that glutathione has a strong influence on disulfide pairing in the tertiary structure of target proteins, including proinsulin [22][23][24][25][26][27][28][29][30]. In ensuring protein folding, oxidized glutathione (GSSG) acts as an oxidant, providing disulfide bond formation, whereas reduced glutathione (GSH) functions as a reducing agent, cleaving mis-bridged disulfide bonds in maturating proteins [23,31]. This means that glutathione deficiency, on the one hand, may be per se a key factor responsible for the formation of unfolded or misfolded proteins in the endoplasmic reticulum (ER), on the other hand, contributes to oxidative stress, causing a vicious circle in disrupting disulfide bond formation and aggregation of misfolded proteins, leading to the activation of unfolded protein response and apoptotic pathways in type 2 diabetes. The unfolded protein response (UPR) is a cellular stress response that is activated as a result of the accumulation of unfolded or misfolded proteins in the ER lumen, with the goal of restoring normal cell function by attenuating protein translation, degrading misfolded proteins, and activating signaling pathways that lead to increased production of molecular chaperones that implement protein folding [32]. Numerous studies have demonstrated that chronic ER stress and associated UPR are characteristic features of type 2 diabetes [17,18,[33][34][35][36][37]. An intriguing fact is that the accumulation of misfolded proinsulin molecules in the ER due to disruption of the formation of disulfide-related complexes is detected early in prediabetes and subsequently exacerbated, leading to apoptosis and beta-cell dysfunction [19][20][21]. Our recent study discovered that polymorphisms in genes encoding glutathione metabolizing enzymes, such as glutathione synthase (GSS) and gamma-glutamyltransferase 7 (GGT7), are associated with susceptibility to type 2 diabetes [16]. Other genes involved in glutathione metabolism have also been found to influence the risk of type 2 diabetes [38][39][40][41][42][43]. Interestingly, these polymorphisms were found to correlate with the tissue-specific expression of genes involved in the UPR pathway (EIF2S2, EIF6, UQCC1, PIGU, ITCH, TRPC4AP, EDEM2, TP53INP2, and MAP1LC3A), glucose-stimulated insulin secretion, and glucose uptake by cells (MYH7B, CPNE1, TRPC4AP, and NCOA6) [16], providing additional evidence for the existence of shared molecular mechanisms involved in the regulation of glutathione metabolism and protein folding, and linking them to the pathogenesis of type 2 diabetes. No studies have so far been performed to investigate the expression levels of genes regulating protein folding, UPR, and glutathione metabolism in pancreatic beta cells in type 2 diabetes. The objective of this study was to examine the transcriptome profile and co-expression of genes regulating glutathione metabolism, protein-folding, and unfolded protein response in the pancreatic beta cells of type 2 diabetics and non-diabetic individuals.

Correlations between Transcription Levels of Differentially Expressed Genes in β-Cells
A correlation analysis was performed to investigate the relationship between expression levels of the studied genes in pancreatic beta cells of T2D patients and non-diabetic controls in both datasets. Correlation matrices of diabetics and controls are shown in Supplementary Table S2. Partial correlations between gene expression levels (p ≤ 0.05) are highlighted in red (Supplementary Table S2). Multiple correlations between expression levels genes of glutathione metabolism, molecular chaperons, and unfolding protein response pathway were observed. The correlations between differentially expressed genes in patients with type 2 diabetes and non-diabetics are of interest. Figure 1 shows a schematic presentation of correlations between DEGs in T2D patients (A) and non-diabetics (B) from Dataset 1. In diabetics, ANPEP correlated with ERP27     In each dataset, cluster analysis of correlation matrices was used to depict subgroups of inter-correlated genes in diabetics and non-diabetic subjects.  In each dataset, cluster analysis of correlation matrices was used to depict subgroups of inter-correlated genes in diabetics and non-diabetic subjects. Figures 3 and 4 illustrate tree diagrams depicting clusters of inter-correlated expression between genes in T2D patients and non-diabetic controls from Datasets 1 and 2, respectively. Tree networks differed across the T2D samples and datasets. As can be seen from the tree diagram of T2D patients ( Figure 3), DEGs such as ANPEP and ERN1, G6PD and DNAJB11, and GCLC and TNFRSF21 are present in several independent clusters. In the T2D patients in Dataset 2 ( Figure 4), DEGs were also present in independent gene clusters. Interestingly, one such cluster comprised several inter-correlated DEGs such as HSP90AB1, HSP90AA1, HYOU1, HSPA8, and ST13 (NFE2L2 and GSK3B are also present in this cluster). Moreover, IDH2 and CTH genes were united by the shared cluster. The remaining DEGs were present separately in independent clusters.

The Impact of Glutathione Biosynthesis Genes on the Expression Levels of Genes for Protein Folding and Unfolded Protein Response in Pancreatic β-Cells
Since the level of intracellular glutathione is of great importance for the formation of disulfide bonds during the formation of the three-dimensional structure of proteins in the endoplasmic reticulum [23,31], it is very important to investigate whether the expression of genes encoding key enzymes of glutathione biosynthesis, such as GCLC, GCLM, and GSS, affects the transcriptional activity of genes encoding molecular chaperones and regulators of the unfolded protein response. Pursuing this interest, linear regression analysis was applied to investigate the impact of the expression of glutathione biosynthesis genes on the expression levels of genes regulating protein folding and the unfolded protein response in pancreatic beta cells of patients with type 2 diabetes and non-diabetic individuals. The beta regression coefficients estimated for the impact of GCLC, GCLM, and GSS genes on the expression of genes regulating protein folding and the UPR pathway in Datasets 1 and 2 are shown in Supplementary Tables S3 and S4, respectively. A schematic representation of the statistically significant regression coefficients at p ≤ 0.05 observed in Datasets 1 and 2 is shown in Figures 5 and 6, respectively. tree diagrams depicting clusters of inter-correlated expression between genes in T2D patients and non-diabetic controls from Datasets 1 and 2, respectively. Tree networks differed across the T2D samples and datasets. As can be seen from the tree diagram of T2D patients (Figure 3), DEGs such as ANPEP and ERN1, G6PD and DNAJB11, and GCLC and TNFRSF21 are present in several independent clusters. In the T2D patients in Dataset 2 ( Figure 4), DEGs were also present in independent gene clusters. Interestingly, one such cluster comprised several inter-correlated DEGs such as HSP90AB1, HSP90AA1, HYOU1, HSPA8, and ST13 (NFE2L2 and GSK3B are also present in this cluster). Moreover, IDH2 and CTH genes were united by the shared cluster. The remaining DEGs were present separately in independent clusters.    uals. The beta regression coefficients estimated for the impact of GCLC, GCLM, an genes on the expression of genes regulating protein folding and the UPR pathway tasets 1 and 2 are shown in Supplementary Tables S3 and S4, respectively. A sche representation of the statistically significant regression coefficients at p ≤ 0.05 obser Datasets 1 and 2 is shown in Figures 5 and 6, respectively.  . Schematic presentation of the impact of genes encoding enzymes for glutathione biosynthesis on the expression levels of genes for protein folding and unfolded protein response in pancreatic β-cells of diabetic patients (A) and non-diabetics (B) from Dataset 1. Schemes were constructed using beta regression coefficients for the impact of GCLC, GCLM, and GSS genes on the expression levels of genes regulating protein folding and the unfolded protein response. Glutathione biosynthesis genes are shown in blue, genes for protein folding are shown in green, and genes of the UPR pathway are shown in red (the impact on these genes was found at FDR ≤ 0.05; genes which were impacted at p ≤ 0.05 and FDR ≥ 0.05 are shown in grey). Red lines indicate positive impacts, and blue lines indicate negative impacts. The line thickness reflects the impact strength. A yEd graph editor version 3.22 was used to draw the figure.
expression levels of genes regulating protein folding and the unfolded protein response. Gluta one biosynthesis genes are shown in blue, genes for protein folding are shown in green, and ge of the UPR pathway are shown in red (the impact on these genes was found at FDR ≤ 0.05; ge which were impacted at p ≤ 0.05 and FDR ≥ 0.05 are shown in grey). Red lines indicate pos impacts, and blue lines indicate negative impacts. The line thickness reflects the impact strengt yEd graph editor version 3.22 was used to draw the figure. Figure 6. Schematic presentation of the impact of genes encoding enzymes for glutathione bio thesis on the expression levels of genes for protein folding and unfolded protein response in p creatic β-cells of diabetic patients (A) and non-diabetics (B) from Dataset 2. Schemes were structed using beta regression coefficients for the impact of GCLC, GCLM, and GSS genes on expression levels of genes regulating protein folding and the unfolded protein respo Figure 6. Schematic presentation of the impact of genes encoding enzymes for glutathione biosynthesis on the expression levels of genes for protein folding and unfolded protein response in pancreatic β-cells of diabetic patients (A) and non-diabetics (B) from Dataset 2. Schemes were constructed using beta regression coefficients for the impact of GCLC, GCLM, and GSS genes on the expression levels of genes regulating protein folding and the unfolded protein response. Glutathione biosynthesis genes are shown in blue, genes for protein folding are shown in green, and genes of the UPR pathway are shown in red (the impact on these genes was found at FDR ≤ 0.05; genes which were impacted at p ≤ 0.05 and FDR ≥ 0.05 are shown in grey). Red lines indicate positive impacts, and blue lines indicate negative impacts. The line thickness reflects the impact strength. A yEd graph editor version 3.22 was used to draw the figure.

Discussion
Despite decades of intensive basic and clinical research, the primary molecular mechanisms underlying the etiology and pathogenesis of type 2 diabetes mellitus remain far from being completed. The question of what causes the loss of pancreatic beta cells, a key pathological process underlying the development and progression of T2D, remains unanswered. It is widely accepted that chronic endoplasmic reticulum stress and activation of the response to unfolded proteins initiate apoptotic pathways that lead to beta-cell death [17][18][19]34,37,44,45]. If the reason for UPR activation in monogenic forms of diabetes is clear-mutations cause structural disorders of the proinsulin molecule, which cannot be folded properly-then, in more common variety type 2 diabetes, the reasons responsible for chronic activation of ER stress remain unknown. We have recently proposed [16] that the deficiency of intracellular glutathione, which is well known to be important in the formation of disulfide bridges in the native tertiary structure of a protein, could be the root cause of proinsulin misfolding in type 2 diabetes. We hypothesize that defective proinsulin folding, a condition that has been linked to type 2 diabetes in numerous studies [19][20][21][46][47][48], may be caused by glutathione deficiency in pancreatic beta cells due to its decreased biosynthesis and/or increased depletion.
It is well known that misfolded or unfolded proteins accumulate in the lumen of the ER due to impaired folding, triggering a cell response to unfolded proteins called the unfolded protein response [32]. UPR is regulated by three ER transmembrane proteins: eukaryotic translation initiation factor 2-alpha kinase 3 (PERK), activating transcription factor 6 (ATF6), and inositol-requiring enzyme 1 (IRE1). PERK, through eIF2α phosphorylation, turns on ATF4, a key transducer that controls the transcription of genes involved in autophagy, apoptosis, amino acid metabolism, and antioxidant responses. ATF6 is transported to the Golgi apparatus, where, as a result of interaction with sites 1 (S1P) and 2 (S2P) proteases, it releases its fragment of the cytosolic domain, ATF6f, which activates genes promoting ER-associated degradation (ERAD) and protein folding. IRE1 and its subsequent auto-transphosphorylation activate the transcription factor XBP1, which controls the transcriptional activity of genes encoding genes involved in protein folding and their quality control as well as ER-associated proteasomal degradation [32,49].
To date, no studies have been conducted to specifically examined the expression of key genes involved in protein folding and the unfolded protein response in pancreatic β-cells of patients with type 2 diabetes. In the present study, we used two independent transcriptomic datasets to look at the expression levels of genes involved in protein folding, UPR, and glutathione metabolism in pancreatic beta cells from diabetic and non-diabetic patients. Compared to non-diabetics, T2D patients were found to have decreased expression of genes regulating protein folding such as BAG3, NDC1, HSPA7, HSPB2, RLN1, and TNFRSF21 in Dataset 1, and HSP90AB1, HSPA1B, HSP90AA1, HSPA8, NUP160, ST13, and RPS19BP1 in Dataset 2. Furthermore, several protein folding regulating genes such as RPS19BP1, RAE1, MTOR, POM121C, and GFER in Dataset 1 and NUP35, NUP58, and NUP107 in Dataset 2 have been found to be up-regulated in beta cells from T2D patients when compared to non-diabetics. Induction of chaperones (e.g., from the Hsp70 and Hsp27 families) is known to enhance the adaptive capacity of the endoplasmic reticulum, reduce ER stress, and restore glucose homeostasis, as it has been demonstrated in a mouse model of type 2 diabetes [50]. BAG3 is a member of the BAG family of co-chaperone proteins such as HSP70 and HSC70 [51]. BAG3 is involved in chaperone-assisted selective autophagy, apoptosis, cell adhesion, and cytoskeleton remodeling [52][53][54], suggesting a protective role for this chaperone against type 2 diabetes. Interestingly, the down-regulation of BAG3 in beta cells was found to increase insulin secretion in response to glucose stimulation [55], and in this context, the decreased levels of BAG3 mRNA in diabetics may be required for facilitating insulin secretion. No associations of type 2 diabetes with the expression level of heat shock proteins such as HSPA1B, HSPA7, HSPB2, HSP90AB1, HSPA1B, HSP90AA1, and HSPA8 have been reported in the literature.
HSPA1B is a member of the human HSP70 family, whose expression was also increased in the beta cells of diabetics (Dataset 2). Several studies have demonstrated the involvement of the Hsp70 chaperone family in insulin resistance in type 2 diabetes. HSPA1B is a molecular chaperone of the Hsp70 family implicated in a wide variety of cellular processes, including regulation of proteostasis through mechanisms such as folding and transport of newly synthesized polypeptides, quality control of proteins in the ER by ensuring the correct folding of proteins, refolding of misfolded proteins, and controlling the targeting of proteins for subsequent degradation [56][57][58]. Increased serum levels of HSP70 chaperones have been revealed in patients with type 2 diabetes, and they were correlated with disease duration [57]. Polymorphisms in the HSPA1B gene have been found to be associated with the severity of diabetic foot ulceration [59] and chronic heart failure [60].
HSP90AA1 and HSP90AB1 are members of the Hsp90 family of chaperones comprising 1-2% of cellular protein in unstressed cells and up to 4-6% in stressed cells [61]. Expression of HSP90AA1 and HSP90AB1 genes are differentially regulated at the transcriptional level [61]. These chaperones are primarily involved in signal transduction via transcription factors that initiate gene expression, kinases that transmit information by post-translationally modifying other proteins, and E3-ligases that target proteins for proteasomal degradation [61]. It has been hypothesized that biological processes associated with HSP90AB1 are more relevant to maintaining viability, whereas the processes associated with HSP90AA1 are more related to adaptation to stress or other specialized functions [61]. Experiments on mouse models of diabetes have shown that treatment of mice with a pan-Hsp90 inhibitor [62] and also knockdown of HSP90ab1 [63] have been found to improve glucose homeostasis and insulin sensitivity, suggesting a role of these chaperones in the pathogenesis of type diabetes.
HSPA7 is a putative heat shock 70 kDa protein 7, a pseudogene whose functions are unknown. However, a recent study reported that HSPA7 is a member of long noncoding RNAs (lncRNAs) identified in human atherosclerotic plaques that play a role in the inflammatory transition of vascular smooth muscle cells stimulated by oxidized low-density lipoproteins [64]. HSPA8 (also known as HSC70) is also a member of the Hsp70 family, which, like HSPA1B, is involved in a wide variety of cellular processes such as protection of the proteome from stress, folding and transport of newly synthesized polypeptides, chaperone-mediated autophagy, activation of proteolysis of misfolded proteins and the formation and dissociation of protein complexes [56][57][58]. It has been revealed that patients with type 2 diabetes exhibit higher levels of the stress-related protein complex HSPA8/Hsp90/CSK2α, along with higher platelet aggregation and thrombogenicity than nondiabetic subjects [65]. HSPB2 is a heat shock protein beta-2 which may regulate dystrophia myotonia protein kinase (DMPK) [66]. HSPB2 may also play a role in the regulation of metabolic and mitochondrial functions [67]. It has been experimentally shown that HspB2 deficiency possesses a protective effect form diet-induced glucose intolerance [68].
We also found changes in the expression of nucleoporin genes such as NDC1, NUP35, NUP58, NUP107, NUP160, and POM121C in pancreatic beta cells of patients with type 2 diabetes mellitus. They are members of the nucleoporin family, structural components of the nuclear pore complex (NPC) [69], which allows for access to the nucleus and controls how proteins and RNA are transported over the nuclear envelope, thereby providing bidirectional molecular movements between the nucleus and cytoplasm [70]. NDC1 (Nuclear-Division-Cycle 1) is a transmembrane nucleoporin (also known as Transmembrane Protein 48 or TMEM48), which plays a key role in de novo assembly and insertion of NPC in the nuclear envelope [69]. Altered expression of NDC1 is observed in ischemic and dilated cardiomyopathy as well as in numerous malignancies [71]. NUP35 also functions as a structural component of NPC, where it plays a role in docking or interacting with partners for transiently associated nuclear transport factors [72]. NUP58 is a nucleoporin, a component of NPC that is involved in nucleocytoplasmic trafficking [73]. NUP107 is a nucleoporin required for the assembly of peripheral proteins into the NPC [74]. NUP160 is a component of NPC involved in the poly (A)+ RNA transport [75]. Interestingly, NUP160 expression was found to be upregulated, which is associated with the inhibition of autophagy and increased inflammatory response in mice with diabetic nephropathy [76]. POM121C is also an essential component of the nuclear pore complex. Interestingly, expression of POM121C in beta cells of T2D patients was up-regulated in Dataset 1, and down-regulated in Dataset 2, suggesting its role in diabetes. POM121C is considered a candidate gene for fasting serum insulin (FSI), and the FSI-increasing allele at rs1167800 is associated with lower POM121C expression, which is associated with systemic insulin sensitivity, adipocyte insulin sensitivity, and adipose hyperplasia [77]. POM121C may contribute to insulin resistance in type 2 diabetes by stimulating adipogenesis and increasing the sensitivity of adipocytes to insulin [77]. Thus, as discussed above, changes in the expression levels of most nucleoporins have not been previously established in type 2 diabetes. However, some studies have demonstrated the importance of nucleoporins in the development of diabetes and its complications.
ST13 is an Hsp70-interacting protein that may contribute to the interaction of HSC70 (chaperone HSPA8) with various target proteins (https://www.uniprot.org/, accession number P50502, accessed on 7 July 2023). RLN1 is a member of relaxins, endocrine, and autocrine/paracrine hormones belonging to the insulin gene superfamily (https://www. uniprot.org/, accession number P04808, accessed on 7 July 2023). TNFRSF21 is a tumor necrosis factor receptor superfamily member 21 that can promote apoptosis, possibly through the activation of NF-kappa-B and BAX pathways and by the release of cytochrome c [78]. TN-FRSF21 has recently been shown to play a pathogenic role in heart disease in patients with type 2 diabetes [79]. Increased expression of TNFRSF21 is considered a urinary biomarker of diabetic nephropathy progression [80]. RAE1 (ribonucleic acid export 1) is an mRNA export factor involved in nucleocytoplasmic transport by attaching cytoplasmic microribonucleoproteins to the cytoskeleton [81]. MTOR is a serine/threonine protein kinase that controls cellular metabolism, growth, and survival in response to hormones, growth factors, nutrition, energy, and stress signals (https://www.uniprot.org/, accession number P42345, accessed on 7 September 2023). The mechanistic target of rapamycin (mTOR) signaling is an important intracellular pathway that integrates local nutrition and systemic energy status at the organismal and cellular levels. mTOR signaling dysregulation is linked to a variety of diseases, including type 2 diabetes [82]. MTOR is known to regulate the adaptation of beta cells to blood glucose, and dysregulation in mTOR signaling may facilitate the development of type 2 diabetes or insulin resistance. It has been revealed that short-term mTORC1 activation enhances beta cell mass and improves glucose metabolism, whereas long-term mTORC1 activation deteriorates beta cell mass and function, as observed in type 2 diabetic beta cells [82,83]. RPS19BP1 is another gene that showed differential expression in the beta cells of diabetics in both analyzed datasets. RPS19BP1 is a ribosomal protein S19 binding protein 1, which is a part of the small subunit processome that directly regulates SIRT1 [84]. It has been demonstrated experimentally in mice that the activation of SirT1 optimizes the organism for metabolic adaptation to insulin resistance by improving hepatic insulin sensitivity and decreasing whole-body energy requirements [85]. This finding suggests that RPS19BP1-mediated activation of SIRT1 plays a role in the pathogenesis of type 2 diabetes. Increased expression of GFER in pancreatic beta cells in T2D patients was found in Dataset 1. The GFER gene (growth factor, augmenter of liver regeneration) encodes a protein, a part of the CHCHD4:GFER complex, also known as a disulfide relay machinery, which catalyzes the oxidation of cysteine residues in precursor proteins in the mitochondrial intermembrane space to form disulfide bonds, determining protein folding and stability [86][87][88]. GFER maintains mitochondrial ROS levels for optimal functioning of complexes III and IV of the electron transport chain through glutathionation and promotes the formation of disulfide bonds in the CHCHD4 chaperone molecule, and polymorphism of the GFER has been found to be associated with the risk of type 2 diabetes [89].
Genes regulating the unfolded protein response that are differentially expressed in the beta cells of patients with type 2 diabetes are of great interest. In particular, In particular, CREB3L2, CREB3L4, and ERP27 were up-regulated, whereas SSR1, DNAJB11, and HER-PUD1 were down-regulated in pancreatic β cells in diabetics when analyzed in Dataset 1. In Dataset 2, decreased expression of BID, DDIT3, HYOU1, MAPK8, and HSPA5, along with increased expression of EIF2S3, were found in T2D samples compared to non-diabetic controls. CREB3L2 and CREB3L4 are ER-localized proteins belonging to the bZIP family and are transcriptional activators that play roles in the unfolded protein response as regulators of cell secretory capacity and cell-specific cargo [90]. Importantly, transcription factors of the CREB3 family were found to regulate high-fat diet-induced obesity and energy metabolism, which are conditions associated with type 2 diabetes [91]. It has been revealed that Creb3l4 knockout mice exhibit glucose tolerance and decreased insulin sensitivity [92], suggesting a role of CREB3L4 in both obesity and type 2 diabetes. ERP27 (a noncatalytic member of the protein disulfide isomerase) is an endoplasmic reticulum-resident protein that plays a role in the unfolded stress response by binding to unfolded proteins and recruiting PDIA3 protein disulfide isomerase to unfolded substrates (ERP27 has a unique affinity for unfolded proteins and may attract protein disulfide I to substrates that are unfolded) [93] and expression of ERP27 in beta cells was up-regulated in diabetic patients [94]. EIF2S3 (eukaryotic translation initiation factor 2 subunit 3) is an eIF2 complex component that functions in the early stages of protein synthesis by building a ternary complex with GTP and an initiator tRNA (https://www.uniprot.org/, accession number P41091, accessed on 7 October 2023). EIF2S3 may play a critical role in human hypothalamo-pituitary development and function and regulation of glucose metabolism [95], and mutations in the EIF2S3 gene may contribute to neonatal hypoglycemia, followed by early onset diabetes and hypopituitarism [96]. The SSR1 (signal sequence receptor subunit 1) gene encodes translocon-associated protein (TRAP) subunit alpha, being a part of a complex that binds calcium to the ER membrane and thereby regulating the retention of ER-resident proteins and may function as a membrane-bound chaperone facilitating folding of translocated proteins (https://www.uniprot.org/, accession number P43307, accessed on 7 October 2023). It is important to note that common single-nucleotide polymorphisms in the human TRAP gene are linked to type 2 diabetes susceptibility, and the associated pancreatic cell dysfunction indicates that impaired preproinsulin translocation and proinsulin trafficking may play a role in T2D pathogenesis [97]. HSPA5 is an endoplasmic reticulum chaperone BiP with a variety of biological functions, including a key role in protein folding and quality control in the ER lumen [98], post-translational transport of small presecretory proteins across the ER [99], provides the correct folding of proteins, and degradation of misfolded proteins [100], and acts as a key repressor of the ERN1/IRE1-mediated arm of the UPR pathway [100]. HSPA5 is also known as GRP78, whose serum concentration was found to be significantly higher in T2D patients than in healthy controls and was positively correlated with HbA1c levels in diabetics [101]. Experiments in mice suggest that decreased GRP78 expression in the liver may induce resistance to insulin by inhibiting AKT activation and may play an important role in the development of type 2 diabetes [102]. The overproduction of GRP78 in pancreatic beta cells was found to protect mice against high-fat diet-induced diabetes [103]. DNAJB11 is known as a co-chaperone for HSPA5 and is required for the proper folding, trafficking, or degradation of proteins [104]. DNAJB11 binds to unfolded proteins, which are ERAD substrates, as well as to nascent unfolded peptide chains [105]. This co-chaperone may help recruit HSPA5 and other chaperones to the substrate and stimulate the ATPase activity of HSPA5 [106]. Notably, decreased expression of DNAJB1 in pancreatic islets has also been observed in mice with experimental T2D [107]. HERPUD1 (homocysteine inducible ER protein with ubiquitin-like domain 1) is a component of the ER quality control system involved in ubiquitin-dependent ERAD degradation of misfolded proteins [108]. Herpud1 has been found to regulate insulin secretion in mice [109]. BID (BH3 interacting domain death agonist (BID) is a protein that induces caspase and apoptosis pathways [110] and allows the release of cytochrome c by mitochondria [111]. An experimental study on mouse islets demonstrated that Bid is essential for death receptor-induced apoptosis of islets [112], suggesting its role in beta cell dysfunction in type 2 diabetes. DDIT3 (DNA damage-inducible transcript 3 protein), also known as CHOP, represents a multifunctional transcription factor playing an essential role in response to a variety of cell stresses and induces cell cycle arrest and apoptosis in response to ER stress [113,114]. It has been revealed that depletion of Ddit3 in β cells alleviates ER stress in mice [115]. Similar results were obtained in a study by Oyadomari with colleagues [116]. HYOU1 is a hypoxia-upregulated protein 1, ER chaperone that plays a key role in cytoprotection under oxygen deprivation and also as a molecular chaperone participating in protein folding [117]. Up-regulation of HYOU1 and HSPA5 was found in the tubular epithelia of patients with diabetic nephropathy compared with control kidneys [118]. MAPK8 is a mitogen-activated protein kinase 8 that is responsive to various cellular stress stimuli and regulates a wide range of biological processes, including cellular proliferation, differentiation, migration, transformation, autophagy, and apoptosis [119]. Like other MAPKs, MAPK8 may affect insulin signaling and, therefore, play a role in the pathophysiology of type 2 diabetes [120,121]. Moreover, MAPKs have been observed to phosphorylate the HSF1 chaperone, inhibiting HSF1-induced transcription [122], which may result in decreased downstream chaperone activation and protein folding. A decrease in MAPK8 gene expression in type 2 diabetes pancreatic beta cells could be due to, on the one hand, a shift in insulin signaling and glucose transport into the cell and, on the other hand, a decrease in chaperone activation and protein folding. Undoubtedly, all the above assumptions require experimental confirmation. Thus, the results obtained from the analysis of both datasets, despite the observation of dataset-specific DEGs, showed that in pancreatic beta cells, there were changes in the expression of genes that regulate protein folding and those involved in the regulation of UPR. These findings are in line with those of studies that observed abnormal folding of proinsulin in patients with type 2 diabetes [19][20][21].
Given the relevance of defective glutathione metabolism in the pathophysiology of type 2 diabetes [10][11][12]15,16], genes encoding enzymes involved in glutathione metabolism regulation may contribute to disease development and are attractive targets for diabetes research. This sparked an interest in studying the expression levels of these genes in pancreatic beta cells from type 2 diabetes patients. We also wanted to know if the expression levels of glutathione metabolism genes are related to the expression levels of genes that govern protein folding and those involved in the unfolded protein response. The current study identified for the first time differentially expressed genes involved in glutathione metabolism in type 2 diabetes pancreatic beta cells. In particular, we found a decreased expression of the GCLC gene in the samples from patients with type 2 diabetes compared with non-diabetic controls in both datasets. The GCLC gene encodes the catalytic subunit of glutamate-cysteine ligase, the first-rate limiting enzyme of glutathione biosynthesis [123], and therefore the decreased expression of GCLC in T2D patients may contribute to a diminished capacity of beta cells to synthesize glutathione. This assumption may be supported by our recent findings that genetic variants in glutamate-cysteine ligase confer protection against type 2 diabetes due to their positive effect on glutathione levels [38]. Moreover, the levels of ANPEP were increased in the beta cells of diabetics as compared with non-diabetic samples, as has been reported previously by Marselli and co-authors [94] and further replicated in a study by Locke and colleagues [124]. Interestingly, transforming growth factor-beta1, whose plasma concentrations are increased in type 2 diabetes and diabetic renal disease [125,126], was found to increase both the expression and activity of ANPEP in a time-and concentration-dependent manner [127]. Moreover, polymorphisms of the ANPEP gene were found to be associated with the risk of T2D [128,129]. Thus, the above findings clearly show that the ANPEP gene can be involved in the pathogenesis of type 2 diabetes. The mechanisms by which ANPEP contributes to T2D susceptibility remain unclear. ANPEP is alanyl aminopeptidase, an enzyme with broad substrate specificity catalyzing the hydrolysis of peptide bonds with the removal of amino acids from the amino terminus of proteins and peptides. It is important to note ANPEP hydrolyzes peptide L-cysteinylglycine to cysteine and glycine-substrates for de novo biosynthesis of GSH [130,131]. As we recently proposed [16] that the increased levels of ANPEP mRNA in pancreatic beta cells and plasma levels of GGT1 (gamma-glutamyltransferase) in T2D patients may be indicative of an adaptive mechanism in response to endogenous glutathione deficiency: up-regulation of these membrane-bound enzymes in diabetes is necessary for the sequential production of three amino acids (AA)-precursors for intracellular GSH biosynthesis. In particular, GGT1 catabolizes extracellular GSH into L-glutamate (first AA) and L-cysteinylglycine [100,101], then the latter is hydrolyzed by ANPEP with the formation of cysteine (second AA) and glycine (third AA) [130,131], which are then transported into the cell for GSH biosynthesis. This hypothesis is attractive in that it explains several seemingly unrelated biochemical abnormalities that are often detected in type 2 diabetes mellitus and, at the same time, explains the commonality of their occurrence as a result of an endogenous deficiency of glutathione.
In addition, we observed decreased levels of PGD in pancreatic beta cells of diabetic patients compared with non-diabetic controls in both datasets. Furthermore, decreased expression of IDH1 and IDH2, increased expression of CTH in Dataset 2, and increased expression of G6PD in Dataset 1 were observed in T2D samples, in contrast to the controls. The PGD gene encodes phosphogluconate dehydrogenase, an enzyme catalyzing the oxidative decarboxylation of 6-phosphogluconate to ribulose 5-phosphate and CO 2 , a reaction accompanying the reduction of NADP to NADPH (https://www.uniprot.org/, accession number P52209, accessed on 7 October 2023). IDH1 (isocitrate dehydrogenase (NADP) cytoplasmic) is an enzyme that catalyzes the NADP+-dependent oxidative decarboxylation of isocitrate to 2-ketoglutarate with generation NADPH [132]. IDH2 is a mitochondrial isocitrate dehydrogenase (NADP), an enzyme associated with the pyruvate dehydrogenase complex and a major source of NADPH, playing a role in energy production [133]. Knockdown of the mitochondrial isocitrate dehydrogenase enzyme in pancreatic beta cells inhibits insulin secretion [134]. Interestingly, transcriptional up-regulation of IDH2 expression has been found to be associated with improved mitochondrial function in endothelial progenitor cells in diabetic patients [135]. A recent study observed that elevated levels of IDH2 in the mouse liver under overnutrition contribute to elevated gluconeogenesis and glycogen synthesis [136]. Thus, it can be assumed that deficiency of IDH1 and IDH2 in pancreatic beta cells may be associated with a decrease in the formation of G6PD, which is necessary for the reduction of GSSG to GSH. G6PD (glucose-6-phosphate dehydrogenase) is an enzyme of the pentose-phosphate shunt, representing an alternative to the glycolysis pathway in the dissimilation of carbohydrates, with the primary function to supply NADPH for many biosynthesis pathways, including the reduction of oxidized glutathione [137,138]. A systematic review and meta-analysis showed that G6PD deficiency may be a risk factor for diabetes [139]. We can hypothesize that increased expression of G6PD in beta cells may favor the production of NADPH, which is necessary for the reduction of oxidized glutathione. CTH (cystathionine gamma-lyase) is an enzyme that catalyzes the last step in the transsulfuration pathway from L-methionine to L-cysteine, which is utilized for glutathione biosynthesis [140]. CTH, an enzyme of the transsulfuration pathway supplying cysteine for GSH synthesis, whose expression level in beta cells in T2D patients was higher than that in non-diabetic controls (dataset 2). CTH is known to promote PERK activity, leading to attenuation of protein translation during the response to ER stress [141]. It can be assumed that the increased expression of the CTH, ANPEP, and G6PD genes in T2D samples compared with non-diabetic ones appear to be adaptive and aimed at replenishing glutathione deficiency in beta cells of diabetics: CTH and ANPEP supply amino acids such as L-cysteine (CTH and ANPEP) and glycine (ANPEP) for de novo biosynthesis of GSH, whereas G6PD generates NADPH for the reduction of GSSG to GSH. The correlations between genes involved in glutathione metabolism confirm the preceding assumptions about potential disruptions in glutathione synthesis in pancreatic beta cells in type 2 diabetes (Figures 1 and 2).
We also found correlations in expression level between genes for glutathione metabolism and genes regulating protein folding and unfolded protein response. In particular, in T2D patients from Dataset 1, ANPEP correlated with ERP27 and RLN1, G6PD correlated with PGD and SSR1, PGD correlated with SSR1 and GFER, HSPA7 correlated with HSPB2, HSPB2 correlated with NDC1, RAE1 correlated with HERPUD1, MTOR correlated with GCLC and TNFRSF21, POM121C correlated with HERPUD1. In diabetics from Dataset 2, GCLC correlated with HSP90AB1, PGD correlated with POM121C, IDH2 correlated with CTH, BID, NUP35 and NUP107, CTH correlated with NUP107 and NUP35, HSP90AB1 correlated with HYOU1, HSP90AA1 and NUP160, HSPA1B correlated with CTH, IDH2 and NUP107, HSP90AA1 correlated with HYOU1, NUP35 correlated with NUP107, NUP58 correlated with GCLC, ST13 and NUP107, NUP107 correlated with POM121C, RPS19BP1 and BID, NUP160 correlated with IDH2 and DDIT3, ST13 correlated with HYOU1, HSP90AA1 and HSPA8, RPS19BP1 correlated with MAPK8. An interesting finding was seen in both datasets: the expression of the insulin gene in β-cells of T2D patients positively correlated with the expression of CALR, a calcium-binding chaperone that promotes folding, oligomeric assembly, and quality control in the ER (https://www.uniprot.org/, accession number P27797, accessed on 7 December 2023), and this cor-relation was not observed in non-diabetic samples. Furthermore, in non-diabetic β-cells, insulin expression correlated positively with cyclic AMP-dependent transcription factors involved in the main arms of the unfolded protein response: ATF4 in Dataset 1 and ATF3 and ATF6 in Dataset 2. These findings suggest that the expression of insulin in type 2 diabetes is co-regulated by gene networks governing protein folding and the unfolded protein response in beta cells. Cluster analysis revealed groups of inter-correlated genes, demonstrating that the differentially expressed genes in beta cells in type 2 diabetes may have a relatively independent regulation. For instance, ANPEP showed a positive correlation with ERP27 (both genes were strongly up-regulated in T2D samples), and these genes were united by a shared relatively independent gene cluster ( Figure 3A).
Taken together, these findings suggest co-expression of genes regulating glutathione metabolism and genes involved in protein folding and the UPR pathway. Such correlations were expected because it is known that glutathione is required to regulate the formation of native disulfide bonds within proteins entering the secretory pathway [22,23,142]; therefore, cross-talk between glutathione and protein folding pathways is native. It was not surprising that correlations in gene expression levels differed in Datasets 1 and 2. We hypothesize that differences in gene signatures in beta cells of diabetics between the datasets may be attributed to a variety of reasons, including the degree of glutathione deficiency determining the redox state due to insufficient intake of amino acid GSH precursors into the cell, the impact of polymorphisms at genes encoding enzymes of glutathione metabolism, and genes whose products process protein folding in the endoplasmic reticulum and regulate the unfolded protein response.
The unfolded protein response is a universal molecular mechanism that maintains proteostasis under stress in both the endoplasmic reticulum and mitochondria, which are the primary energy metabolism and protein biosynthesis centers in cells [36]. In type 2 diabetes, dysregulated UPR signaling occurs in both the endoplasmic reticulum and the mitochondria, which are considered to crosstalk via the PERK signaling pathway [36]. In pancreatic β-cells, mitochondrial activity is essential for glucose-stimulated insulin secretion and mitochondrial dysfunction to the pathophysiology of type 2 diabetes [143]. According to a recent study, mitochondrial metabolic activity is communicated to the ER via the relay of redox metabolites, and increasing in the oxidized glutathione may influence oxidative protein maturation and protect from ER stress [138]. Glutathione is known to maintain the proper mitochondrial redox environment by acting directly or as a cofactor in reactions catalyzed by other mitochondrial enzymes, as well as preventing oxidative modifications that can lead to mitochondrial dysfunction [144]. Given the foregoing, it is reasonable to assume that metabolic pathways like glutathione metabolism, protein folding, and unfolded protein response are coordinately regulated and that identifying these common pathways will allow us to better understand the nature of metabolic disorders underlying type 2 diabetes. Biological processes underpinning the link between glutathione metabolism, protein folding, and cellular response to unfolded proteins are subjects of great interest. First of all, since glutathione provides the main redox buffer against ERgenerated oxidative stress by maintaining ER oxidoreductases in a reduced state, it has been implicated in the formation of native disulfide bonds in proteins in the endoplasmic reticulum through a complex process involving not only disulfide-bond formation but also the isomerization of non-native disulfide bonds [22]. Protein disulfide isomerase (PDI) is capable of catalyzing both processes [145]. For disulfide bond formation, the ratio of reduced to oxidized glutathione ((GSH):(GSSG)) is optimal [146,147]. Disulfidebond formation in the ER can lead to the formation of ROS, and their detoxification by GSH results in an increase in GSSG [22]. Notably, increased ROS generation is linked to protein misfolding [148,149]. Decreased total glutathione levels were found to be associated with increased formation of disulfide bonds, which were simultaneously accompanied by increased formation of non-native disulfide bonds [22]. These findings could be explained by the fact that PDI's ability to isomerize non-native disulfide bonds is limited as this enzyme becomes more oxidized. Thus, it could be suggested that the elevated levels of glutathione in the ER are required for both the elimination of ROS and the facilitation of native disulfide-bond formation in folding proteins. This means that a deficiency of glutathione in the cell can disrupt these processes, causing, on the one hand, oxidative stress and, on the other hand, disruption of protein folding with subsequent activation of the UPR pathway. Notably, activation of unfolded protein response has been found to be associated with an increase in glutathione biosynthesis. It has been observed that the ATF4 and Nrf2 transcription factors link ER stress to a broader cellular response that boosts glutathione synthesis via enhancing amino acid metabolism [150,151]. Transcription factor Nrf2 plays a key role in the response to oxidative stress of the NFE2L2/NRF2 pathway [152]. Interestingly, the positive impact of NFE2L2 (gene encoding Nrf2) on the expression of the glutamate-cysteine ligase modifier subunit (GCLM) was observed only in the β-cells of individuals without type 2 diabetes, whereas no such effect was observed in diabetic samples from both datasets. This finding suggests that the NFE2L2/NRF2 pathway does not activate glutathione biosynthesis in the beta cells of patients with diabetes. The activation of PERK results in a general decrease in protein synthesis and an increase in glutathione synthesis [22]. Moreover, oxidized glutathione may enhance chaperone activity, as it has been demonstrated for the folding of alpha-crystallin [153]. Taken together, these data clearly indicate that optimal levels of glutathione are vitally important in endoplasmic reticulum for efficient protein folding and protection against oxidative stress.
Glutathione is necessary for the formation of disulfide bonds in proteins, the process in which ROS are produced as by-products, leading to impaired redox homeostasis in the ER and oxidative stress [154,155], a pathological condition prevented by glutathione. ROS have been shown to activate the UPR, which can be reversed by N-acetylcysteine by replenishing the glutathione pool. It was shown that oxidative stress may explain the activation of ER stress [154]. This means that protein folding in the ER is dependent on redox homeostasis, and oxidative stress can disrupt the mechanisms of effective protein folding, enhancing the production of misfolded or unfolded proteins and causing ER stress. Given glutathione's critical role in maintaining redox homeostasis and ensuring the process of protein folding in the endoplasmic reticulum, there is reason to believe that glutathione deficiency and the associated defective ER folding of proteins may be the cause of the increased amount of misfolded proinsulin molecules found in type 2 diabetes. Due to the fact that pancreatic beta cells have extremely high protein-synthetic activity (they produce almost 1 million insulin molecules per minute [156]), this ability makes the pancreas extremely sensitive to ER stress. Studies of type 2 diabetes mellitus show that the process of proinsulin folding is closely related to the regulation of ER stress and unfolded protein response, whose prolonged activation triggers the apoptosis of pancreatic beta cells [44].
The study has several limitations. Firstly, the range of selected genes involved in protein folding and unfolded proteins is incomplete since we have selected the main genes presented simultaneously in three metabolic databases such as KEGG, Reactome, and WikiPathways. For example, we did not look at many genes that encode enzymes from the protein disulfide isomerase gene family, which is important in enzyme-mediated disulfide bond formation in proteins. In this regard, it cannot be ruled out that alterations in the expression of genes not studied in the current study may exist. Secondly, the transcriptome data samples of pancreatic beta cells from T2D patients and non-diabetic controls used in the study were small and characterized by significant variability in gene expression across samples. Third, the comparison results obtained from the transcriptome datasets generated using different methodologies resulted in different scales of gene expression, thereby neutralizing similar trends in changes in gene expression in diabetes between the datasets. This means that establishing alterations in gene expression in type 2 diabetes relative to controls in both datasets was extremely challenging. Fourth, when interpreting the results, one should also take into account the fact that differentially expressed genes may reflect disease-induced rather than disease-causing changes in the transcriptome, as was noted by Porcu with co-workers [157]. In addition, it cannot be ruled out that changes in gene expression may be associated with other pathological processes present in patients. Furthermore, the correlations between expression levels of the studied genes should be interpreted with caution since they were inflated due to the small number of investigated pancreatic β-cell samples. Further studies on larger samples obtained at an early stage of the disease, as well as the use of unified methodological techniques for gene expression analysis, would allow a more objective assessment of the nature of gene expression changes in pancreatic beta cells in type 2 diabetes.

Datasets
Gene Expression Omnibus (GEO) datasets on gene expression profiles in pancreatic islets of Langerhans from type 2 diabetes patients and non-diabetic subjects were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/ (accessed on 14 February 2022)). The keywords 'type 2 diabetes' (All Fields) AND 'Homo sapiens' (Organism) AND 'beta cells' (All Fields) OR 'islet cells' (All Fields) OR 'pancreatic cells' (All Fields) were used to search for gene expression datasets. The list of found datasets was manually evaluated to identify datasets with transcriptome profiles of islet beta cells from T2D patients and non-diabetic controls. As a result, two datasets, GSE20966 (Dataset 1) and GSE81608 (Dataset 2), with transcriptomic profiles of pancreatic islets from patients with type 2 diabetes and those without diabetes, were identified and selected for the present study. Dataset 1 includes gene expression profiles of beta cells obtained by laser capture microdissection from cadaver pancreases of 10 subjects with T2D and 10 non-diabetics of European descent [94]. Extracted RNA samples were amplified, biotinylated, and hybridized to the GeneChip Human X3P Array (Affymetrix, Santa Clara, CA, USA). Dataset 2 represents the transcriptomes of 1492 single human pancreatic islet cells such as α-, β-, δand pancreatic polypeptide (PP cells) that were obtained from cadaveric pancreas of 6 T2D and 12 non-diabetics (70% were Europeans, 30% were Asians or Hispanics) [158]. Pancreatic islet cell types were hybridized with mRNA probes for human GCG, INS, SST, and PPY genes to estimate bimodal expression distribution of the genes. Isolated RNA samples from each cell type were analyzed using single-cell RNA sequencing on Illumina HiSeq2500 (Illumina, San Diego, CA, USA). RNA-Seq data analysis was performed after the quality control of gene expression. Expressed genes were defined by ≥ 1 RPKM. The transcriptome data of single cells (β-cells) obtained from each donor were averaged.

Selection of Genes
The following groups of genes or pathways were subjects of interest for our study: (1) genes encoding molecular chaperones or those implementing protein folding; (2) genes involved in the unfolded protein response (UPR); and (3) genes encoding enzymes regulating glutathione metabolism. The three pathway databases KEGG (https://www.genome. jp/kegg (accessed on 24 April 2022)), Reactome (https://reactome.org (accessed on 27 April 2022)), and WikiPathways (https://www.wikipathways.org (accessed on 24 April 2022)) were used to select genes representing each of the above groups. In particular, pathways "Glutathione metabolism" (KEGG, WikiPathways), "Cysteine and methionine metabolism" (KEGG), "Glutathione synthesis and recycling" (Reactome), "Gamma-glutamyl cycle for the biosynthesis and degradation of glutathione, including diseases" (WikiPathways) were used to select genes regulating glutathione metabolism. "Protein processing in endoplasmic reticulum "(KEGG) and "Unfolded Protein Response" (Reactome, WikiPathways) pathways were used for selection of genes encoding molecular chaperones or involved in protein folding and unfolded protein response. The criterion for gene selection was the representation of a gene in at least two of the three pathway databases. The CTH gene was added to the list of glutathione metabolism genes because it encodes cystathionine gamma-lyase, which catalyzes the last step in the trans-sulfuration pathway from L-methionine to L-cysteine, which is utilized by cells for biosynthesis of glutathione (https://www.uniprot.org/, accessed on 3 July 2023). The DNAJC3 [45], EDEM2 [16], and GFER [89] genes were also selected for this study because they were found to be associated with T2D susceptibility in earlier studies. Moreover, the insulin gene (INS) was included in the gene list to investigate the correlation between its expression and the studied genes. As a result, raw expression data on 142 genes were selected from whole transcriptomic datasets 1 and 2, with 86 genes encoding molecular chaperones or those involved in protein folding (PF), 36 genes representing the three arms of the UPR pathway, and 20 genes encoding glutathione metabolism enzymes (GME). The full list of genes selected for the study is presented in Supplementary Table S1.

Data Analysis
The limma R software package version 4.3 was used to investigate differentially expressed genes (DEGs) in pancreatic beta cells of T2D patients and non-diabetic controls from Dataset 1 (GSE20966) [159]. Initially, outlier values in the gene expression data were removed from the entire dataset. Expression data were log2-transformed and normalized. Dataset 2 (GSE81608) was processed with the same approach using GREIN (http://www.ilincs.org/apps/grein/, accessed on 6 June 2023), an interactive web platform for exploring and analyzing GEO RNA-seq data to obtain differentially expressed genes [160]. The log fold-change (LogFC) measure, the log-ratio of a gene's expression values in two analyzed samples, was used to describe the quantity of changes in gene expression existing between T2D patients and non-diabetic controls. The false discovery rate (FDR) to control for multiple comparisons (FDR < 0.05) was calculated by FDR online calculator (https://www.sdmproject.com/utilities/?show=FDR, accessed on 6 June 2023) to identify differentially expressed genes among the 143 selected genes. The expression levels of 143 genes were extracted from both the transcriptome datasets for further statistical analysis using STATISTICA Ultimate version 13.0 (TIBCO, Palo Alto, CA, USA). Pearson's correlation coefficients were calculated to measure the relationships between expression levels of the genes separately in diabetics and non-diabetics for each dataset. Hierarchical cluster analysis of correlation matrices (Ward's method) was applied to visualize the groups of co-expressed genes. Linear regression analysis was used to evaluate whether the expression of genes for glutathione metabolism has an impact on the expression levels of genes regulating protein folding and unfolded protein response in pancreatic beta cells of patients with type 2 diabetes and non-diabetic individuals.

Conclusions
The accumulation of misfolded and unfolded proteins in the endoplasmic reticulum is known to cause ER stress by activating a signaling pathway known as the unfolded protein response [161,162], focusing on restoring homeostasis by inhibiting the transcription of genes, enhancing the folding capacity of the ER, and degrading unfolded proteins [163]. However, chronic activation of UPR makes ER fail to recover protein homeostasis, resulting in cellular dysfunction and apoptosis [164]. Numerous studies have demonstrated that chronic ER stress and associated UPR are characteristic features of type 2 diabetes, and accumulation of misfolded proinsulin molecules in β-cells is detected early in prediabetes and subsequently exacerbated, leading to apoptosis [19][20][21]. The primary mechanisms underlying ER stress activation and subsequent induction of beta cell apoptosis in type 2 diabetes mellitus remain unknown.
The present study is the first to discover differentially expressed genes encoding enzymes for glutathione metabolism, molecular chaperones, and the UPR pathway in pancreatic beta cells of patients with type 2 diabetes. Although differentially expressed genes varied across the studied transcriptome datasets, the common findings were a decrease in the expression levels of two genes for glutathione metabolism, GCLC and PGD, as well as differentially expressed genes regulating protein folding, such as POM121C and RPS19BP1. Three other DEGs regulating glutathione metabolism, ANPEP (Dataset 1) and IDH2 and CTH (Dataset 2), were observed in pancreatic beta cells of diabetic patients. Furthermore, beta cells of T2D patients showed a dataset-specific gene signature for DEGs involved in the regulation of protein folding. In particular, beta cells of diabetics exhibited decreased expression of genes encoding molecular chaperones belonging to the heat shock families Hsp70 (HSPA1B, HSPA7, and HSPA8) and Hsp90 (HSP90AA1 and HSP90AB1) as well as co-chaperones (BAG3 and ST13) and nucleoporins (NDC1, NUP160, and POM121C). T2D samples also exhibit DEGs regulating the unfolded protein response: decreased expression of BID, a protein that induces caspase and apoptosis pathways (Dataset 2), and increased expression of CREB3L4 and ERP27 (Dataset 1). The expression of molecular chaperones and genes involved in the UPR pathway in T2D patients was influenced by the expression of genes responsible for de novo biosynthesis of glutathione, GCLC, GCLM, and GSS, which supports the hypothesis that dysregulation of glutathione biosynthesis might be a key condition responsible for the impaired folding of proinsulin, triggering an unfolded protein response, and ultimately leading to β-cell apoptosis and type 2 diabetes [16]. Undoubtedly, these findings warrant further investigation. The discovery of gene networks involved in the coordinated regulation of gene expression modulating glutathione metabolism, protein folding, and unfolded protein response will allow for a better understanding of the molecular mechanisms governing redox homeostasis in type 2 diabetes, as well as unravel the cross-talk between impaired glutathione metabolism and proteostasis, linking them to well-recognized pathological processes underlying disease development and progression.