Integrated Bioinformatics Analysis Revealed Immune Checkpoint Genes Relevant to Type 2 Diabetes

Introduction

Diabetes has become one of the most serious public health problems in the world with its increasing morbidity and mortality. There are about 537 million adults worldwide with diabetes, and approximately 90% of them have T2D.¹ Increasing evidence showed that chronic low-grade inflammation and innate immune system activation are contributed to the pathogenesis of T2D² and some anti-inflammatory agents may have effects on T2D, such as aspirin³ and interleukin-1 receptor antagonist.⁴ However, the specific mechanisms about the activation and regulation of this inflammatory state in T2D remain unclear. Therefore, figuring out how the pancreatic islet inflammation modulates the progression of T2D may provide a new insight into the pathogenesis of T2D.

Immune checkpoints, most of which are inhibitory pathways in the immune system, are critical to keep self-tolerance and regulate type, intensity, and duration of immune response.⁵ Immune checkpoint blockade (ICB) therapy with the use of inhibitors, which target different immune checkpoints, has emerged as a promising therapy in the treatment of several types of cancer, such as gastric cancer,⁶ hepatocellular carcinoma,⁷ and pancreatic cancer.⁸ As potential treatment targets for ICB therapy, immune checkpoint genes (ICGs) are involved in immune pathways, and the expression levels of ICGs are important biomarkers for ICB therapy;⁹ however, it is not well investigated in the field of T2D. Several studies reported that immune checkpoint inhibitors induced diabetes or aggravated pre-existing T2D.^10,11 Also, a recent study showed that CD24 -Siglec-E interaction is a protective immune checkpoint against metabolic dysfunction and metaflammation,¹² suggesting that immune checkpoint molecules may be involved in the pathogenesis of T2D. Meanwhile, Sun et al found that PD-1 expression levels on immune cells were significantly lower in patients with T2D but were positively associated with insulin levels and duration of diabetes.¹³ However, the roles of immune checkpoints in T2D especially in the progress of pancreatic immune cell infiltration and β cell functions has not yet been addressed.

In this study, we screened out the differential expressions of ICGs between T2D patients and normal people in pancreatic islets from the transcriptome sequencing gene data of GEO. A risk model based on these ICGs was built. The key pathways and proteins were identified from the analysis of GOs, KEGG pathways, and GSEA. The hub genes were obtained via combined WGCNA and PPI networks. Further, we assessed immune cell infiltration in T2D by CIBERSORT. This study offered a new clue to reveal the pathogenesis of T2D and potential therapeutic targets for T2D.

Materials and Methods

Microarray Data Acquisition

The gene matrix was obtained from the GEO database. We conducted screening was according to the following criteria: (1) the gene sets must include experimental and control groups, (2) T2D patients, (3) the samples of datasets more than 10, and (4) the type of samples are pancreatic islets. Finally, GPL570 dataset GSE76894 including 19 cases of T2D and 84 control cases was downloaded from the GEO datasets for main analysis. GPL16791 dataset GSE164416 including 39 cases of T2D and 18 cases of control was selected as validation dataset. GSE76894 contained the islet transcription of islets isolated by enzymatic digestion from 103 organ donors.¹⁴ GSE164416 contained samples collected from 133 metabolically phenotyped pancreatectomized patients.¹⁵

ICGs Selection and Analysis

To analyze the expression of ICGs, we screened 92 ICGs after an extensive literature review, as Hu et al did.¹⁶ Then, after intersection with the expression profiles of GSE76894, 42 genes were left: BTN2A1. BTN2A2, BTN3A1, CD209, CD274, CD276, CD40, CD44, CD47, CEACAM1, HAVCR2, HLA-A, HLA-B, HLA-C, HHLA2, HLA-DMA, HLA-DMB, HLA-DOA, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQB1, HLA-DRA, HLA-E, HLA-F, HLA-G, ICOSLG, ICOSLG.1, NRP1, LAIR1, LGALS9, PVR, SIRPA, TDO2, TMIGD2, TNFRSF14, TNFRSF18, TNFRSF4, TNFSF15, TNFSF4, TNFSF9, VTCN1.

