Identifying and Validating GSTM5 as an Immunogenic Gene in Diabetic Foot Ulcer Using Bioinformatics and Machine Learning

Introduction

Patients with newly diagnosed diabetes or a history of diabetes who develop peripheral arterial disease, peripheral neuropathy, infection, ulcers, neuro-osteoarthropathy, gangrene, or amputation of any kind in the foot are considered to have diabetes-related foot diseases. Diabetic foot ulcer (DFU) is one of the most serious complications. Patients with diabetic foot ulcers usually have one or more risk factors, among which diabetic peripheral neuropathy and peripheral arterial disease play important roles.¹ DFU represents one of the most catastrophic outcomes of diabetes and stands as a considerable healthcare challenge. It imposes serious economic, physical, and psychological burdens on patients.² Extensive alterations related to diabetes, such as neuropathy and vascular conditions, frequently worsen the advancement of DFU.³ The lifetime occurrence of DFU among individuals with diabetes varies between 15% and 25%. Furthermore, in cases where there is a prior history of foot injuries or infections, the occurrence of DFU in diabetic patients rises to a range of 19% to 34%, with approximately 17% of these cases ultimately necessitating amputation.⁴ In the year 2019, worldwide direct healthcare costs associated with diabetes amounted to around $700 billion. Projections suggest that by 2030, these expenses are anticipated to rise to $825 billion, with the medical costs linked to DFU constituting a third of the overall expenditure for managing diabetes.⁵ The development of DFU involves a notably intricate process, frequently triggered by distal neuropathy and peripheral vascular disorders in the lower limbs. In over half of DFU instances, infections are also present, adding complexity and difficulty to the treatment process.⁶ Therefore, it is crucial to study new biomarkers of DFU and inhibit the development of chronic wounds in order to improve treatment efficiency and prognosis for DFU patients.⁷

The slow onset and delayed resolution of inflammatory response are among the major factors contributing to the difficulties in wound healing.⁸ Immune cell infiltration is also a pivotal factor in the onset and progression of DFU. When viewed from a molecular standpoint, the process of wound healing takes place following a disruption in the skin’s protective barrier. This process is typically facilitated by growth factors and cytokines released by specialized cells that become activated as part of the immune response.⁹ Sawaya et al observed that in cases of DFU, the transcription factors FOMI1 and STAT3, which support the survival of immune cells, are suppressed. This ultimately hinders the process of wound healing in diabetic patients.¹⁰ Through this study, our objective is to clarify the potential of immune-related biomarkers in enhancing the effectiveness of DFU treatment.

The high-throughput transcriptome sequencing data and clinical annotations provided by the DFU Initiative enable us to investigate the changes in transcriptional patterns and associated molecular pathways relevant to DFU in biological research.¹¹ Several research investigations have made use of gene expression data extracted from the Gene Expression Omnibus (GEO) to explore the molecular markers associated with the progression of DFU.^7,12 Due to the progress in high-throughput sequencing technologies and the integration of machine learning in the medical field, novel methodologies have surfaced for examining molecular targets in diverse diseases.^11,13 Through the application of bioinformatics techniques, we acquired a gene expression dataset for DFU from the GEO database. In order to pinpoint biomarkers linked to DFU, we employed two machine learning algorithms: LASSO and SVM-RFE. Following that, candidate genes strongly associated with immune infiltration were validated using a separate, independent validation dataset. The glutathione S-transferase (GST) gene family is one of the major xenobiotic detoxification enzymes that can protect cells from the harmful effects of toxic drugs and environmental electrophilic agents. GSTM5 has been found to be dysregulated in various tumors.¹⁴ However, there is no relevant research indicating its role in DFU. Therefore, in this study, we constructed a mouse model of DFU and validated the GSTM5 gene. The flowchart of this study is shown in Figure 1.

Figure 1 Research framework.

Materials and Methods

Raw Data

Retrieval of microarray data on mRNA expression profiles related to DFU from the GEO database. This study utilized the GEO datasets GSE68183, GSE80178, and GSE134431, with their respective platform files GPL16686 and GPL18573. GSE68183 and GSE80178 are used as training datasets, while GSE134431 is used as a validation set. Detailed information about the aforementioned datasets can be found in Table 1.

