The gene matrix for DEP was obtained from the GEO database. Screening was according to the following criteria: (1) the gene sets must include experimental and control groups, (2) DEP patients, and (3) the samples of datasets more than 10. Finally, GPL570 dataset GSE76896, which included 32 non-diabetic samples and 20 DEP samples, was selected as the test set; GPL16791 dataset GSE164416, which included 18 non-diabetic samples and 35 DEP samples, was selected as validation set. In addition, DEP patients were all from metabolically phenotyped pancreatectomized patients with < 1 year duration of diabetes, and their general information is presented in Supplementary Table S1 . We normalized the GSE76896 microarray using the robust multiarray averaging (RMA) algorithm. The “limma” package in the R software was used to screen the DEGs in GSE76896 (13). The parameters were set to adj.P.Val < 0.05 and |log2FC| ≥ 0.5.

WGCNA is a common method for systems biology analysis, which identifies key genes by linking gene networks to clinical traits. The top 25% of genes with the largest variance in all samples were selected as the input dataset for the subsequent WGCNA. First, to guarantee that the constructed network was valid and reliable, we removed the outlier cases by hierarchical clustering analysis. Then, the appropriate soft threshold was determined by scale independence and mean connectivity, and the correlation matrix was converted into an adjacency matrix. Subsequently, the adjacency matrix was further processed into a topological overlap matrix (TOM). The dynamic tree cutting method was performed by hierarchical clustering of genes to identify modules. The minimum size cutoff of 30 and a deepSplit value of 2 were chosen as the distance measure for the dendrogram. Eigengene was calculated to merge comparable modules. Finally, the WGCNA package was used for analysis and preservation (12).

To select meaningful modules, Pearson’s correlation test was applied to evaluate the correlation between features and gene modules. Then, we defined the module–gene correlation as module membership (MM) and the trait–gene correlation as the gene significance (GS). Finally, we further analyzed the correlation between MM and GS in the selected meaningful modules.

The intersection of the DEGs was taken with the genes in our selected module and as a new set of genes. The results were visualized by using Venn diagram. The PPI network according to all Co-DEGs was built through the online STRING. Next, we imported the interaction results into Cytoscape software (v3.7.2) to construct the PPI network for better visualization. Then, the Minimal Common Oncology Data Elements (MCODE) was performed to extract significant gene clusters and acquire cluster scores (screening criteria: max depth=100; k-core=2; node score cutoff=0.2; degree cutoff=2). Gene clusters with scores >5 were focused. In addition, hub genes in the network was identified by cytoHubba (14). Three algorithms were used, namely, Maximal Clique Centrality (MCC), Maximum Neighborhood Component (MNC), and Degree, to identify the top 10 core genes. Finally, the results of the three algorithms were crossed to get the final hub genes.

To further explore the characteristics and functions of Co-DEGs, we used Gene Ontology (GO) annotations to characterize biological properties, including molecular function, cellular component, and biological process (15). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was conducted to determine functional attributes (16). The threshold for significant terms was set as p < 0.01.

CIBERSORT is an important tool for measuring immune infiltration and can be utilized to assess the abundance of member cell types in mixed cell populations based on gene expression (17). Immune cell signatures for each sample were obtained from GSE76896 using CIBERSORT with the number of permutations set to 100. Wilcoxon rank-sum test was then used to compare the proportion of islet immune cells in the normal and DEP groups, and the result was presented in a box plot.

The hub genes that have not been studied in DEP were further validated in another dataset GSE164416 downloaded from the GEO database. In the GSE164416 dataset, “EdgeR” was used for linear correction of expression profiles. Then, the corrected gene matrix was extracted by “Limma Voom.” In addition, ROC curves for selected genes were calculated using “pROC” package via R software to assess the ability of hub genes to diagnose DEP.

miRNAs are small non-coding RNAs that have been shown to regulate gene expression by inhibiting mRNA translation or promoting mRNA degradation (18). Therefore, miRNet (version 2.0) was used to predict target miRNAs of hub genes. The Human microRNA Disease Database (HMDD) is a database that collects experimental evidence for the associations of human miRNA with disease (19). We acquired DM-associated miRNAs and intersected with miRNAs predicted by hub genes. To explore the function of miRNAs, GO BP analysis was performed using mirPath v3.0 software in DIANA Tools (20). A p-value <0.01 was used to screen for significant GO terms. The miRNA–mRNA interaction network according to the association between miRNAs and mRNAs was constructed by using Cytoscape.

The R software was used for statistical analysis, and the ggplot2 package was used for box plots. Data were analyzed using SPSS 23.0, and Student’s t-test or Wilcoxon rank-sum test was performed to compare the differences between the two groups. A two-tailed p < 0.05 was assumed statistically significant.

After matching the sample traits with the expression matrix, the sample cluster tree is shown in Figure 2A ; meanwhile, 5 was selected as the soft-threshold power according to the scale independence and mean connectivity values ( Figure 2B ), which indicates that a scale-free network has been successfully constructed. A total of 12 modules were divided by WGCNA, which were indicated by different colors, and the larger the area of the module, the higher the number of genes included ( Figure 2C ). In addition, based on TOM, the network heat map plot between all genes is shown in Figure 2D , indicating that these modules were independent of each other ( Figure 2D ).