We normalized the GSE76894 using the robust multiarray averaging algorithm. To analyze the effects of expression of ICGs on T2D, the R package limma was used to screen the DEGs in GSE76894.¹⁷ The thresholds were |log2 (fold change)| (|log2FC|) ≥0.5 and adjusted P value <0.05. The genes with log2FC ≥0.5 and adjusted. P Values <0.05 were upregulated, and the genes with log2FC ≥0.5 and adjusted.P.Value <0.05 were downregulated. The downregulated or upregulated DEGs were shown in volcano plots. To analyze the expressions of ICGs between the T2D and control groups, the box plots of the two groups were plotted with the two groups compared using Wilcoxon rank sum test. The significantly different genes between the T2D and control groups were involved in the following analysis.

To explore the correlations between ICGs in all patients, the Pearson correlations were calculated. The absolute values of correlation coefficient larger than 0.3 and P<0.05 indicated the presence of correlation.¹⁸

To explore the interactions between the chosen ICGs-related proteins, we used the online STRING to build the PPI network of the chosen ICGs at the confidence level of 0.4.¹⁹ Then, the PPI network was further analyzed and visualized on Cytoscape 3.9.1.²⁰

Risk Model Based on ICGs

T2D and control samples may express different levels of ICGs, which may affect T2D. We built a risk model based on ICGs. We imported the chosen ICGs into ridge regression on the R package glmnet,²¹ and got the optimal λ and the coefficients of the chosen ICGs. The coefficients of the chosen ICGs were multiplied by the expression level of each gene and added together to get a T2D risk score. Then, the T2D risk scores of the two groups were compared using Wilcoxon rank sum test and presented as dot plots. Then, we plotted receiver’s operating characteristic curve (ROC) and calculated the area under curve (AUC) using IBM SPSS Statistics, Version 26 (IBM Corporation, Armonk, New York, USA).

Then, the risk model was employed in GPL570 dataset GSE156993 including PBMC (peripheral blood mononuclear cell) samples of 12 T2D cases and 6 control cases to validate its clinical application potential.

Functional Enrichment Analysis

To further explore the different biological characteristics between the T2D samples and the control samples, we used GO annotations to characterize biological properties of DEGs, including molecular function, cellular component, and biological process.²² KEGG pathway enrichment analysis was conducted to determine functional attributes of DEGs.²³ To study the different biological states between the T2D samples and the control samples, we subjected all DEGs to GSEA. We downloaded reference gene sets “c5.all.v2022.1.Hs.entrez.gmt” from database MSigDB.²⁴

Identification of Gene Clusters and Hub Genes

WGCNA can identify clusters of correlated genes and find key genes by connecting gene network to clinical traits. All DEGs were chosen as input dataset for WGCNA using the R package WGCNA.²⁵

To explore the gene associations of the T2D and control group, the protein–protein interaction (PPI) network was built according to the selected module of WGCNA by the online STRING at the confidence level of 0.4.¹⁹ We then imported the result into Cytoscape 3.9.1 for further analysis and visualization.²⁰

We calculated the network attributes of nodes and identified core genes in the network by cytoHubba.²⁶ The hub genes were identified as the top 10 nodes ranked by the betweenness centrality. These genes were highly correlated with others and may play major roles in biological processes, which are important in the progression of T2D. To further explore the correlation of ICGs and hub genes, we calculated the Pearson correlations between them. The absolute values of correlation coefficient larger than 0.3 and P<0.05 indicated the presence of correlation.

Analysis of Immune Infiltrating Cells

Beta cells of pancreatic islets are in a microenvironment consisting of other endocrine cells, vascular endothelial cells, extracellular matrix, and immune cells. Immune microenvironment composed of immune molecules, immune cells and cytokines is an important part of islet microenvironment. Immune checkpoint molecules are critical in modulating tissue microenvironment. Therefore, it is necessary to analyze the immune infiltrating cells in islets. CIBERSORT is an algorithm for analyzing the immune infiltration and estimating the abundance of 22 types of immune cells.²⁷ The CIBERSORT of R language²⁷ was used to calculate the abundance of immune cells between the T2D group and the control group. Wilcoxon rank sum test was used to compare the proportion of islet immune cells in the normal and T2D groups at the significant level of P <0.05, and the result was visualized as a box plot.