Table 1 Dataset Information

Data Filtering and Processing

Transform the downloaded probe matrix into a gene expression matrix by referencing the probe annotation file. In cases where a gene corresponds to multiple probes, calculate the mean value of those probes to represent the ultimate expression level of the gene.¹⁵ After standardizing datasets GSE68183 and GSE80178, the SVA package is used to address batch effects between GSE68183 and GSE80178. Principal component analysis (PCA) is employed to evaluate the removal of batch effects. Subsequently, the limma package is utilized to analyze the differential expression matrix comparing the control group and DFU. To identify genes related to immune infiltration in DFU patients, we employ the criteria of |log FC| > 1 and an adjusted p-value < 0.05¹⁶ as thresholds for detecting differentially expressed genes (DEGs).

Machine Learning

The LASSO method is used to identify essential features from large-scale data and screen gene models. The “glmnet” R package is utilized to implement this method in the process.¹⁷ SVM-RFE is implemented through the “e1071”¹⁸ R package and has been widely used as an effective feature selection algorithm for the selection of central diagnostic genes.¹⁹ The genes derived from the aforementioned two machine learning models can be mapped and employed for subsequent analysis. The Receiver Operating Characteristic (ROC) curve evaluates the diagnostic utility of these biomarkers.

Pathway and Functional Enrichment Analysis

Apply Kyoto Encyclopedia of Genes and Genomes (KEGG)²⁰ and Gene Ontology (GO)²¹ enrichment analysis using the R package clusterProfiler.²² Additionally, we performed Gene Set Enrichment Analysis (GSEA) to identify potential pathways.²³ Gene Set Variation Analysis (GSVA) is an unsupervised, non-parametric analytical technique primarily employed to evaluate the enrichment of gene sets in sequencing data. GSVA enables the assessment of potential alterations in pathway activity within each individual sample.²⁴

Immunocyte Infiltration Analysis

The relative abundance of 22 lymphocyte subtypes in each DFU and normal sample was calculated based on the CIBERSORT algorithm and the immune cell LM22 gene set.²⁵ The Wilcoxon test was used to compare the differences in 22 immune cell types between DFU and normal samples.

Correlation Analysis Between Diagnostic Biomarkers and Immune Cells

Performing Spearman correlation analysis to evaluate the correlation between diagnostic biomarkers and infiltrating immune cells.²⁶

Establishment and Treatment of DFU Mouse Model

Male C57BL/6 mice of adult age (21–25 g, 8–10 weeks old) were procured from Shanghai Model Organisms Center Inc. These mice were accommodated at the Experimental Animal Center of Shanghai University of Traditional Chinese Medicine in an environment that met specific pathogen-free standards. They were maintained on a 12-hour light/dark cycle and provided unrestricted access to both water and food.²⁷ The mice were allocated randomly into distinct groups (n=4), which included the non-diabetic group treated with phosphate-buffered saline (PBS), and the DFU group. For a duration of three weeks, the DFU group of mice were fed a high-fat diet constituting 60% of their total calorie intake (FB-D12451, obtained from Wuxi Fan Bo Biotechnology Co., Ltd.), and within one week, they were subjected to intraperitoneal injections of STZ (40 mg/kg, Cat. No. 2196GR001, BioFRoxx).²⁸ Until the level of fasting blood glucose surpassed 16.7 mmol/L, it was classified as type 2 diabetes.²⁹ Following that, a full-thickness skin wound measuring 1*1cm2 in size, penetrating down to the fascial layer, was generated on the dorsal surface of the mouse by lifting the skin with forceps.³⁰ The wound was covered with a 1cm² gauze saturated with a wound dressing solution and replaced daily. To assess the condition of the dorsal wounds in mice, photographs of the wounds were taken on days 0, 7, and11 using a digital camera (Nikon, To to).³¹ Quantification of wound closure rate was performed using ImageJ software (Bethesda, MD) and calculated as follows: t-t₀/t₀ × 100% (t: evaluation of wound healing, t0: initial injury).²⁷

Hematoxylin and Eosin and Masson Staining

