{'fraction_non_alphanumeric': 0.03480013436345314, 'fraction_numerical': 0.02418542156533423, 'mean_word_length': 4.9496402877697845, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 5, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bioinformatics Analysis Identifies Potential Ferroptosis Key Gene in Type 2 Diabetic Islet Dysfunction

Background Islet β cells dysfunction (IBCD) is a cortical component in pathogenesis of type 2 diabetic mellitus (T2DM). However, the relationship of ferroptosis and IBCD remains unknown. This study was aimed to screen potential ferroptosis key genes to reveal latent physiological and pathological process of IBCD in T2DM. Methods Firstly, T2DM key genes were screened by combining with differentially expressed genes (DEGs) analysis and WGCNA. Then, ferroptosis-related genes (FRGs) in IBCD of T2DM were identified by taking the intersection between T2DM key genes and FRGs. Finally, T2DM-FRGs were validated in another T2DM dataset as well as islet single-cell RNA sequencing dataset and the miRNA regulated T2DM-FRG was predicted by using four miRNA databases. Results 89 T2DM key genes were identified between DEGs and WGCNA. Then, 3 T2DM-FRGs were screened by taking the intersection of T2DM key genes and FRGs, namely ITGA6, MGST1 and ENO2. At last, MGST1 were validated as the T2DM-FRG in another T2DM islet issues dataset and islet single-cell RNA sequencing dataset. Conclusion MGST1 may be the potential ferroptosis key gene of IBCD in T2DM.


INTRODUCTION
According to International Diabetes Federation (IDF) reports, there are approximately 536.6 million adult diabetic peoples worldwide and almost half of patients were unaware of their diabetes status (1,2). Some 90% of diabetic individuals have type 2 diabetes mellitus (T2DM), which is currently the most common type of diabetes mellitus and increases the burden on the state. These patients have higher risks of cardiovascular damage, nerve damage, bacterial infection and other complications (3). The typical features of T2DM are glucotoxicity, lipotoxicity, insulin resistance and islet b cells dysfunction (IBCD) (4). IBCD occurs in progressive stage of T2DM and influences the prognosis of T2DM. Understanding the clear mechanism of IBCD can help us to early identify IBCD and perform targeted therapy in T2DM. At present, studies show that the mechanisms of IBCD in T2DM include mitochondrial dysfunction, inflammation, autophagy, oxidative stress (5). But the specific pathogenesis of IBCD in T2DM is still not completely clear and failure to discerned early.
Ferroptosis, a form of nonapoptotic regulated cell death (RCD), is caused by the lipid peroxide accumulation dependent on iron and antioxidant system dysfunction (6). It is worth noting that islet b cells express low level of antioxidant enzymes such as superoxide dismutase (SOD), glutathione (GSH) peroxidase and catalase, so they are susceptible to suffer from oxidative stress (7). Interestingly, lipid peroxide and iron accumulation were observed in type 2 diabetic islet b cells (8,9). In addition, the islet b cells significantly died and its capacity of glucose-stimulated insulin secretion reduced, when they were treated with ferroptosis inducer erastin in vitro (10). These adverse effects improved in ferroptosis inhibitors Ferrostatin-1-and desferrioxamine-treated islets b cells (10). Thus, IBCD in T2DM probably related to ferroptosis. Probing the underlying molecular mechanism between IBCD in T2DM and ferroptosis may be a breakthrough for its early diagnosis and treatment.
Here, based on the above researches, we screened the key genes by differential expression genes (DEGs) analysis and weighted gene co-expression network analysis (WGCNA) (11) in order to explore the effect of ferroptosis on type 2 diabetic IBCD. These T2DM key genes were intersected with the ferroptosis-related genes (FRGs) to acquire T2DM-FRGs in IBCD. We also validated the T2DM-FRGs using another T2DM islet issues dataset as well as islet single-cell RNA sequencing dataset and predicted microRNAs (miRNAs) regulated T2DM-FRG.

Differentially Expressed Genes Analysis
DEGs were identified by "limma" R package in R software. We considered DEGs as p-value <0.05 and |log 2 fold change (FC)| >log 2 1.2. DEGs were visualized by the "volcano" and "heatmap" R packages. The common up-regulated and downregulated DEGs between GSE38642, GSE41762 and GSE76894 were identified and visualized by Venn diagram in online tool Hiplot (https://hiplot.com.cn/).

Functional Enrichment Analysis
Online tool DAVID (https://david.ncifcrf.gov/) (19) was used to perform the gene ontology (GO) (20) and kyoto encyclopedia of genes and genomes (KEGG) analyses (21). p-value <0.05 was considered as the threshold. GO analysis including biological process (BP) and molecular function (MF). The records of GO and KEGG analysis were visualized as chord chart and bubble diagrams by "ggplot2" R package respectively.

Weighted Gene Co-Expression Network Analysis (WGCNA)
WGCNA is a well-established method for analyzing the relationship between diseases phenotype and genes (11). The dataset GSE76894, which had a large sample size, was chosen to find key genes related type 2 diabetic phenotype by "WGCNA" R package. A scale-free network was constructed by calculating different power value (range from 1 to 30), and the suitable power value was chosen when scale free R 2 was 0.85. Subsequently, we performed the cluster analysis to recognize the highly similar modules and the min module size set as 30. The connection between modules and traits were calculated to identify the key trait-module genes. In the module-trait analysis, the module key genes were chosen with module membership (MM) >0.8 and gene significance (GS) >0.2.

Statistical Analysis
GraphPad Prism 8.0 was used to perform statistical analysis and drawing. We selected unpair t-test as the statistical method. pvalue <0.05 was considered as statistically significant. Data were presented as mean ± standard deviation (SD).

Identification of DEGs Between T2DM and Normal Islet Issues
Three datasets GSE38642, GSE41762 and GSE76894 were analyzed by "limma" package in R with p <0.05 and  Figure 1A). In dataset GSE41762, we extracted 712 DEGs, which contains 324 upregulated genes and 388 down-regulated genes ( Figure 1C). 2220 DEGs, comprised of 937 up-regulated genes and 1283 down-regulated genes were identified in dataset GSE76894 ( Figure 1E). The top 200 DEGs in these datasets were respectively performed clustering and shown in the form of heatmaps ( Figures 1B, D, F). 44 up-regulated DEGs and 160 down-regulated DEGs from intersection of three datasets was chosen by online tool Hiplot (https://hiplot.com.cn/) and shown in Venn diagram (Figures 2A, B). The percentages and gene names of DEGs in each dataset were shown in Supplementary Material Table 1.

Results of Functional Enrichment Analysis
DEGs between T2DM and Normal islet tissues were mainly enriched in the following BP terms based on enrichment analysis: gene expression, apoptosis, insulin secretion, cellcell adhesion, gluconeogenesis, ion transmembrane transport, exocytosis, cellular response to hypoxia. MF of DEGs were primarily enriched for transmembrane transporter activity, magnesium ion binding, integrin binding, SNARE binding, calcium channel activity ( Figure 2C). KEGG analysis shown that DEGs were mainly participated in PI3K-Akt signaling pathway, MAPK signaling pathway and insulin secretion ( Figure 2D). The concrete GO and KEGG terms of DEGs in each dataset were shown in Supplementary Material Table 2.

Construction of Co-Expression Module and Identification of Key Trait-Module Genes
WGCNA, a good method to find the relationship of genes and traits, was performed to excavate key genes significantly related to traits by combining with DEGs. Sample clustering tree and phenotypic heatmap were shown in Figure 3A. The darker the color, the higher the expression in phenotypic heatmap, and the grey represented missing values. No sample was removed in dataset GSE76894 for no outlier. The top 75% genes of the median absolute deviation (MAD) in this dataset were screened for WGCNA. When the soft-thresholding power was equal to 6, the evaluation parameter R 2 was greater than 0.85 ( Figures 3B, C). A scale-free network was constructed with the soft-thresholding power of b = 6 and the genes were classified into different modules according to expression profiles. We identified a total of 24 modules via the average linkage clustering ( Figure 3D). Grey indicated genes that had not been classified into modules. Green module and turquoise module had the most significant gene significance (GS) in T2DM and there was no conspicuous module in other traits. ( Figure 3E). Finally, we respectively extracted 73 genes in green module and 397 genes in turquoise module with MM >0.8 and GS >0.2 (Figures 3F, G, Supplementary Material Table 3). These 470 genes were considered as the trait-module key genes.
Screening and Validation of T2DM-FRGs 88 genes were elected from the union set between DEGs and trait-module key genes in WGCNA, which was visualized in Venn diagram ( Figure 4A). These genes were intersected with ferroptosis-related genes contained 599 genes subsequently, and we obtained 3 T2DM-FRGs (ITGA6, MGST1 and ENO2) eventually, shown in Figure 4B. MGST1 and ITGA6 increased in T2DM group, while ENO2 decreased. GSE76895 and E-MTAB-5061 datasets were used to verified T2DM-FRGs. In single-cell sequencing dataset E-MTAB-5061, we selected the b cell samples to perform validation. The differential expression of MGST1, not ITGA6 and ENO2, was proved to be statistically significant in two validation datasets ( Figures 4C-H).

Prediction of miRNA That Regulated Hub Gene
As shown in Figure 4I, has-miR-3619-5p may regulated MGST1 by prediction in Targetscan

DISCUSSION
IBCD is considered as the important factor in the progression of T2DM and is also another direction to treat T2DM. Longterm hyperglycemia, hyperlipidemia and uncontrol of blood glucose in patients with T2DM caused the islet b cells damage and even death (4). Some ways of islet b cells death in T2DM, such as apoptosis and autophagy were confirmed in studies (22,23). But the early diagnosis and treatment of IBCD in T2DM is still a difficult and unsolved problem. As a newly identified iron-dependent form of oxidative cell death, ferroptosis remains unclear in T2DM. The features that type 2 diabetic islet b cells appear low antioxidant activity and iron overload indicated us the possible link between IBCD and ferroptosis (7,9). We hope to discover the association between the two and detect the target for early diagnosis and treatment in IBCD of T2DM by our study. MAPK and PI3K-Akt signaling pathways, which have been proved as the pathways regulated ferroptosis, were significantly enriched in KEGG of DEGs in IBCD of T2DM (24,25). However, no significant ferroptosis-related terms were found in GO of DEGs. We considered that serious ferroptosis may have not occurred in IBCD at pre-middle stage of T2DM, but changes in ferroptosis-related pathways, such as MAPK and PI3K pathways. The differential genes between T2DM and Normal islet tissues were screened by combining DEGs and trait-related genes in WGCNA, which ensured the reliability of data. We obtained 3 T2DM-FRGs (ITGA6, MGST1 and ENO2) from the intersection between differential genes and ferroptosis-related genes. We found MGST1 and ITGA6 is up-regulated, while ENO2 is downregulated in T2DM. MGST1, not ITGA6 and ENO2, was considered as the key T2DM-FRG after validation in another type 2 diabetic islet issues dataset and single-cell sequencing dataset. This study mays provide reference for research the connection between IBCD and ferroptosis in T2DM.
MGST1 (Microsomal Glutathione S-transferase 1) is a membrane-bound transferase regulates oxidative stress. It mainly located in mitochondria, ER plasma membrane and peroxisome (26). Researchers found that it fighted ferroptosis in pancreatic cancer cells. MGST1 is regulated by Nrf2 (Nuclear Factor erythroid 2-Related Factor 2) and partly combined with ALOX5 (Arachidonate Lipoxygenases 5) to inhibit ferroptosis (27). In addition, MGST1 is up-regulated during ferroptosis (27). We prove that mRNA expression of MGST1 increases in type 2 diabetic IBCD. Thus, it is conceivable to hypothesis that MGST1 mays contribute to alleviate ferroptosis in type 2 diabetic IBCD, which needs to be proved experimentally in the future. ENO2 (Enolase 2 or gamma enolase, also known as neuron-specific enolase, NSE), is a glycolytic enzyme which catalyzes 2-phosphoglycerate to form PEP (Phosphoenolpyruvate). It not only participates in glycolysis and gluconeogenesis, but also is related with inflammation and ferroptosis (28,29). ENO2 mays influence ferroptosis by neuroglobin (29). A study indicates that ENO2 is positively correlated with insulin secretion and negatively correlated with HbA1c (30). ENO2 is decreased in cultured INS-1 cells treated with high glucose in 12 hours and in type 2 diabetic islet cells (30). These results in vivo and in vitro suggest that ENO2 is a key molecule protect the function of islet b cells. However, mRNA expression of ENO2 did not change in islet single-cell RNA sequencing samples. This may be related to the small single-cell sequencing sample size. Whether ENO2 regulates islet b cells in T2DM through ferroptosis pathway remains to be studied.
In addition, the study about islet issue sequencing has certain limitations. There is no way to know whether the target gene changes are present in b cells because islet issues include a variety of cells, such as a cells, b cells and PP cells.
There are many studies about islet single-cell RNA sequencing of T2DM at present (16,31,32). Therefore, we used the islet  single-cell sequencing dataset to check the target genes again. Single-cell sequencing, though, is more specific to get a lot of b cells and its sequencing data, but contains a few samples in patients for its high cost, which will cause the inability to exclude individual differences. Finally, we selected MGST1 as the potential ferroptosis-related gene in IBCD of T2DM by combing large patient samples of islet issue sequencing datasets and specificity of islet single-cell sequencing dataset. This study reveals the potential relationship between ferroptosis and IBCD in T2DM based on bioinformatics analysis. It deeply adds to our understanding of type 2 diabetic islet b cells impairment and provides reference for developing targeted therapy for T2DM.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
HY designed the study and wrote manuscript, HY and RW conducted the bioinformatic analysis. JW and YW took part in a discussion. XZ and LW amended the manuscript and co-wrote the final draft with HY. All authors contributed to the article and approved the submitted version.