Furthermore, to investigate the correlations of the immune infiltrating cells or the correlations of selected ICGs and the hub genes with the immune infiltrating cells, we calculated the Pearson correlations between them. The absolute values of correlation coefficient larger than 0.3 and P<0.05 indicated the presence of correlation.

The Investigation of Disease Gene Expression Levels

T2D-related disease genes, insulin secretion- related genes, and β cell function-related genes were downloaded from the GeneCards database (https://www.genecards.org/). We analyzed the level of genes and the T2D and control groups were compared using Wilcoxon rank sum test. The significant level was P < 0.05. The correlations between ICGs, T2D-related disease genes, insulin secretion-related genes or β cell function-related genes were analyzed via Pearson correlation. The absolute values of correlation coefficient larger than 0.3 and P<0.05 indicated the presence of correlation.

Druggability Evaluation

We searched the interaction between the identified ICGs and drugs using DGIdb,²⁸ ChEMBL,²⁹ and DrugBank³⁰ databases to further assess whether these ICGs can be potential therapeutic targets. We documented information on names, the development process, and the indication of drugs.

RT-PCR Validation of the ICGs

To confirm the findings from the bioinformatics analysis, pancreatic islet tissue from 5 male leptin receptor-deficient db/db mice mimicking the conditions of T2D and 5 wild-type mice were harvested for RT-PCR validation. Total RNA from islet tissue was extracted with RNeasy Lipid Tissue Mini Kit (QIAGEN GmbH, Germany). RNA samples from total RNA were reverse transcribed to cDNA using High Capacity RNA-to-cDNA Kit (Thermo Fisher Scientific, USA), and RT-PCR was carried out using Power Up SYBR Green Master Mix (Thermo Fisher Scientific, USA). Actin was used as an internal reference. Relative mRNA expression was calculated using the 2^−ΔΔCt method. The sequences of primer used here are provided in Table S1.

Statistical Analysis

Data were analyzed and processed on R 4.2.1 or IBM SPSS Statistics, Version 26. Continuous variables in normal distribution were compared between groups using independent Student t test, those with equal variances not assumed were compared using Welch’s t test, and those in non-normal distribution were compared using Wilcoxon rank sum test. The correlation coefficients between genes were calculated via Pearson correlation analysis. All P values are two-sided. P < 0.05 indicates significance.

Results

ICGs Selection and Analysis

The procedure of this study was conducted based on the flow diagram (Figure 1). To analyze the difference between the T2D and control group, we found 458 DEGs, including 120 upregulated genes and 338 downregulated genes (Figure 2A). A heat map of DEGs was presented using gene clustering to differentiate the two groups (Figure 2B). Histogram of ICGs showed CD44, CD47, HAVCR2, SIRPA, TNFSF9, and VTCN1 were significantly different between the two groups (Figure 2C). The correlation analysis of the selected ICGs showed that CD44, CD47, SIRPA, TNFSF9, and VTCN1 were highly associated (Figure 2D). A PPI network of the selected ICGs showed that these 6 ICGs were interactive (Figure 2E).

Figure 1 Flow diagram presenting the main plan and process of the study.

Figure 2 Selection and analysis of ICGs and DEGs. (A): Volcano plots of DEGs, x-axis: log2FoldChange, y-axis: -log10 (adjust P-value); red, gray and blue nodes indicate the differentially expressed genes are upregulated, insignificant, and downregulated, respectively. (B): Heat map of DEGs, blue: control group, red: type 2 diabetes group. (C): Histogram of expressions of ICGs in type 2 diabetes group and control group; x-axis: ICGs, y-axis: gene expression level; red: type 2 diabetes group, blue: control group. *: P <0.05, **: P <0.01. (D): Correlations among ICGs. Colors indicate correlations. A redder color means stronger positive correlation while a bluer color means stronger negative correlation. *: P <0.05. (E): PPI network of ICGs.

Risk Model Based on ICGs

To assess the effect of the ICGs on T2D, we used ridge regression and obtained the coefficients of these 6 ICGs (Figure 3A and B). Then, we calculated a T2D risk score and plotted a dot plot. Results showed that the T2D risk score of the case group was significantly higher than the control group (Figure 3C). Also, we plotted an ROC curve and calculated the AUC (Figure 3D).

Figure 3 Diagnostic model based on ICGs. (A-B): ICGs identified by ridge regression. (C): Type 2 diabetes risk score of case group and control group. Red: type 2 diabetes group, blue: control group. ***: P <0.001. (D): ROC curve of type 2 diabetes risk model. AUC: area under curve; CI: confidence interval.

Functional Enrichment Analysis

The GO analysis showed DEGs were mainly enriched in bioprocess of signal release, peptide secretion and transport, insulin secretion, and hormone secretion, cellular components of presynapse, transport vesicle, neuronal cell body, exocytic vesicle, and synaptic vesicle, and molecular function of metal ion transmembrane transport activity, SNARE binding, cation channel activity, and ion channel activity (Figure 4A, Table S2). DEGs were enriched in KEGG pathways of insulin secretion, Ras signaling pathway, cytokine–cytokine receptor interaction, PI3K-Akt signaling pathway, and NOD-like receptor signaling pathway (Figure 4B, Table S3).

Figure 4 Functional enrichment analysis. (A): GO enrichment results, x-axis: gene count, y-axis: GO terms. (B): KEGG enrichment results. (C): GSEA results.

Then, all genes between groups were subjected to GSEA (Table S4). Results showed that the following biological processes were significantly enriched. In the T2D group, bioprocess including response to cytokine, regulation of response to stress, cytokine mediated signaling pathway; biological process involved in interspecies interaction between organisms, defense response to other organism, abnormal joint morphology, defense response, innate immune response; anatomical structure formation involved in morphogenesis, tube morphogenesis, positive regulation of gene expression, and epithelium development were activated, while positive regulation of cytosolic calcium ion concentration was inhibited (Figure 4C, Table S3).

Identification of Gene Clusters and Hub Genes

To explore the associations between DEGs, we imported the DEGs to WGCNA (Figure 5A). One coexpression gene module was identified (Figure 5B). Then, only a set of 457 key genes was selected and involved in subsequent analysis (Figure 5B-D).

Figure 5 WGCNA. (A): WGCNA of threshold screening. (B): Coexpression gene clustering. (C): Correlations between gene clusters and type 2 diabetes patients; (D): Correlation analysis between most significant gene clusters and type 2 diabetes.

Visualizing through Cytoscape, the network contains 1093 edges and 324 nodes (Figure 6A). The color shades of the node indicated its level of betweenness centrality (BC). Specifically, LI1B is closely connected with 20 DEGs. The function interactive subnetworks including 10 genes were extracted on CytoHubba (Figure 6B). We then plotted an ROC curve of the 10 hub genes and the results showed that these genes can well differentiate the two groups of samples (Figure 6C).

Figure 6 PPI network of DEGs. (A): PPI network of key genes; (B): Subnetwork of PPI network. (C): ROC curve of hub genes. (D): Correlation analysis between ICGs and hub genes. Colors indicate correlations. A redder color means stronger positive correlation while a bluer color means stronger negative correlation. *: P <0.05.

We further calculated the correlations of ICGs and hub genes. Results showed that CD44, CD47, SIRPA, TNFSF9, and VTCN1 were closely associated with the expression of several hub genes (Figure 6D).

Analysis of Immune Infiltrating Cells

CIBESORT (Figure 7A) showed the concentrations of eosinophils, T cells CD4 naive, and T cells regulatory were significantly higher in T2D patients. Compared with the control group, the level of T cells CD4 memory resting and monocytes were lower in the patients (Figure 7A). We also calculated the correlations of the selected ICGs and hub genes with immune cell contents. We found that the expressions of several ICGs were associated with B cells naïve, plasma cells, T cells CD4 resting, T cells regulatory, and macrophages M2 (Figure 7B). T cells CD4 memory resting and dendritic cell activated were positively correlated with several hub genes, and T cells regulatory were negatively correlated with some hub genes (Figure 7C). Also, we calculated the correlations of immune cell contents between the control group and T2D group. In the T2D patients, the B cells naïve, plasma cells, T cells CD8, NK cells, and neutrophils were significantly correlated with several types of immune cells (Figure 7D). In the control group, the content of T cells CD4 naive was significantly correlated with other types of immune cells (Figure 7E).

Figure 7 Analysis of Immune Infiltration. (A): Histograms of immune cell concentrations; x-axis: immune cells; y-axis: cell concentration; red: type 2 diabetes samples; blue: control samples. (B): Correlations between immune cells and ICGs. Colors indicate correlations. A redder color means stronger positive correlation while a bluer color means stronger negative correlation. *: P <0.05. (C): Correlations between immune cells and hub genes. (D–E): Correlations of immune cell concentrations in the test group (D) and control group (E); red: positive correlation, blue: negative correlation. A redder color means stronger positive correlation while a bluer color means stronger negative correlation. *: P <0.05.

The Investigation of Disease Gene Expression Level

A total of 13,237 T2D-related disease genes, 3184 insulin secretion-related genes, and 1666 β cell function-related genes were obtained from the GeneCards database. The levels of genes between T2D and control group were analyzed (Figure 8A). The expression of ICGs was significantly related to the expression of several T2D-related genes, insulin secretion-related genes, and β cell function-related genes (Figure 8B-D).

Figure 8 Correlation between the level of ICGs and several T2D-related genes. *: P <0.05, **: P <0.01, ***: P <0.001. (A): Box plots displaying the expression of the top 20 genes related to T2D. (B): Correlation analysis between ICGs and top 20 type 2 diabetes-related genes. A redder color means stronger positive correlation while a bluer color means stronger negative correlation. *: P <0.05. (C): Correlation analysis between ICGs and top 20 insulin secretion-related genes. A redder color means stronger positive correlation while a bluer color means stronger negative correlation. *: P <0.05. (D): Correlation analysis between ICGs and top 20 β cell function-related genes. A redder color means stronger positive correlation while a bluer color means stronger negative correlation. *: P <0.05.

Druggability Evaluation

In druggability evaluation, we found that CD44, CD47, SIRPA, HAVCR2, and TNFSF9 have been targeted for drug development (Table S5). Most of the drugs are developed for the treatment of cancer and autoimmune diseases.

RT-PCR Validation of the ICGs

The results showed that the relative expression levels of 4 ICGs including Cd44, Cd47, Vtcn1, and Sirpa were consistent with the previous findings. Havcr2 and Tnfsf9 showed no statistically significant difference (Figure 9).

Figure 9 RT-PCR validation of the ICGs between db/db and wild-type mice (n=5). x-axis: ICGs, y-axis: gene expression level; red: db/db mice, blue: wild type. **: P <0.01, ***: P <0.001.

Validation of Identified ICGs and Risk Model

We employed the risk model in T2D PBMCs dataset of GEO. The results showed that the risk score of T2D group was higher than control though not significantly and the AUC = 0.708 (95% CI 0.457–0.959) (Figure S1-2). Also, another independent GEO dataset was used to validate the identified ICGs. The results showed that CD44, HAVCR2, and SIRPA were significantly upregulated in T2D group (Figure S3).

Discussion

The normal function of β cells depends on the islet microenvironment consisting of endocrine cells, immune cells, neurons, and vascular endothelial cells³¹ and immune microenvironment is an important component of islet microenvironment. It has been reported that chronic low-grade tissue inflammation with immune cell infiltration and increased inflammatory cytokines is detrimental to insulin signaling and metabolic homeostasis.³² Early diagnosis and treatment can prevent T2D complications, which is the leading cause of deaths in patients with T2D. Therefore, we used bioinformatics method to investigate the roles of ICGs in the pathogenesis of T2D for the first time.

In the current study, we selected GSE76894 containing pancreatic islets and compared the expression levels of ICGs between T2D patients and normal individuals. We found 458 DEGs, among which 120 genes were upregulated and 338 genes were downregulated and 10 hub genes which were associated with T2D. Six significantly different ICGs between T2D and control group were identified, including CD44, CD47, HAVCR2, SIRPA, TNFSF9, and VTCN1. Functional enrichment analysis showed that the biological process of cytokine–cytokine receptor interaction, defense response, response to cytokine, cytokine mediated signaling pathway and innate immune response is activated in T2D samples, which further implies the potential role of immune pathway in T2D. We proposed a risk model of T2D at gene level based on ICGs, which can predict T2D well. Moreover, close correlations between ICGs and immune cells, several hub genes, T2D-related disease genes, insulin secretion-related genes and β cell function-related genes were observed suggesting that ICGs might get involved in the progression of T2D via regulating infiltrated immune cells and β cell function.

CD44, a cell membrane receptor widely expressed in lymphocytes, fibroblasts, and smooth muscle cells, regulates the adhesion of lymphocytes, T cell activation, and recruitment of leukocytes,³³ thus playing a critical role in inflammatory signaling pathways. CD44 deleted pancreatic β cell line Min6 cells presented a better capacity of insulin secretion as well as insulin content.³⁴ Therefore, CD44 may be important in the normal function of pancreatic β cell and a potential therapeutic target for T2D. In the current study, we found CD44 gene expressions were significantly higher in pancreatic islets of T2D; meanwhile, CD44 was negatively correlated with T cells CD4 memory resting and positively with T cells CD4 memory activated though not significant. Furthermore, CD44 was also correlated with several T2D-related genes. Higher expression of CD44 may promote the progression of activation of T cell resulting in the imbalance of islet microenvironment and worsening function of pancreatic β cell.

CD47 is a multifunctional transmembrane glycoprotein of the immunoglobulin superfamily, and its receptor is signal regulatory protein alpha (SIRPα), an immune inhibitory receptor encoded by SIRPA, broadly expressed across different cells in the body.³⁵ CD47/SIRPα axis is a vital immune checkpoint and an important “marker of self” which prevents immoderate self-phagocytosis by macrophages.³⁵ In the current study, the expressions of CD47 and SIRPA were higher in patients with T2D and CD47 was significantly correlated with SIRPA, suggesting that CD47/SIRPα axis may also be involved in the progression of T2D. There were significant associations between CD47 and T cells CD4 memory resting, SIRPA and T cells CD4 naïve in the current study. A study using nonobese diabetic mouse suggested CD47/SIRPα axis was associated with T cell proliferation and type 1 diabetogenicity.³⁶ Thus, CD47/SIRPα axis may also disturb pancreatic islet immune cell infiltration through its impact on T cells in the progression of type 2 diabetes. We also found a significant correlation between CD47 or SIRPA and several insulin secretion-related and β cell function-related genes, indicating that pancreatic β cell is potentially regulated by CD47/SIRPA axis.

Hepatitis A virus cellular receptor 2 (HAVCR2) encoding T cell immunoglobulin and mucin domain-containing protein 3 (Tim-3), an inhibitory immune checkpoint protein, is a crucial regulatory molecule in immune tolerance.³⁷ A recent study showed Tim-3 on NK cells was overexpressed in T2D patients, and the abundance of NK cells was significantly reduced in the peripheral blood. Another study found that the cytotoxic function of NK cells was impaired in T2D.³⁸ Similarly, our results showed that NK cells activated were lower in T2D patients compared with the control group, whereas NK cells resting were comparatively higher though not significantly. This probably explains why T2D patients are more susceptible to infectious diseases, which contribute to the system or islet local chronic low-grade inflammation to a certain extent. Tumor necrosis factor superfamily 9 (TNFSF9) mainly expressed in antigen presenting cells (APCs), also known as CD137L, is a type II transmembrane glycoprotein with broad and multiple functions.³⁹ Binding to its receptor CD137, TNFSF9 is a bidirectional inflammatory signal.⁴⁰ A study found that CD137 deficient mice fed an HFD presented improved weight gain, glucose tolerance, fatty liver disease, and lowered adipose tissue inflammatory responses. Furthermore, the infiltration of macrophages and T cells into adipose tissue was significantly reduced.⁴¹ In the current study, we found that the expression of TNFSF9 was markedly higher in the islet of T2D patients, implying that it may affect the microenvironment of pancreatic islets as well. However, the inconsistent results of RT-PCR in diabetic mice and bioinformatics analysis on the two ICGs demonstrate that more studies are warranted to identify the role of HAVCR2 and TNFSF9 in T2D.

V-set domain containing T cell activation inhibitor 1 (VTCN1), an inhibitory molecule of the B7/CD28 superfamily, is important in regulating the inducement of tolerance in autoimmune disease by inhibiting T cell activation and function. VTCN1 has been well studied in type 1 diabetes, resulting from the autoimmune destruction of insulin-producing β cells in the pancreas induced by cytotoxic T cells.⁴² However, whether VTCN1 plays a role in T2D, in which T cells were not considered to be the effective cells, was still unclear. VTCN1 protein expression is colocalized with insulin in β cells, and its expression is significantly higher in T2D patients,⁴³ suggesting it may be involved in the process of type 2 diabetes through different mechanisms from those of type 1 diabetes. It was notable that in the current study, VTCN1 was inversely correlated with T2D (coefficient of ridge regression: −0.088) and correlated with insulin secretion-related and β cell function-related genes, suggesting that VTCN1 may play a crucial role in T2D by regulating β cell function and insulin secretion. The exact role of VTCN1 in the pathogenesis of T2D is warranted to be further researched.

In druggability evaluation, we found that CD44, CD47, SIRPA, and HAVCR2 were targeted for drug development, though most of them were developed for cancer or autoimmune disease treatment. However, some drugs may aggravate or induce diabetes. Atkins et al reported that combination avelumab and utomilumab (targeting TNFRSF9) can induce diabetic ketoacidosis,⁴⁴ suggesting that more attention should be paid to the effects on the metabolism system when applying combination immunotherapies. Furthermore, preclinical studies showed that anti-CD44 monoclonal antibody or blocking the HA-CD44 interaction improved insulin sensitivity and reduced body weight.^45,46 Also, our study found that overexpression of B7-H4 (VTCN1) in pancreatic β cell can improve HFD-induced impaired glucose tolerance, and the relevant results are being organized and submitted for publication. Our research indicated that repurposing or developing therapy targeting these ICGs for the treatment of metabolic disorders, such as T2D, may be promising.

In addition to the role of immune regulation, ICGs were also found to directly regulate functions and behaviors of cells that express these molecules.⁴⁷ We have previously reported that Vtcn1 can enhance the differentiation of murine leukemia-initiating cells via the PTEN/AKT/RCOR2/RUNX1 pathways.⁴⁸ Another immune checkpoint molecule B7-H3 has been shown to be positively associated with glycolysis in several studies.^49,50 Therefore, in the current study, we further explored the relationship between ICGs, T2D-related genes, insulin secretion-related genes, and β cell function-related genes. We found that ICGs were associated with several insulin secretion-related or β cell function-related genes, such as GCK, ISL1, and SLC30A8, indicating that ICGs may also affect the progression of T2D through direct regulation of β cell.

To the best of our knowledge, our research is the first study to investigate the associations of ICGs and T2D based on bioinformatics analyses. However, the present study had several limitations. First, our sample size was small and the selected ICGs are needed for further validation in a larger population. Second, our analyses were based on bioinformatic methods, and further research may include clinical samples or external experiments. Third, CIBERSORT tool was developed outside of pancreatic islet context, and the results should be validated in further studies. Finally, the abundance of infiltrated immune cells was calculated via CIBERSORT, which present the relative expressions, so it may be insufficient to reflect the islet microenvironment.

In conclusion, using bioinformatics methods, we compared islet DNA-seq data between T2D and control groups. Six ICGs relevant to T2D were identified. Significant correlation between ICGs, hub genes, T2D-related genes, insulin secretion-related genes, and β cell function-related genes. Our findings provide new clues to the pathogenesis, diagnosis and treatment of T2D. Further experiments are warranted to validate the role of ICGs and pathways in T2D.