The wound tissue was preserved in 4% paraformaldehyde (PFA, catalog number P0099, Beyotime) for 48 hours. Following standard paraffin embedding and sectioning procedures, the wound tissue underwent staining with hematoxylin and eosin (H&E, catalog number C01105M, Beyotime, China) as well as Masson’s trichrome stain (catalog number G1340, Solarbio Life Sciences, China). The stained tissues were subsequently examined utilizing a digital slide scanning system (point M8, Precipoint).

Immunohistochemistry Staining

The ulcerated tissues from both groups of mice were gathered and preserved in 4% paraformaldehyde for a duration of 24 hours. Subsequently, these tissues were encased in paraffin wax and later sliced into sections with a thickness of 5μm. Immunohistochemistry (IHC) analyses were conducted on these paraffin-embedded sections. The sections were subjected to an overnight incubation at 4°C with GSTM5 antibodies. Afterward, the slides designated for IHC staining underwent incubation with the secondary antibody (catalog number GB23303, Servicebio, diluted at 1:200). Visualization was achieved using DAB chromogen, followed by counterstaining with hematoxylin. Images were recorded using an optical microscope (Olympus BX41, Shanghai, China).

Statistical Analysis

Statistical assessment of disparities in gene expression between the two groups was conducted utilizing the ggpubr package within R software. All gene expression data and statistical analyses were executed using R version 4.3.1. If the data conforms to a normal distribution, a two-sample independent t-test is employed to compare the two groups. In cases where the data does not conform to a normal distribution, the Wilcoxon test is utilized for between-group comparisons. A p-value below 0.05 is regarded as statistically significant.

Ethics Approval and Consent to Participate