The relationship between modules and trait indicated that several modules were associated with DEP, and it clearly showed that the black module (r=0.71, p=3e−07) was most significantly associated with DEP ( Figure 3A ). In addition, DEP was combined with modules to redo the clustering tree and to make a heat map of Mes correlations. The shades of the overlapping colors of the heat map suggested the magnitude of the correlation between the different modules ( Figure 3B ). Moreover, Figure 3C shows that these genes in the black module had a high correlation with DEP, and the different modules were independent of each other (cor=0.58, p=2.4e−84) ( Figure 3C ).

To further identify genes with differences in the black module, we used the limma package processing dataset GSE76896 to identify DEGs. Compared with non-diabetic samples, a total of 1,129 DEGs were identified in the DEP samples; volcano plot was used to visualize these DEGs ( Figure 3D ). The top 20 upregulated genes and 10 downregulated genes identified in the DEGs are shown in Supplementary Table S2 . Next, The Co-DEGs were obtained from the intersection of the black module and the DEGs, including 373 genes, and the result is shown in Figure 3E .

We uploaded 373 Co-DeGs to STRING for PPI analysis, resulting in a total of 255 nodes and 1,636 edges. The result of interactions was then further analyzed using Cytoscape as shown in Figure 4A . Moreover, the color shades of the node indicated its level of neighborhood connectivity, and the magnitude of the node indicated the value of degree. Network data was processed with MCODE to identify gene clusters, and clusters with scores >5 were focused. A total of three gene clusters with high scores were screened, and their scores were 23.643, 13.500, and 6.556, with cluster 1 being the most important ( Figure 4B ). We used the cytoHubba in Cytoscape to determine hub genes and identified five hub genes by crossing the results of the three algorithms, including Degree, MCC, and MNC ( Figure 4C ). Details of the five hub genes are shown in Table 1 .

In order to explore the signatures of Co-DEGs, we performed GO annotation and KEGG pathway analyses ( Figures 5A, B ). The results of GO annotation suggested that these differential genes were mostly enriched in processes involved in immune effector process, regulation of macrophage migration, and toll-like receptor signaling pathway. Considering the KEGG pathway, the majority of genes were enriched in cell adhesion molecules (CAMs), cytokine–cytokine receptor interaction, and Toll-like receptor signaling pathway. For each gene cluster and hub genes, we chose several keywords to describe their major biological functions. Among them, hub genes were mainly related to macrophage activation ( Figure 5C ).

We detected the infiltration of immune cells in the islet microenvironment by CIBERSORT. The bar chart showed the proportion of immune cells in all samples ( Figure 6A ). Obviously, macrophages accounted for the majority of all infiltrating cells. Then, we compared the differences in immune-infiltrating cells between the normal group and the DEP group ( Figure 6B ). Five types of immune cells, namely, memory B cells, M1 macrophages, eosinophils, follicular helper T (Tfh) cells, and activated natural killer (NK) cells, were differentially expressed (all p<0.05). We focused on the correlation between TLR4 and immune-infiltrating cells ( Figure 6C ). Among them, six immune cells were positively correlated with TLR4, and seven immune cells were negatively correlated with TLR4. Also based on the GSE76896 dataset, we found that MYD88 expression was elevated in DEP (p<0.05), while TRIF expression was unchanged ( Figures 6D, E ). What is more, we detected the correlation between several other hub genes and immune-infiltrating cells, and both of them were positively correlated with M1 macrophages ( Supplementary Figure S1 ).

To verify screened genes, we analyzed the expression of hub genes in the GSE164416 dataset. As shown in Figure 7 , the expression of TLR4, ITGAM, ITGB2, PTPRC, and CSF1R was significantly upregulated (all p < 0.05) in the islets of DEP compared to the control group in datasets GSE76896 and GSE164416 ( Figures 7A, B ). What is more, we plotted the ROC curves of hub genes and calculated the area under the curve (AUC) to distinguish DEP from the control group. Among the five hub genes, TLR4 has the highest diagnostic value. The diagnostic values of hub genes in GSE76896 are follows: TLR4 (AUC, 0.8578), ITGAM (AUC, 0.7797), PTPRC (AUC, 0.7766), ITGB2 (AUC, 0.7594), and CSF1R (AUC, 0.7531) ( Figure 7C ). In addition, they also have high diagnostic value in the GSE164416 dataset ( Figure 7D ).

By searching the HMDD database, we found that 71 miRNAs were related to DM ( Supplementary Table S3 ). Figure 8A shows that there were 15 common miRNAs in diabetes-related and predicted miRNAs ( Figure 8A ). Then, the signatures of these 15 miRNAs were investigated. The enrichment analysis revealed that three biological processes, namely, cellular nitrogen compound metabolic process, cellular protein modification process, and neurotrophin TRK receptor signaling pathway, were regulated by all 15 miRNAs. Interestingly, these results also included “toll-like receptor 4 signaling pathway,” which suggested that those common miRNAs involved in DM could also regulate the TLR4 signaling pathway ( Figure 8B ). This confirmed our analysis again.

Based on the prediction results of miRNet database, an miRNA and mRNA co-expression regulatory network was constructed by Cytoscape, which consisted of 44 nodes and 44 edges ( Figure 9 ). Among these miRNAs, the blue color represents the miRNA associated with DM. Combined with the enrichment results of miRNA, TLR4 signaling pathway may be regulated by has-miR-155-5p, has-miR-27a-3p, and has-miR-21-5p. Thus, we hypothesized that has-miR-155-5p/has-miR-27a-3p/has-miR-21-5p-TLR4 might lead to TLR4 signaling pathway activation in DEP.