Single-cell transcription data were obtained from the GEO database; construction of the GSE165816 dataset was based on the GPL24676 platform. The Seurat package (version 4.0.2) in R software was taken for cell cluster, following t-distribution randomized neighbor embedding (t-SNE) analyses plus principal component analysis (PCA). After cells with <200, >2,500, or >5% mitochondrial genes were filtered out, 32,574 filter cells were chosen for study. Logarithmic normalization was performed to normalize the gene expression, which was further scaled. The “VST” method was then used on every specimen to identify 2000 hypervariable genes (HUGs). Following that, principal component analysis was carried out to identify important principal components (PC), with the JackStraw and ScoreJackStraw functions used to visualize p-value distribution. Batch correction was achieved using the R package “Harmony” (edition .1.0). Based on these results, 18 PCs were used for the t-SNE testing. The cells were grouped using the FindCluster function according to 0.2, 0.3, 0.4, and 0.5, respectively. Based on the results of every group, the final resolution was determined as 0.4, and the cells were classified into 21 distinguished clusters. The differentially expressed genes (DEGs) for every cluster were found using the FindAllMarkers function of Logfc.Threshold = 0.25 as a guide. The cell types were then manually annotated in accordance with other research.

Screening the DEGs of every cluster revealed IRGs from the ImmPort database (https://www.immport.org/shared/home); 1,509 IRGs inside the DEGs were then chosen for IRG rating using AUCell (Version 1.12.0). Based on a gene set enriching study, the AUCell R program ranked the pathways for every cell. To assess the percentage of the chosen 497 IRGs that are strongly expressed in every cell, gene expression rankings for every cell were created based on the area under the curve (AUC) scores of those IRGs. Cells with higher AUC values were those that expressed more genes from the gene set. The threshold for identifying cells with active gene sets was determined using the “AUCell explore Thresholds” tool. The ggplot2 R package (Version 3.3.5) was then used to map every cell’s AUC score to the UMAP embedding in order to display the active clusters.

GEOquery software was used to obtain the raw data for GSE80178 from the GEO database (Version 2.58.0). Limma software was used to compute DEGs (edition 3.46.0). Significantly dysregulated genes were defined as DEGs with an absolute logFC more than 1 and an adjusted p-value less than .05. With the help of the ggplot2 program, volcano and heatmap plots were created (edition 3.3.5).

Intersection for the M1 macrophage cluster in GSE165816 and DEGs in GSE80178 yielded 131 differentially expressed IRGs, which were then introduced into the online bio-informatics website (https://maayanlab.cloud/Enrichr/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) for Gene Ontology (GO) analysis, and the first 10 pathways were screened out according to the P-score order .

Using STRING (https://string-db.org), a network analysis of protein–protein interactions (PPIs) was conducted. The TSV format file was download, which was imported into Cytoscape software for visualization, and the Cytoscape program was employed to filter hub TFs using the cytoHubba plug-in.

By using the CIBERSORT algorithm, the expression data for GSE80178 were analyzed, and the proportion of 22 types of immune cells was determined. These immune cells included immature B cells, memory B cells, plasma cells, CD8+T cells, immature CD4+T cells, resting memory CD4+T cells, memory cells with activated function, follicular helper T cells, regulatory T cells, γ-incremental T cells, and resting NK cells. The R software’s vioplot tool was employed to compare the different immune infiltration levels of every immune cell between the two cohorts.

The DFU cohort of 22 patients with DFUs admitted to our hospital between March 2021 and January 2022, comprising 12 men and 10 women with a mean age of 56.12 ± 12.93 years and an average diabetes duration of 3–24 years. All patients had no heart, liver, kidney, or malignant tumors, and none had ever had a cerebral hemorrhage, a cerebral infarction, or an acute myocardial infarction (AMI). The healthy control (HC) group consisted of 12 healthy people who received physical exams at our hospital within the same time frame. The average age of this group of nine men and three women was 56.12 ± 12.93 years. There was no discernible difference between the three groups in terms of age or gender ratio (all p more than .05). A heparin anticoagulant tube containing 10 ml of every subject’s fasted blood was collected from their cubital vein in the morning. The serum was then separated and kept at −80°C after being centrifuged for 20 min at 3,500 rotations per minute. Six cases of foot skin tissue from DFUs and six HC cases were chosen and fixed with formalin. Immunofluorescence detection was then used to identify the disease.

After washing, the skin tissue blocks were cut into small pieces and placed in a homogenizing tube for thorough homogenization. The specimen tube was taken out after homogenization, and ice bath and shock were carried out to ensure complete cleavage of the tissue. After centrifugation, the supernatant was collected to form the total protein solution. We prepared separation glue, added TEMED, shook well, and poured the glue. After solidification, electrophoresis began, and the membrane was then transferred. The primary antibody was added and incubated overnight. The secondary antibody was added after washing. Following incubation and washing, chemical luminescence, development, and fixation were performed. The optical density score of the target band was examined using an alpha program processing system—the protein expression level.

Using the TRIzol reagent, total RNAs from the plasma were extracted and reverse-transcribed into cDNA. The PCR reaction comprised of 40 cycles of denaturation at 95°C for 10 min, annealing at 60°C for 1 min, and extension at 95°C for 15 s. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was utilized as the internal standard. The levels of related factor expression were determined using a PCR reaction with 40 cycles of denaturation at 94°C for 15 min, annealing at 55°C for 30 s, and extension at 70°C for 1 min. GAPDH was used as an internal reference. Each experiment was conducted a minimum of three times. Using the 2^−ΔΔCT (cycle threshold (CT)) formula, the relative concentrations of NFE2L2, REL, ETV6, MAF, and NF1B were calculated. All the primers needed were purchased from Servicebio.

Using 13 antibodies coupled to either cyanine 3 (cy3) or cyanine 5 (cy5) and 4′,6-diamidino-2-phenylindole (DAPI), eight staining cycles were performed on FFPE tissue microarrays. The antibodies employed were NRF2, ETV6, RELB, and NF1B GB113808, GB112082, GB11985 and bs-11899R are the numbers of antibody primer information of NRF2, ETV6, RELB, and NF1B, respectively, provided by Servicebio Laboratories, Inc. provided fluorescent-labeled secondary antibodies for use in the research. The paraffin sections were removed, put in water, and repaired in a microwave-safe repair box using citric acid (PH6.0) antigen repair solution. The primary and secondary antibodies were applied to the sections after they had been gently dried, and a circle was drawn around the tissue with a tissue chemical pen. The sections were then incubated for 50 min at room temperature in the dark. The slides were submerged in PBS (PH 7.4) and subjected to three 5-min decolorization washings. After 5 min, the ring was treated with autofluorescence quencher and washed under running water for 20 min. Anti-fluorescence quenched tablets were then used to seal the DAPI counterstained nuclei. Afterward, the slices were examined using a fluorescence microscope, and pictures were recorded and gathered.

The ScRNA sequencing dataset (GSE165816) was analyzed; after filtering, 22,009 feature analysis data of 32,574 cells were retained. Figure 1A shows the precise expression of every specimen. With a correlation coefficient of 0.83, the nCount RNA— representing the number of different molecular identifiers (UMI)—and the nCount_RNA— representing the number of genes—exhibited a positive association (Figure 1B). The top 10 HVGs were determined (Figure 1C). The top two HVGs—IGKC and IGHG1—are tiny calcium-binding proteins that are highly expressed in inflammatory situations. Using JackStrawPlot, all 20 PCs were detected by PCA with a p-value less than 0.05 (Figure 1D). Based on the previous research, the PCs were clustered with a selected number with the resolution determined, and all cells were finally divided into clusters. The differences of clusters before the annotation were analyzed, with the markers of all clusters determined. Eventually, 10 PCs were used to identify 16 clusters, and the first 10 DEGs for every cluster were presented (Figure 1E).

T-SNE analysis allowed the 16 clusters to be visualized (Figure 2A). The proportion of M1-Macro clusters in the DFU group decreased in the ellipse-designated M1-Macro clusters (Figure 2B). The dot diagram (Figure 2C) plus violin chart shows the expression of cell-type marker genes (Figure 2D).

The DEG of every cluster was screened according to the imported database to generate IRG and explore the IRG expression characteristics of skin cell groups in patients with DFUs. The IRGs in published studies were summarized in this database, and 212 IRGs were involved from the DEG of every cluster in the skin cell group. Using the AUCell R program, the IRG activity of every cell line was assessed (Figure 3A), and the AUC score of M1-Macro and M2-Macro was found to be significantly higher than other sub-populations of cells. Cells with more gene expression showed higher AUC values, which were mainly located in M1-Macro and M2-Macro in yellow in Figure 3B. Every specimen’s cell count and percentage are shown in detail in Figure 3C. Compared to the healing group, M1-Macro clusters were significantly reduced in the non-healing group (20% vs. 71%). In the subsequent analysis, we focus on M1-Macro clustering was mainly focused on (Supplementary Table S1). We also conducted GO and KEGG enrichment assessments on M1-Macro differentially expressed IRGs from skin tissue, which is mainly concentrated in immune-related pathways (Figures 3D,E).

In terms of the expression characteristics of DFUs skin tissue, the high-volume RNA sequencing dataset GSE80178 from six foot-skin DFUs, three forearm-skin DFUs, and three control subjects were analyzed; the immune-associated biomarkers of DFUs were screened out by differential analysis. Those with | logfc | > 1 and corrected p-values less than 0.05 DEG (Supplementary Table S2) were chosen. A total of 577 up-regulated and 2,174 down-regulated DEGs in total were kept in DFUs (Figure 4A). The top 30 up-regulated DEGs and the top 30 down-regulated DEGs are shown as heat maps to display a relative consistency among groupings (Figure 4B). Interestingly, the top 10 GO keywords of DEGs in DFU skin large data were likewise centered on immune responses, which is congruent with the expression features of the M1-Macro cluster in DFU single cells (Figures 4C,D). The outcomes suggest a potential role for M1 macro-cells in foot skin tissue from DFU patients.

By running the CIBERSORT algorithm, the difference in immuno-osmosis between DFU and healthy skin tissues was assessed within 22 immune cell subtypes. Compared with healthy tissues, the number of M1 macrophages in DFU specimens were significantly low, which confirmed the conclusions of single-cell analysis (Figure 4E, p < 0. 05).

We subsequently investigated the co-expression features of IRGs, taking into account the GO for IRGs in IRGs and DFU skin—which is mostly focused on immunological response in the skin tissue cell M1 macrophage cluster and DFU skin. Some 1,665 TFs were obtained from HumanTFDB to study the transcriptional regulatory activity of IRGS. A total of 26 TFs were acquired by intersecting with the MI macrophage cluster and IRGs in DFU skin (Figure 5A). The default arguments for the FindAllMarkers function were only.pos is true, min.pct is equal to 0.25, and logfc.Threshold is equal to 0.25. The expression of the 26 TFs in GSE165816 and GSE80178 is shown in Figures 5B,C, among which three factors were up-regulated in DFUs, and all 26 TFs were expressed in M1 macrophages. The different IRGs of the scRNA analysis portion of the M1macro and the bulk RNA-seq part were intersected to obtain 25 common IRGs (Figure 5D and Supplementary Table S3). Their expression is shown in Figure 5E. The 26 co-expressed TF were input into STRING database to construct a PPI network. The cytoHubba plug-in in Cytoscape was analyzed by hub, with the different methods ranked as follows: MultiClass Classification, DMNC, MN, Degree, EPC, EcIcentrity of Bottle Necks, Radiality of Closeness, Betweenness, the Stress factor, and ClusteringCoefficient. The top five with a relatively greater frequency were selected hubTFs for constructing the PPI network. NFE2L2, REL, ETV6, MAF, and NF1B may play a key role in the transcriptional regulation of IRGs as central genes (Figure 5F).

Plasma tissues from 22 DFUs patients and 12 healthy subjects were clinically collected to further support the hypothesis that the genes screened previously are associated with the immune microenvironment of DFUs. The expression levels of NFE2L2, REL, ETV6, MAF, and NF1B genes were first discovered in the blood by WB and qPCR. We noticed that the blood of DFU patients had significantly higher expression of NFE2L2, REL, ETV6, MAF, and NF1B (Figures 6A,B). qPCR analysis of patients’ blood pressure also showed that immune-related factors were significantly higher than those of healthy people, with a statistical significance (Figure 6C and Supplementary Table S4).

A number of immune-associated genes were also identified in skin tissue taken from six DFU patients and healthy subjects using immunofluorescence. The expression of NFE2L2, REL, ETV6, and NF1B in skin tissues of healthy subjects and DFU patients were analyzed by labeling under a fluorescence microscope. The findings demonstrated that NFE2L2, REL, ETV6, and NF1B transcription in DFU skin was significantly higher than that in healthy human skin. This was consistent with the results of our single-cell sequencing, which showed that IRGs may be involved in the immune function of DFUs (Figure 7).