All animal experiments were approved by the Review Board of Animal Ethics Committee of Shanghai University of Traditional Chinese Medicine (Approval No.: PZSHUTCM220711028). The study is reported in accordance with ARRIVE guidelines (https://arriveguidelines.org). The animal experimentation in this study was conducted in full compliance with the National Institutes of Health guide for the care and use of Laboratory animals (NIH Publications No. 8023, revised 1978).

Results

GEO Data Processing

We have integrated two DFU datasets, GSE68183 and GSE80178, which collectively contain 6 normal skin samples and 12 DFU samples. As shown in Figure 2, the gene expression levels and PCA of each sample before (Figure 2A) and after (Figure 2B) eliminating batch effects.

Figure 2 Dataset preprocessing and difference analysis. (A) Principal component analysis (PCA) analysis before batch effect removal; (B) PCA analysis after batch effect removal; (C) Heatmap of differential analysis between diabetic foot ulcer (DFU) group and control group; (D) Volcano plot of differential analysis between DFU group and control group.

Identifying DEGs Results in DFU

Differential analysis results indicate that there are a total of 234 DEGs in this study, with 59 upregulated and 175 downregulated. DEGs can be visualized through heatmaps (Figure 2C) and volcano plots (Figure 2D).

Screening Disease-Specific Genes

In this study, we employed LASSO and SVM-RFE models to identify DFU feature genes. In the LASSO model, the λ value of 7 was selected as the optimal variable selection criterion, resulting in the identification of genes LGR5, GSTM5, LCE3D, TRMT1L, XDH, and SRY (Figure 3A and B) as the feature genes. Based on the SVM-RFE algorithm, genes GSTM5, CBR4, KLHDC1, ZBED8, ENPP5, ZNF248, CCDC15, TRMT1L, MATN2, CSTA, TMEM45B, NRN1, and FLRT2 (Figure 3C and D) were identified as the feature genes. Subsequently, the above two gene models were mapped, resulting in the final identification of genes GSTM5 and TRMT1L (Figure 3E).

Figure 3 Construction of machine learning models. (A) Construction of Lasso model; (B) Variation of hub gene coefficients with different logarithmic λ values; (C) and (D) Construction of model SVM-RFE; (E) Mapping between Lasso and SVM-RFE models.

Analysis Disease-Specific Genes

In the training set, the violin plot reflects the expression levels of genes GSTM5 and TRMT1L in the DFU group and the control group. As shown in the figure, both genes GSTM5 (Figure 4A) and TRMT1L (Figure 4B) are down-regulated in the DFU group. The line graph illustrates the expression profiles of genes GSTM5 and TRMT1L across various samples (Figure 4C). Furthermore, the ROC curve demonstrates that both genes GSTM5 (Figure 4D) and TRMT1L (Figure 4E) have good diagnostic performance.

Figure 4 Differences in the intersection model within the training set. (A) The differential expression of gene GSTM5 between the DFU group and the control group; (B) The differential expression of gene TRMT1L between the DFU group and the control group; (C) The expression profiles of genes GSTM5 and TRMT1L in each sample; (D) The ROC curve of GSTM5; (E) The ROC curve of TRMT1L. Note: ***P<0.001.

In the validation group, only GSTM5 (Figure 5A) showed significant down-regulation in the DFU group. However, there were no differences observed between the two groups in terms of the TRMT1L gene (Figure 5B). Furthermore, the ROC curve also indicates that the diagnostic performance of GSTM5 (Figure 5C) is superior to TRMT1L (Figure 5D).

Figure 5 Differences in the intersection model within the validation set. (A) The differential expression of gene GSTM5 between the DFU group and the control group; (B) The differential expression of gene TRMT1L between the DFU group and the control group; (C) The ROC curve of GSTM5; (D) The ROC curve of TRMT1L. Note: ***P<0.001.

Analysis of Gene Differential Expression and Enrichment Analysis in the Model

Based on the high and low expression of GSTM5, we conducted differential analysis of the genes in the training set. The results of the differential analysis were visualized as a heatmap (Figure 6A) and correlation heatmap (Figure 6B). Enrichment analysis indicates that the genetic biological processes of differentially expressed genes are concentrated in cell-substrate adhesion (GO:0031589), epithelial cell proliferation (GO:0050673), and negative regulation of leukocyte apoptotic process (GO:2000107) (Figure 6C). KEGG enrichment analysis indicates that the differentially expressed genes are primarily involved in NF-κB signaling pathway (hsa04064), Nitrogen metabolism (hsa00910), and GnRH signaling pathway (hsa04912) (Figure 6D). These results indicate that GSTM5 is closely associated with inflammation. GSEA analysis indicates that in the high-expression group of GSTM5, butanoate metabolism, complement and coagulation cascades, and drug metabolism cytochrome p450 are active (Figure 6E). In the low-expression group of GSTM5, insulin signaling pathway, MAPK signaling pathway, and neurotrophin signaling pathway are active (Figure 6F). GSVA analysis also indicates a close association between GSTM5 and terpenoid backbone biosynthesis and complement and coagulation cascades (Figure 6G).

Figure 6 Differential analysis and enrichment analysis of genes related to GSTM5. (A) Differential analysis heatmap; (B) Correlation analysis heatmap; (C) Differential gene GO analysis; (D) Differential gene KEGG analysis; (E) GSEA analysis of downregulated genes; (F) GSEA analysis of upregulated genes; (G) GSTM5-related GSVA analysis.

Analysis of Immune Function and Immune Cell Infiltration

The immune function analysis suggests that GSTM5 may have a certain relationship with B cells and Cytolytic activity (Figure 7A). In DFU samples, B cells naive and Dendritic cells activated are higher compared to the control group, while T cells CD4 memory resting, T cells regulatory (Tregs), T cells follicular helper, Macrophages M1, and Mast cells resting are lower compared to the control group (Figure 7B). The distribution of immune cells and the PCA analysis are shown in the Figure 7C. The distribution of immune cells is shown in Figure 7D. Further investigation of the relationship between 22 immune cells in all samples revealed that Macrophages M1 and T cells CD4 memory resting show a strong positive correlation, while B cells naive and Macrophages M2 exhibit a strong negative correlation (Figure 7E). The correlation analysis between gene A and immune cells is shown in Figure 8A. We found that GSTM5 is significantly positively correlated with T cells regulatory (Tregs) (r  = 0.66, p=0.019) (Figure 8B) and T cells follicular helper (r  = 0.59, p=0.041) (Figure 8C).

Figure 7 Immune analysis. (A) Immune functional analysis of the GSTM5 gene; (B) DFU immune infiltration analysis; (C) PCA analysis of immune infiltration results; (D) Immune cell distribution; (E) Immune cell correlation analysis. Notes: *p<0.05.

Figure 8 Immune cell correlation analysis of GSTM5. (A) Correlation analysis; (B) GSTM5 correlation analysis with T cells follicular helper; (C) GSTM5 correlation analysis with T cells regulatory (Tregs).

Gene Expression in DFU

The results revealed that in the control group, the scabs on the wound fell off, the wound size decreased, and there was redness. In the DFU group, the wound surface appeared yellow-brown scabs with a small amount of thin purulent discharge (Figure 9A), and the line graph of wound size measurement is found in Supplementary Material 1. Histological staining with HE (Hematoxylin and Eosin) and Masson’s trichrome revealed that in the control group, there was epidermal regeneration with a small amount of fibroblasts and collagen, and a large number of capillaries in the dermis. In the DFU group, the skin appeared thinner with disordered and reduced epidermal cell arrangement, scattered fibroblasts, fewer newly formed blood vessels, reduced lumens, and extensive infiltration of inflammatory cells (Figure 9B). In conclusion, these results indicate the successful establishment of the DFU model. The IHC results demonstrated that the control group exhibited positive expression of GSTM5, while the DFU group showed significant inhibition (Figure 9C).

Figure 9 Animal experimental validation. (A) Optical pictures of the wound closure rate in the Control group and DFU group at days 0, 7, and 11 after the skin operation; (B) H&E and Masson staining images of wound tissues in the Control group and DFU group at day 7; (C) Representative images of immunohistochemical staining of GSTM5. Note: ***P<0.001.

Discussion

DFU is one of the most common complications of diabetes, causing not only burdens for patients but also posing challenges for overall health, nursing practices, and the social environment.³² Patients with diabetes often experience weakened immune systems, and when coupled with the vulnerability of the feet as one of the body’s vital organs, infections can exacerbate inflammation and ultimately play a role in the onset of DFU.³³ This study, from an immunological perspective, identified GSTM5 as a potential key target for DFU using machine learning and bioinformatics methods. The low expression of GSTM5 in both the validation set and animal experiments among the DFU group suggests the reliability of the findings in this study.

Based on the conclusions of this study, GSTM5 is significantly underexpressed in DFU tissues, suggesting that it may be a potential protective factor in DFU patients. GSTM5 is a member of the Mu class within the glutathione S-transferase (GST) family, and these family members are distributed across a 97 kb region on chromosome 1p13.³⁴ The GST gene family is one of the major xenobiotic detoxification enzymes, protecting cells from the harmful effects of toxic drugs and environmental electrophiles.^35,36 Being a member of the glutathione S-transferase family, GSTM5 is linked to multiple types of cancers, such as breast cancer, prostate cancer, ovarian cancer, and colorectal cancer.³⁷ Earlier research has suggested that GSTM5 expression is reduced in bladder cancer tissues when compared to normal tissues. Furthermore, it has been demonstrated that elevated GSTM5 levels can lead to a reduction in the proliferation, migration, and colony-forming capacity of bladder cancer cells. These findings imply that GSTM5 may effectively serve as a tumor suppressor in bladder cancer, offering a substantial protective effect.³⁸ Furthermore, studies have also indicated that GSTM5 has low expression levels in breast cancer,³⁸ prostate cancer,³⁹ and Barrett’s adenocarcinoma.⁴⁰ GSTM5 encodes a protein belonging to the glutathione S-transferase family, and these proteins are involved in limiting oxidative damage to tissues.⁴¹ Oxidative damage is a key driving factor in the development of DFU.⁴² In the current study, GSTM5 has shown a significant protective effect in most diseases, and our IHC testing also indicates its low expression in DFU tissues.

The inflammatory response, encompassing early reactions following an injury,⁴³ holds a pivotal role in facilitating the healing process. The recruitment of immune cells to the injury site aids in orchestrating multifaceted healing responses among cells and serves to thwart infections.^44,45 The process of acute wound healing comprises several sequential and overlapping stages, including hemostasis, inflammation, proliferation, and remodeling.⁴³ The inflammatory response is often regarded as the “driving force” that initiates the repair process, and because of its capacity for potential harm, it is closely controlled and regulated. The results of immune infiltration indicate that in the DFU group, B cells naive and dendritic cells are highly expressed, while in the control group, T cells CD4 memory resting, T cells regulatory (Tregs), T cells follicular helper, Macrophages M1, and Mast cells resting are highly expressed. The ssGSEA results also suggest that GSTM5 regulates immune functions related to B cells and cytolytic activity, and the immune cells T cells regulatory (Tregs) and T cells follicular helper were found to be positively correlated with GSTM5 in this study. Studies suggest that the timing and effectiveness of wound healing can be deliberately manipulated through the direct application of mature, unactivated B cells. This strategy can successfully modify the equilibrium of mature immune cell populations within the wound microenvironment, thereby expediting the healing process.⁴⁶ Tregs, a distinctive and vital subset of T cells, are identified by their expression of Foxp3 and are classified as CD4+CD25+ cells. They play a unique and essential role in immune surveillance.^47,48 Studies have indicated that T cells migrate to the site of injury during the proliferative phase.⁴⁹ According to studies, dendritic epidermal T cells have been demonstrated to enhance impaired wound healing in diabetic mice. They achieve this by promoting the migration and proliferation of keratinocytes while reducing keratinocyte apoptosis in an environment resembling diabetes.⁵⁰ Additionally, studies have indicated that factors secreted by T cells may help drive macrophage polarization during the wound healing process.⁵¹

The enrichment analysis suggests that GSTM5 may be correlated with NF-κB, GnRH, and MAPK signaling pathway. Elevated expression of NF-κB genes or proteins can result in increased oxidative stress, potentially impacting the wound healing process in diabetic patients.^8,52 Current understanding suggests that in diabetic wounds, the persistent elevation of TNF-α triggers NF-κB activation, which in turn causes the release of numerous inflammatory cytokines, thereby sustaining an ongoing inflammatory response. Consequently, various chemotactic factors and cytokines essential for wound healing lose their effectiveness, resulting in prolonged healing durations.⁵³ As per available information, the use of GnRH agonists is linked to an increased risk of diabetes.⁵⁴ Research suggests that in diabetic models, mice treated with GnRH agonists exhibit more delayed wound healing compared to untreated control mice, and treatment with GnRH antagonists inhibits the damage caused by GnRH agonist treatment.⁵⁵ The activation of the MAPK signaling pathway acts as a critical junction for numerous trauma-induced signaling pathways, facilitating wound healing. Nevertheless, during the wound healing process in diabetes, the local hyperglycemic condition can disrupt normal repair mechanisms. Adequate keratinocyte migration is essential for wound healing, yet high glucose levels hinder keratinocyte migration by suppressing the p38 MAPK/autophagy pathway. This interference results in challenges in wound healing.⁵⁶ Studies have indicated that the hyperglycemic condition in individuals with diabetes activates the MAPK pathway through both protein kinase and non-protein kinase pathways. An upregulation of MAPK is associated with the progression of DFU. Inhibiting the MAPK pathway can offer therapeutic benefits.⁵⁷

Indeed, there is still a relative lack of research on biomarkers in the context of DFU. Recently, the combination of bioinformatics and machine learning analysis has emerged as a valuable tool for discovering novel biomarkers in DFU.⁵⁸ A previous study suggested that MAPK3 may be a biomarker associated with iron death in DFU.⁷ It is worth noting that there is currently a lack of research in the literature concerning the immune system and its relevance to DFU. This study combines machine learning and immune infiltration methods to re-identify GSTM5 while also providing insights into its related functions and mechanisms. Additionally, the low expression of GSTM5 in both the validation set and animal experiments among the DFU group suggests the reliability of the findings in this study. While encouraging results have been found, it is important to acknowledge the limitations of this study. To begin with, the clinical data utilized in this study were sourced from public databases, and the available clinical information for the samples was incomplete. This included limitations such as the unavailability of clinical pathological characteristics for the GSE series. Secondly, the data employed in this study relied on RNA sequencing. While animal experiments corroborated the findings of the bioinformatics analysis, it is imperative to validate the reproducibility and generalizability of these results in future investigations using clinical samples. However, this was not feasible within the constraints of the limited timeframe, as clinical samples were not readily accessible. Lastly, even if diabetic rats were created, the wound model created in these rats is an acute wound model, not a chronic wound.

Conclusion

In conclusion, this study identifies GSTM5 as a potential immune-related key target in DFU using machine learning and bioinformatics methods. Subsequent validation through external datasets and animal experiments also confirms GSTM5 as a critical biomarker for DFU.