Single-cell RNA sequencing reveals a new mechanism of endothelial cell heterogeneity and healing in diabetic foot ulcers

Abstract

Diabetic foot ulcers (DFU) are a common and severe complication among diabetic patients, posing a significant burden on patients’ quality of life and healthcare systems due to their high incidence, amputation rates, and mortality. This study utilized single-cell RNA sequencing technology to deeply analyze the cellular heterogeneity of the skin on the feet ofDFU patients and the transcriptomic characteristics of endothelial cells, aiming to identify key cell populations and genes associated with the healing and progression of DFU. The study found that endothelial cells from DFU patients exhibited significant transcriptomic differences under various conditions, particularly in signaling pathways related to inflammatory responses and angiogenesis. Through trajectory analysis and cell communication research, we revealed the key role of endothelial cell subsets in the development of DFU and identified multiple important gene modules associated with the progression of DFU. Notably, the promoting effect of the SH3BGRL3 gene on endothelial cell proliferation, migration, and angiogenic capabilities under high glucose conditions was experimentally verified, providing a new potential target and theoretical basis for the treatment of DFU. This study not only enhances the understanding of the pathogenesis ofDFU but also provides a scientific basis for the development ofnew therapeutic strategies.

Introduction

Diabetic foot ulcers (DFU) are a common and severe complication affecting millions of people with diabetes worldwide. As the number of individuals with diabetes increases, so does the incidence of DFU, reaching up to 34% in patients with a history of foot injury or infection, and it is projected to affect 693 million people by 2045 [1]. The occurrence of DFU is primarily associated with hyperglycemia, peripheral neuropathy, and vascular disease, which not only impair wound healing but can also lead to amputation and death [2]. Studies have shown that approximately 20% of DFU patients eventually require amputation, and 10% face the risk of death within a year after initial diagnosis, significantly impacting the quality of life and healthcare systems [3]. Current treatment methods for DFU include wound debridement, pressure relief, blood glucose control, and infection prevention, yet the effectiveness of existing therapies remains unsatisfactory [4]. Emerging treatments such as negative pressure wound therapy and hyperbaric oxygen therapy are gradually being introduced into clinical practice but require further research to enhance their efficacy [5]. Therefore, investigating new biomarkers and therapeutic strategies to improve treatment outcomes and prognosis for DFU patients is crucial [6]. Effective management of DFU requires a comprehensive approach, including blood glucose control, improved blood circulation, debridement, pressure relief, and infection control, to reduce amputation rates and improve the quality of life for patients.

Diabetic foot ulcers are a common and severe complication among diabetic patients, significantly affecting their quality of life and imposing a substantial burden on healthcare systems due to their high incidence, amputation rates, and mortality. The healing process of DFU is influenced by various factors, among which the function of endothelial cells (ECs) is particularly crucial. Endothelial cells form the inner lining of blood vessels, are responsible for regulating the exchange of substances between the blood and surrounding tissues, and play a key role in immune and inflammatory responses [7]. Studies have shown that diabetes is characterized by endothelial dysfunction and reduced angiogenesis, especially in the presence of peripheral arterial disease, where the prognosis for DFU patients with foot ulcers is often poor, likely closely related to the biological damage caused to endothelial cells by hyperglycemia and hypoxia [8, 9]. Although the role of endothelial cells in DFU healing is vital, research on the transcriptomic characteristics of endothelial cells in non-healing DFU remains relatively scarce. Notably, endothelial cells from DFU patients exhibit significant differences in transcriptomic features compared to healthy individuals, indicating the need for further exploration of the functional characteristics of DFU-derived endothelial cells. During the healing process of DFU, an increase in cellular ATP levels can promote angiogenesis and collagen synthesis, thereby accelerating wound healing [10]. However, epigenetic changes under high glucose conditions may affect the expression of key factors such as NF-kappaB, thereby disrupting the normal function of endothelial cells [11]. Therefore, a deeper understanding of the role and potential mechanisms of endothelial cells in diabetic foot ulcers is of great significance for advancing the development of new therapeutic strategies and improving the quality of life for patients.

The application of single-cell sequencing technology in the study of diabetic foot ulcers is rapidly evolving, providing new perspectives for understanding its pathogenesis. This technology’s high accuracy and specificity make it an ideal tool for studying cellular heterogeneity and dynamic changes, enabling high-throughput, unbiased analysis even with minimal sample sizes. Through single-cell RNA sequencing (scRNA-seq), researchers can delve into the transcriptomic characteristics of different cell types in the DFU microenvironment, thereby identifying key cell populations associated with healing and non-healing DFU [12, 13]. Additionally, scRNA-seq can help identify novel biomarkers related to DFU treatment, providing potential targets for the development of personalized therapeutic strategies [14]. This study aims to explore the heterogeneity of endothelial cells in DFU and their role in the healing process using single-cell RNA sequencing technology. We downloaded data from 33-foot skin samples, including healthy non-diabetic subjects, diabetic patients without DFU, patients with healed DFU, and patients with non-healed DFU. The study found that endothelial cells from DFU patients exhibit significant transcriptomic differences under various conditions, particularly in signaling pathways related to inflammatory responses and angiogenesis. Furthermore, through trajectory analysis and cell communication research, the study revealed the key role of endothelial cell subsets in the development of DFU and identified multiple important gene modules associated with DFU progression. Ultimately, the study experimentally validated the promotive effect of the SH3BGRL3 gene on endothelial cell proliferation, migration, and angiogenic capabilities under high glucose conditions, providing a new potential target and theoretical basis for the treatment of DFU.

Materials and methods

Sources of raw data

We downloaded single-cell RNA sequencing data from 33-foot skin samples in the GSE165816 dataset, which includes 11 cases from healthy non-diabetic subjects, 8 cases from diabetic subjects without DFU, 9 cases from subjects with healed DFU, and 5 cases from subjects with non-healing DFU [15]. To screen for genes associated with DFU, the GSE80178 cohort, which includes 6 DFU patients and 3 diabetic patients without DFU, was utilized for machine learning analysis [16]. This dataset is a gene expression array generated by GPL16686 (Affymetrix Human Gene 2.0 ST Array). Preprocessing of the raw chip data was conducted in R, including background correction, normalization, and log2 transformation [17].

Raw processing of single-cell RNA sequencing

Using the “Seurat” R software package (version 4.2.0), the raw gene expression matrix was imported into R software [18]. The “DoubletFinder” R software package was employed to remove doublets [19]. The “DecontX” R software package was utilized to eliminate environmental-free RNA contamination from single-cell transcriptomes [20]. A Seurat object was generated using the “Seurat” R package, and cells with gene expression levels ranging from 300 to 7000, mitochondrial genes below 20%, and hemoglobin genes below 10% were filtered out. The “FindVariableFeatures” function was used to identify the top 2000 highly variable genes. Data normalization was performed using the “ScaleData” function. Subsequently, the remaining single-cell transcriptome expression matrices were integrated using the “harmony” R software package [21]. High-variance genes were selected for principal component analysis (PCA), and the first 30 significant principal components were used for Uniform Manifold Approximation and Projection (UMAP) dimensionality reduction. The “FindAllMarkers” function was employed to determine differentially expressed genes (DEGs) in each cellular subset, and cell types and subtypes were annotated based on the expression of established canonical marker genes for each cell type. In this study, we employed the “miloR” R software tool for gene set enrichment analysis of single-cell data. MiloR is a k-nearest neighbor (KNN) graph-based differential abundance analysis method that can identify differentially expressed gene sets from single-cell datasets [22].

Trajectory analysis and inference of cellular stemness

In this study, we utilized the “CytoTRACE” R package to infer the cellular differentiation states from single-cell RNA sequencing data [23]. This package predicts the differentiation potential of endothelial cells by analyzing gene count signatures, independently of specific time scales or continuous developmental processes. We employed Monocle3 for the analysis of single-cell data, beginning with the dimensionality reduction and clustering of the Seurat object, followed by its conversion into a Monocle object to reconstruct the differentiation trajectories [24]. Based on the results from CytoTRACE and selecting the starting point from normal skin cell populations, we traced the pseudo-time trajectories of endothelial cells from the healthy group toward the DFU group. Subsequently, the Monocle2 algorithm was applied to analyze the developmental trajectories of two inferred endothelial cell subsets. After dimensionality reduction and unit ordering, cell trajectories were inferred using default parameters [25].

Pathway enrichment analysis and prediction of transcription factors

To explore the biological significance of DEGs, we utilized the “clusterProfiler” R package for Gene Ontology (GO) and KEGG pathway enrichment analysis [26]. Metabolic activities were assessed using the “scMetabolism” R package [27]. “irGSEA” integrates multiple enrichment analysis methods based on individual cell expression levels, including AUCell, UCell, singscore, ssGSEA, JASMINE, and Viper, and evaluates the results of differential analysis using the rank aggregation algorithm (robust rank aggregation, RRA) to screen for differential gene sets that show similar enrichment levels across these methods [28]. We employed the SCENIC (Single-Cell ENrichment of INteractive gene Clusters) tool to infer gene regulatory networks (GRNs). SCENIC constructs and infers gene regulatory networks for cell types by analyzing single-cell RNA-seq data in conjunction with transcription factor (TF) motifs and gene expression data [29].

Cellular communication

In this study, we employed the MultiNicheNet framework to analyze cell–cell communication within single-cell transcriptomic data across multiple samples and conditions. The primary objective of MultiNicheNet is to infer differentially expressed ligand-receptor pairs under various conditions and predict the potential downstream target genes for these pairs. MultiNicheNet is based on the NicheNet framework and utilizes the same prior knowledge network. We can visualize ligand activity for specific group-receiver combinations and display the predicted ligand-target links as well as the expression of predicted target genes across different samples [30, 31].

HdWGCNA analysis

High-density weighted gene co-expression network analysis (hdWGCNA) was used to identify key genes associated with DFU endothelial cells in DFU samples [32]. Endothelial cell subsets were filtered from the scRNA data, and gene expression correlation matrices, weighted gene co-expression networks, and module detection were performed. Module-trait relationship analysis identified modules significantly associated with DFU and determined hub genes within important modules based on connectivity within the modules. The top 15 hub genes were considered to be key genes related to the progression of DFU. For the GSE80178 cohort, the least absolute shrinkage and selection operator (LASSO) regularization regression algorithm in the “glmnet” R package was used on the 30 key genes to determine genes associated with the progression of DFU.

Cell culture and cytological experiments

Human umbilical vein endothelial cells (HUVECs) purchased from the American Type Culture Collection (ATCC) (Manassas, USA), which are immortalized, were cultured in Dulbecco’s Modified Eagle Medium (DMEM; Thermo Fisher Scientific, Waltham, MA) supplemented with 10% Fetal Bovine Serum (FBS, Gibco) at 37 °C and 5% CO2. DMEM containing 5.5 mM glucose was considered the normal glucose (CON) condition, while the high glucose (HG) condition contained 30 mM glucose, achieved by the addition of extra glucose. In this study, lentiviral-mediated gene overexpression technology was used to achieve stable expression of the target gene in human umbilical vein endothelial cells under high-glucose conditions. A lentiviral packaging kit (Genechem, Shanghai, China) was used to transfect cells according to the manufacturer’s instructions. Briefly, the Lenti-Easy Packaging Mix and the target gene lentiviral plasmid were co-transfected into 293 T packaging cells. After filtering the viral supernatant, the target cells (HUVECs) were infected for subsequent experiments. The SH3BGRL3 overexpression plasmid sequence is as follows: Forward Primer (5’—3’): TTTCGACATTTAAATTTAATATGGTCATCCGCGTGTTCAT; Reverse Primer (5’—3’): ATTCCTGCAGCCCGTAGTTTCTAAGGTTCTGCCTTTGATG. Normal glucose and high-glucose control groups were established as controls, and the role of the target gene in diabetic wound mechanisms was further assessed in combination with functional experiments (such as angiogenesis assays, scratch assays, and flow cytometry).

5-Ethynyl-2’-deoxyuridine (EdU) Assay: Cells were cultured in a 24-well plate at 60% confluence and incubated with 10 μM EdU solution. After fixation with 4% paraformaldehyde and permeabilization with 0.3% Triton, cells were stained with the click reaction mixture and DAPI. Representative images were captured using an inverted fluorescence microscope.

Scratch Assay: Cells were seeded at a density of approximately 5 × 10^5 cells/well and grown to confluence. Two perpendicular scratches were made with a sterile pipette tip, and photographs were taken at different time points to assess scratch healing. Statistical analysis was performed using ImageJ software.

Transwell Assay: 5 × 10^4 cells/well were added to the upper chamber of a transwell insert with 8 μm pores, while 600 μL of 10% serum-containing medium was added to the lower chamber. After 24 h, non-migrated cells were removed, and the migrated cells were fixed with 4% paraformaldehyde and stained with crystal violet. The number of migrated cells was counted under a light microscope.

Angiogenesis Assay: Matrigel matrix was spread evenly into a 96-well plate and solidified at 37 °C for 30 min. Cell suspensions (3000 cells/well) were added and incubated for 6 h. Tubule formation was observed and photographed using an inverted microscope, and analyzed using ImageJ software.

Western blot

To assess protein expression in HUVECs following different treatments, we performed Western blot analysis. Total protein from each group was extracted using RIPA lysis buffer containing protease and phosphatase inhibitors. Protein concentration was quantified using a BCA protein assay kit. Equal amounts of protein (20–30 µg) were loaded onto SDS-PAGE gels and separated by electrophoresis, followed by transfer to PVDF membranes. The membranes were blocked at room temperature with TBS-T containing 5% non-fat dry milk for 1 h and then incubated overnight at 4 °C with primary antibodies specific to the proteins of interest. After washing with TBS-T, the membranes were incubated with HRP-conjugated secondary antibodies at room temperature for 1 h. Protein bands were visualized using an enhanced chemiluminescence (ECL) detection system and quantified using ImageJ software. GAPDH was used as a loading control for normalization.

Mouse wound healing experiment

To investigate the role of SH3BGRL3 in diabetic wound healing, we used a well-established mouse model of type 1 diabetes induced by streptozotocin (STZ). A total of 20 healthy 10-week-old mice were randomly assigned into two groups (n = 10 per group): the AAV-SH3BGRL3 overexpression group and the AAV-NC control group. The mice were selected based on their age and health status to ensure consistency across the study. Adeno-associated viruses (AAV-SH3BGRL3 and AAV-NC) were purchased from Genechem (Shanghai, China) and administered via a single local subcutaneous injection on the back skin 21 days before wounding, with a total viral dose of 1 × 10^12 viral genomes (VG) per mouse (50 μl). To induce diabetes, the mice were fasted overnight and then intraperitoneally injected with streptozotocin (STZ, 50 mg/kg, Sigma, Missouri, USA) for five consecutive days. Blood glucose levels were continuously monitored one week later using a glucometer to confirm the induction of diabetes. On the day of wounding, mice were anesthetized via intraperitoneal injection of sodium pentobarbital (50 mg/kg). A full-thickness circular wound with a diameter of approximately 1 cm was created on both sides of the mouse’s back using a circular mold. The wound healing process was observed and photographed on days 0, 3, 7, 10, and 14. The wound area was measured using ImageJ software, and the data were statistically analyzed.

Tissue immunofluorescence

Paraffin sections were placed in an oven at 60 °C for 1 h to melt the paraffin. Then, the sections were deparaffinized with two 15-min washes in xylene, followed by 5-min washes in 100%, 95%, 85%, and 75% ethanol. Finally, the sections were rinsed three times for 5 min each with PBS. The sections were immersed in a citrate antigen retrieval solution (Beyotime, Shanghai, China) in a plastic staining jar and microwaved at high power (p-100) for 5 min to boil; and then at medium power (p-60) for 5 min. After cooling to room temperature, the sections were washed twice with tap water for 5 min each. A hydrophobic barrier was drawn around the tissue with an immunohistochemistry pen (Beyotime, Shanghai, China). The sections were then blocked with bovine serum albumin for 30 min. After shaking off the excess and without washing, the primary antibody was applied and incubated overnight at 4 °C for 18 h. The next day, the sections were left at room temperature for 1 h before being washed with PBS for 5 min three times. The corresponding secondary antibody was applied and left to stand at 37 °C for 1 h. After three 5-min washes with PBS, the sections were mounted with DAPI. Photographs were taken using an inverted fluorescence microscope, and the data were analyzed using ImageJ. Antibodies used were: Ki67 (1:500 dilution, Cat No. HA721115, HUABIO); TNFA (1:50 dilution, Cat No. HY-P80914, MCE); Collagen I (1:200 dilution, Cat No. ab34710, abcam); Collagen III (1:200 dilution, Cat No. ab7778, Abcam).

Statistical analysis

Statistical analysis was conducted using R software version 4.2.1. This included Student’s t-tests or Wilcoxon tests, with results expressed as the mean ± standard deviation. The false discovery rate (FDR) method was employed to adjust p-values. The “ClusterGVis” and “scRNAtoolVis” R packages were utilized to enhance the heatmaps and volcano plots of this study (available at https://github.com/junjunlab/ClusterGVis and https://github.com/junjunlab/scRNAtoolVis). The “SCP” R package, designed as an end-to-end solution for single-cell data analysis, assisted in the bioinformatics visualization of this study (available at https://github.com/zhanghao-njmu/SCP). All p-values were calculated using a two-tailed approach, with statistical significance defined as p < 0.05.

Results

Single-cell landscape of foot skin

Figure 1A succinctly outlines the bioinformatics analysis and biological experimental workflow of this study. We analyzed a publicly available single-cell dataset (GSE165816) to elucidate the heterogeneity of various cell types in diabetic foot ulcers (DFU) and to explore the pathogenesis of inflammatory progression. We conducted a detailed analysis of the single-cell RNA sequencing (scRNA-seq) results from the foot skin of normal individuals without diabetes, diabetic patients without skin ulcers, and patients with either healing (Healing DFU) or non-healing (Non-Healing DFU) foot ulcers. After stringent quality control measures, including the removal of doublets and correction of RNA sequencing contamination, and subsequent data normalization and unsupervised clustering, we obtained the transcriptomes of 46,297 cells (healthy individuals: 29,861; diabetic patients: 10,371; Healing DFU: 23,358; Non-Healing DFU: 15,184). Unsupervised Uniform Manifold Approximation and Projection (UMAP) clustering revealed 27 cell clusters with comparable expression patterns (Supplementary Fig. 1A, B). Drawing on several previous studies [15, 33, 34], we utilized canonical markers to verify the distinct cell populations within each cluster: fibroblasts (expressing FGF7, DCN, and CFD), endothelial cells (expressing PECAM1, CDH5, and CLDN5), smooth muscle cells (expressing TAGLN, ACTA2, and TPM2), keratinocytes (expressing KRT5, KRT14, KRT10, and KRT1), lymphatic endothelial cells (expressing CCL21 and LYVE1), sweat gland cells (expressing DCD and AQP5), Schwann cells (expressing S100B and CDH19), melanocytes (expressing MLANA, PMEL, and TYRP1), T cells (expressing CD3E and CD3D), NK cells (expressing GNLY, NKG7, and CCL5), macrophages (expressing CSF1R, CD68, and CD14), dendritic cells (expressing CD1C and CD1A), mast cells (expressing TPSAB1 and MS4A2), B cells (expressing MS4A1 and CD79A), and plasma cells (expressing MZB1 and IGHG1) (Fig. 1B, Supplementary Fig. 1C, D). The cellular lineage profiles from different samples and groups exhibited varying relative cell ratios (Fig. 1C, D). We observed an increased proportion of macrophages and B cells, along with downregulation of endothelial and Schwann cells in the Healing DFU and Non-Healing DFU groups, which appears to be associated with cellular apoptosis and inflammatory lesions in diabetic neuropathy [35]. Heatmaps displayed the principal marker genes of 13 major cell subtypes (Fig. 1E). UMAP plots illustrated the expression of specific marker genes across different cell types (Fig. 1F).

Fig. 1

Single-cell atlas of foot skin from individuals without diabetes, diabetic patients without skin ulcers, Healing DFU, and Non-Healing DFU. A An overview of the study design and workflow. B UMAP plot of single-cell features colored by 13 major cell types identified in this study. C, D Bar graphs showing the relative proportions of each major cell type in each sample and group. E Heatmap displaying the expression levels of major marker genes across all main cell types. Red indicates upregulated gene expression in that cell type. F UMAP plots featuring typical marker genes of each cell type

Heterogeneity and distribution characteristics of vascular endothelial cells in foot skin.

Vascular endothelial cells play a crucial role in the development and healing of diabetic foot ulcers by participating in processes such as angiogenesis, inflammatory responses, endothelial dysfunction, oxidative stress reactions, and tissue regeneration [36]. We performed a re-clustering of endothelial cells and identified 11 major clusters across four groups using UMAP (Fig. 2A, B). A recent study classified vascular endothelial cells in the skin into venous/lymphatic, capillary, and arterial/arteriolar states, as well as inflammatory conditions; we collected all relevant markers and annotated the 11 endothelial cell subsets [34, 37]. We found that EC-7 and EC-11 exhibit capillary states, EC-3 is more likely to be arterial/arteriolar, and EC-1 and EC-10 have distinct inflammatory characteristics. EC-11 also shows a transcriptional signature of mesenchymal activation (Fig. 2E). Notably, compared to the normal and diabetic groups, EC-2, EC-5, and EC-8 have higher expression proportions in the DFU group, while EC-1, EC-3, EC-7, and EC-10 have lower expression proportions in the DFU group. Interestingly, EC-9 is specifically present in diabetic patients but shows no expression in the normal and DFU groups (Fig. 2C, D). To further explain the differential abundance of endothelial cell subsets in different disease states, we used moliR to maintain FDR error control in batch effects. The study results show that, whether comparing the diabetic group with the Healing DFU group or the healthy group with the Non-Healing DFU group, EC-5 and EC-8 both exhibit significantly higher cell abundance in DFU (both healing and non-healing). Specifically, EC-5 tends to be enriched in the Non-Healing DFU group, while EC-8 is more enriched in the Healing DFU group. In contrast, EC-7 has higher cell abundance in the healthy and diabetic groups compared to DFU (Fig. 2F). Finally, we presented the specifically highly expressed and lowly expressed genes of the 11 endothelial cell subsets in the form of volcano plots (Fig. 2G, H).

Fig. 2

Heterogeneity of vascular endothelial cells. A, B UMAP plots of 11 vascular endothelial cell clusters obtained from further clustering analysis. C, D Bar graphs showing the relative proportions of each endothelial cell subtype in each sample. E Dot plot displaying the expression of vascular marker genes across 11 endothelial cell subtypes. F Comparison of endothelial cell subtype abundance between different groups. G, H Volcano plots show all upregulated and downregulated genes in the 11 endothelial cell subtypes

Signal enrichment and transcriptional characteristics of endothelial cell subsets

We further explored the biological characteristics and enriched signaling pathways of these endothelial cells with distinct transcriptional profiles. By conducting enrichment analysis on the marker genes of the 11 endothelial cell subsets, we found that EC-5 is primarily associated with the regulation of acute inflammatory responses in GO_BP and with the TNF signaling pathway in KEGG. EC-8, which has healing characteristics, is mainly associated with the AMPK and PPAR signaling pathways. EC-9, which is specifically expressed in the diabetic group, is related to extracellular matrix organization and the PI3K − Akt signaling pathway (Fig. 3A). HALLMARK enrichment analysis showed that both EC-5 and EC-8 are positively correlated with epithelial-mesenchymal transition. EC-5, which has non-healing characteristics, mainly exhibits inflammatory responses and activation of the TNFA signaling pathway, while EC-8, which has healing characteristics, is associated with pathways such as oxidative phosphorylation, hypoxia, and the PI3K/AKT/MTOR signaling pathways (Fig. 3B). Further transcription factor analysis of the subsets revealed that EC-5 is significantly associated with the expression and regulatory activity of ZNF281, while EC-8 exhibits stronger ATF4 transcription factor activity (Fig. 3C). Given the different expression patterns of these endothelial cell subsets in metabolism, we compared their activity scores across all metabolic pathways in detail and assessed the metabolic scores of each cell using the AUCell method. EC-5 showed relatively high activity in phosphate metabolism, while EC-8 did not have specific metabolic characteristics. EC-9, which specifically appears in the diabetic group, exhibited very high fructose metabolism characteristics, and EC-11 showed high activity in all metabolic pathways (Fig. 3D).

Fig. 3

Enrichment analysis of vascular endothelial cells. A GO-BP/KEGG enrichment analysis results for differential genes and their corresponding transcription factors (TFs) and cell surface proteins (CSPA) in 11 vascular endothelial cell subtypes. B HALLMARK enrichment analysis of 11 vascular endothelial cell subtypes. C SCENIC analysis of transcription factor expression and activity in 11 vascular endothelial cell subtypes. D Dot plot of metabolic pathway activity across 11 vascular endothelial cell subtypes

Signaling of endothelial cell subsets’ differentiation

We further explored the biological characteristics and enriched signaling pathways of these endothelial cells with distinct transcriptional profiles using CytoTRACE to assess the differentiation potential of different endothelial cell subsets. Interestingly, EC-3 and EC-4 exhibited the highest stemness features, while EC-2, EC-5, EC-8, and EC-9 had the lowest stemness (Supplementary Fig. 2A, B). Based on the previous moliR analysis, we found that EC-5 and EC-8 were more abundant in the DFU group, so we chose EC-5 and EC-8 as the endpoints for cell differentiation and conducted Monocle3 analysis. The pseudo-time trajectory revealed a developmental starting state located at EC-3 and another at EC-10, indicating two directions of differentiation for arterial endothelial cells and venous endothelial cells (Supplementary Fig. 2C). Monocle2 was used to determine the cell trajectories and pseudotime distributions of endothelial cells. We observed a total of nine cell states during the development process. EC-5, which has non-healing characteristics, mainly appeared in the early stages of cell development, while EC-8, which has healing characteristics, was mainly distributed in the later stages of cell development (Supplementary Fig. 2D). Along with the process of cell development, the proportion of endothelial cells in the Healing DFU group increased (Supplementary Fig. 2E). In the pseudo-time study, we identified the top 50 highly variable genes (Supplementary Fig. 2F). We conducted GO enrichment analysis in different cellular differentiation states, and the results showed that the late stage of endothelial cell differentiation was mainly associated with biological processes such as cytoplasmic translation (Supplementary Fig. 2G).

Endothelial cell-related signaling communication and downstream signal transduction

Intercellular dialogue forms the fundamental means by which cells communicate and cooperate within multicellular organisms, often prompting cells to change their status and function. Compared to other cell communication workflows, we selected the most appropriate MultiNicheNet method to handle the complex skin samples. We investigated the ligand-receptor binding activity and target gene expression that may be involved in these processes by assessing the top cell–cell interactions with EC-5 and EC-8 as recipients, where EC-5 and EC-8 play different roles in cell signal transduction. As shown in Supplementary Fig. 3A, the highest priority ligands between fibroblasts and EC-5, which has a non-healing phenotype, include SLIT3, APOE, SERPINF1, and JAG1, genes that are mainly associated with inflammatory responses and cell death processes. The signals between EC-8, which has a healing phenotype, and fibroblasts mainly include CXCL12-CXCR4, COL14A1, and COL6A6, which are related to cell adhesion, migration, and cell cycle progression. Further investigation of the mechanisms of these genes and signaling pathways in different cellular phenotypes may help develop new therapeutic strategies to promote tissue healing or inhibit pathological angiogenesis and inflammatory responses. Subsequently, we further predicted the specific downstream affected target genes related to the ligand-receptor links associated with EC-8 (Supplementary Fig. 3B). Not surprisingly, CXCL12 in EC-8 specifically upregulated the expression of downstream genes such as KRT14, KRT17, KRT6A, MT2A, and CXCR4, and these regulatory networks play a crucial role in endothelial cell adhesion, translocation, and vascular reconstruction in ischemic tissue growth.

Identification of key gene modules associated with DFU Using high-dimensional weighted gene co-expression network analysis

We employed HdWGCNA to identify the key molecular signatures most closely associated with the progression of diabetic foot ulcers (DFU). During the construction of the co-expression network, we observed that the scale-free topology fit index reached 0.9 at a soft-thresholding power of 8 (Fig. 4A). This analysis led to the identification of 14 gene modules (Fig. 4B). We calculated module connectivity to determine the connectivity based on the characteristic genes (Fig. 4C). The correlations between the 11 modules are depicted in Fig. 4D. Subsequently, we assessed the module scores for the 11 endothelial cell subsets and found that the blue and green-yellow modules were highly activated in EC-5 and EC-8, respectively (Fig. 4E, G). We extracted the top 15 significant genes from both modules, which were considered key genes in the progression of DFU. Enrichment analysis results indicated that these genes are associated with protein folding and VEGF signaling pathways (Fig. 4F).

Fig. 4

Differentiation trajectories of vascular endothelial cell subtypes. A, B CytoTRACE analysis to assess stemness in 11 endothelial cell subtypes, with higher CytoTRACE scores indicating increased stemness; statistical significance was assessed using the Kruskal–Wallis test. CytoTRACE scores were mapped onto individual cells’ UMAP to more intuitively represent changes in stemness across different cell subtypes. C Monocle3 pseudotime analysis showing the differentiation paths of endothelial cells from initial to terminal states. D Monocle2 pseudotime analysis displaying the cellular development trajectories and states of two endothelial cell subtypes. E Monocle2 pseudotime analysis showing the proportions of different groups during endothelial cell development. F Heatmap displaying the expression of differential genes in pseudotime analysis. G Heatmap showing GO enrichment analysis results and marker genes for each state of endothelial cell differentiation

SH3BGRL3 and DFU healing association

To further screen for key genes affecting the progression of DFU, we downloaded the GSE80178 cohort, which includes 6 DFU patients and 3 diabetic patients without foot ulcers. By performing GSVA scoring on 30 module genes, patients in the DFU group had significantly higher scores (Fig. 5A). Using LASSO regression analysis, we ultimately identified 2 genes (SH3BGRL3 and TMSB10) that characterize DFU (Fig. 5B, C). UMAP displayed the expression of SH3BGRL3 and TMSB10 in endothelial cells at single-cell resolution (Fig. 5D). Interestingly, SH3BGRL3 expressed the highest in the Healing DFU group, while TMSB10 showed no significant difference in expression between the Healing DFU group and the Non-Healing DFU group (Fig. 5E). This suggests that SH3BGRL3 (SH3 domain-binding glutamic acid-rich-like protein 3) may play a crucial role in the healing process of DFU. Endothelial cells are crucial for the formation of new blood vessels and participate in the regulation of angiogenesis during wound healing [38]. To test the impact of SH3BGRL3 under high glucose (HG) conditions on HUVECs, we employed lentiviral overexpression. It is well known that the V-akt murine thymoma viral oncogene homolog 1 (AKT) signaling, extracellular signal-regulated kinase (ERK) signaling, VEGF-A, and HIF-1α work together to promote the wound healing process by stimulating the proliferation, migration of endothelial cells, increasing vascular permeability, and regulating the expression of genes associated with wound healing. The expression of VEGF-A and HIF-1α, as well as the phosphorylation of AKT and ERK, are significantly inhibited under high glucose (HG) conditions. However, the overexpression of SH3BGRL3 can significantly mitigate the downregulation of these signaling pathways (Fig. 5F, G).

Fig. 5

LASSO regression analysis identifies key genes in the progression of DFU, including SH3BGRL3. A Comparison of GSVA scores for 50 candidate hub genes between the DFU group and the control group. B Changes in the coefficients of selected features using the LASSO algorithm. C Feature selection using the LASSO algorithm. D UMAP plots of SH3BGRL3 and TMSB10 at the single-cell sequencing level. E Expression of SH3BGRL3 and TMSB10 in different groups from single-cell sequencing. F, G Western blotting was used to detect the expression levels of SH3BGRL3, VEGF-A, HIF-1α, ERK, and AKT in HUVECs under high-glucose (HG) conditions and after overexpression of SH3BGRL3. Data are presented as mean values ± standard deviation, n = 3, ** p < 0.01, # p < 0.05, ## p < 0.01, ### p < 0.001

SH3BGRL3 promotes angiogenesis and migration of HUVECs under high glucose conditions

The improvement in cell migration suggests accelerated wound healing. To assess the migratory capacity of HUVECs, we conducted various functional assays to evaluate changes in endothelial function under HG conditions with SH3BGRL3 overexpression. A significant migration delay was observed under HG conditions, and SH3BGRL3 overexpression enhanced the migratory ability of endothelial cells (Fig. 6A, B). EdU assays indicated that SH3BGRL3 restored the proliferative capacity of endothelial cells in a high-glucose environment (Fig. 6C, D). We then stained HUVECs with crystal violet, and the transwell migration assay showed an increase in cell migration ability after SH3BGRL3 overexpression treatment (Fig. 6E, F). Additionally, SH3BGRL3 increased the formation of total tube length, which was significantly reduced under HG conditions (Fig. 6G, H). These findings suggest that SH3BGRL3 plays a significant role in promoting angiogenesis and migration of HUVECs under high glucose conditions, potentially contributing to improved wound healing processes.

Fig. 6

SH3BGRL3 promotes angiogenesis, proliferation, and migration of HUVECs under high-glucose conditions. A, B Wound healing assays indicating that overexpression of SH3BGRL3 upregulated the motility of HUVECs. C, D EdU incorporation assays to evaluate DNA synthesis in HUVECs. Red fluorescence represents EdU-positive cells, and blue fluorescence represents total cells. Representative images and quantitative analysis of the EdU assay. E, F Representative images and quantitative analysis of HUVEC transwell migration assays. G, H Tube formation assays and quantitative analysis of HUVECs (branch points and tube length). HG = High-glucose conditions; NC = Normal control. Data are presented as mean values ± SD, n = 3, *** P < 0.001

SH3BGRL3 accelerates wound healing in diabetic mice

To further assess the role of SH3BGRL3 in the diabetic wound healing process, we established type I diabetic mice by intraperitoneal injection of STZ and created two full-thickness skin wounds on the back of each mouse, observing the healing process for a total of 14 days (Fig. 7A, B). We overexpressed SH3BGRL3 using adeno-associated viruses to evaluate its function in the context of diabetic wound healing. The wounds in diabetic mice treated with the adeno-associated virus showed a significant improvement in delayed healing, as evidenced by the reduction in wound area and residual wounds from day 7 to day 14 post-injury (Fig. 7C-E). H&E staining and Masson’s trichrome staining analysis confirmed that SH3BGRL3 has a more pronounced advantage in promoting re-epithelialization, granulation tissue formation, and improved collagen deposition in wound healing (Fig. 7F-I). Immunofluorescence results indicated that SH3BGRL3 promoted proliferation at the wound site (Fig. 8A, F), enhanced angiogenic capacity (Fig. 8B, J), reduced inflammation levels (Fig. 8C, H), and promoted the deposition of type I collagen (Fig. 8D, I) and type III collagen (Fig. 8E, J), while also reducing the ratio of type I to type III collagen (Fig. 8K). These findings suggest that SH3BGRL3 plays a significant role in accelerating wound healing in diabetic mice.

Fig. 7

SH3BGRL3 accelerates wound healing in diabetic mice. A, B Schematic diagram of the skin wound model in diabetic mice and the timeline for wound photography collection. C Representative images of wound areas at days 0, 3, 7, 10, and 14 post-surgery for different treatments. D Quantitative analysis of remaining wound areas. E Line graph of wound healing over time. F, G HE staining images of wound samples on day 14, and quantitative analysis of wound area length. H, I Masson staining images of wound samples on day 14, and quantitative analysis of collagen area fraction in wounds. HG = High-glucose conditions; NC = Normal control. Data are presented as mean values ± SD, n = 3, *** P < 0.001

Fig. 8

Immunofluorescence of wound tissue. A, F Representative images and quantitative analysis of Ki67 on day 14. B, G Representative images and quantitative analysis of α-SMA on day 14. C, H Representative images and quantitative analysis of TNFA on day 14. D, E, I, J Representative images and quantitative analysis of Collagen I and Collagen III on day 14. (K) The ratio of Collagen I to Collagen III. HG = High-glucose conditions; NC = Normal control. Data are presented as mean values ± SD, n = 3, *** P < 0.001

Discussion

Diabetic foot ulcers (DFUs) are a common and severe complication in patients with diabetes, with complex pathogenesis involving multiple factors. The normal wound healing process includes inflammation, angiogenesis, and extracellular matrix (ECM) remodeling, with cells such as vascular endothelial cells, fibroblasts, keratinocytes, mononuclear macrophages, and lymphocytes playing a role in healing [39]. Cytokines such as Transforming Growth Factor (TGF)-β1, Vascular Endothelial Growth Factor (VEGF), Soluble Vascular Cell Adhesion Molecule-1 (VCAM-1), Platelet-Derived Growth Factor (PDGF), and Epidermal Growth Factor (EGF) play crucial roles in wound healing, with angiogenesis being essential and endothelial cells being key participants [40, 41]. Under normal conditions, endothelial cells line the inner surface of blood vessels and control vasoconstriction and vasodilation by regulating the levels of vasomotor factors, such as endothelial Nitric Oxide Synthase (eNOS) [42]. During wound healing, angiogenesis at different stages is primarily regulated by VEGF. In the inflammatory phase, VEGF increases vascular permeability and promotes the migration of leukocytes to the site of injury; in the proliferative phase, VEGF significantly stimulates the proliferation and migration of endothelial cells; and in the remodeling phase, it promotes the assembly of endothelial cells to form vascular lumens [38]. However, in a diabetic environment, hyperglycemia can impair endothelial cell function. Studies have shown that arterial endothelial cell integrity is reduced under hyperglycemic conditions, making them more prone to apoptosis and shedding, thus affecting angiogenesis [43]. Additionally, reduced levels of Nitric Oxide Synthase (NOS) due to peripheral neuropathy and peripheral arterial disease can cause peripheral vasoconstriction and reduced blood flow. The lack of endothelial progenitor cells (EPCs) at the wound site further inhibits new blood vessel formation, delaying wound healing [44]. Therefore, these changes in endothelial cell biology play a key role in the pathogenesis of DFUs, and a deeper understanding of their mechanisms could aid in the development of new treatment strategies for DFUs.

This study conducted an in-depth analysis of single-cell RNA sequencing data from the foot skin of patients with diabetic foot ulcers (DFU), revealing the critical role of endothelial cells in the occurrence and healing of DFU. Our findings indicate that the foot skin of DFU patients exhibits significant cellular heterogeneity, particularly an increase in the proportion of macrophages and B cells, while the proportion of endothelial cells and Schwann cells significantly decreases. These changes in cellular composition may be closely related to the exacerbation of cellular apoptosis and inflammatory responses in diabetic neuropathy, providing a new perspective for understanding the pathogenesis of DFU. Increasing evidence from single-cell transcriptomic analysis suggests that endothelial cell function exhibits tissue specificity across different sample sources of skin. Although a previous study detailed the heterogeneity of dermal fibroblasts and immune cells in DFU, the characterization of vascular endothelial cells under diabetic conditions is limited [45]. After re-clustering endothelial cells, we identified 11 major subsets and observed significant differences in their expression patterns in the DFU group. Notably, EC-5 and EC-8 showed higher abundance in the DFU group, with EC-5 tending to be enriched in non-healing DFU groups and EC-8 primarily in healing DFU groups. This result suggests that endothelial cell subsets may play different functions in the wound-healing process, providing potential targets for further research.

Further signal enrichment analysis revealed that EC-8 is associated with metabolic pathways such as the AMPK and PPAR signaling pathways. AMPK (AMP-activated protein kinase) plays a crucial role in cellular energy homeostasis, promoting catabolic processes that generate energy and inhibiting anabolic processes that consume energy. In wound healing, the activation of AMPK can enhance the migration and proliferation of endothelial cells, thereby promoting angiogenesis and tissue repair [46]. The PPAR (Peroxisome Proliferator-Activated Receptor) signaling pathway is involved in lipid metabolism and inflammation regulation, and its activation can reduce inflammation and oxidative stress, making it an important factor in promoting healing [47]. Therefore, the association of EC-8 with these signaling pathways suggests that it may promote wound healing by enhancing energy metabolism and reducing inflammation. HALLMARK analysis showed that EC-8 is related to oxidative phosphorylation and the PI3K/AKT/MTOR signaling pathway. Oxidative phosphorylation is an important process in cellular energy metabolism, and the PI3K/AKT/MTOR signaling pathway plays a key role in cell proliferation, survival, and metabolism. The activation of these signaling pathways may promote the positive role of EC-8 in wound healing by enhancing cellular energy production and promoting cell proliferation [48].

In addition, EC-5 exhibits non-healing characteristics and is mainly associated with the inflammatory response and the TNFA signaling pathway. The TNFA (Tumor Necrosis Factor Alpha) signaling pathway is a key mediator of inflammation and immune responses. In diabetic foot ulcers (DFU), chronic inflammation is a major obstacle to healing, leading to tissue damage and impaired angiogenesis. In the healing process of diabetic foot (DFU), inflammation and immunity are closely related. First, the inflammatory response is a normal stage of wound healing, but it is abnormal in DFU. In the inflammatory phase of normal wound healing, inflammatory cells aggregate and secrete inflammatory factors, while in DFU patients, immune cell function is impaired [49]. This indicates that in the DFU healing process, endothelial cells are not only involved in the inflammatory response but may also promote healing by regulating metabolic activities. Vascular endothelial cells (ECs) modulate inflammation by regulating the transport, activation status, and function of immune cells. They also have tissue-specific and vessel-type-specific immune regulatory effects and are known as “immunomodulatory ECs” [7]. Vascular endothelial cells play an important role in immune regulation and wound healing. They are not only involved in the regulation of the inflammatory response but also directly affect the angiogenesis and healing process of the wound. Proper immune regulation is crucial for wound healing, and its dysfunction can lead to delayed wound healing [50]. Interestingly, both EC-5 and EC-8 are positively correlated with epithelial-mesenchymal transition (EMT). EMT is the process by which epithelial cells undergo a transition to a mesenchymal cell phenotype, endowing cells with enhanced migratory capacity and promoting their spread and relocalization within tissues. In wound healing, EMT is a key step in tissue repair. However, in diabetic foot ulcers (DFU), EMT may be dysregulated, leading to chronic inflammation and healing disorders [51]. Additionally, EC-9, specific to the diabetic group, is associated with extracellular matrix organization and the PI3K − Akt signaling pathway, which may imply a mechanism of endothelial cell dysfunction in a diabetic environment. EC-11 shows high activity in all metabolic pathways, indicating that it may play a broad role in cellular metabolism and function and may contribute to multiple aspects of wound healing. Further transcription factor analysis showed that EC-5 is significantly associated with the expression and regulatory activity of ZNF281, while EC-8 exhibits stronger ATF4 transcription factor activity. ZNF281, a zinc finger protein involved in gene expression and cell proliferation regulation, may be related to the inflammatory response and cell proliferation in EC-5, further supporting the role of EC-5 in non-healing DFU [52]. This suggests that in the DFU healing process, endothelial cells play multiple roles in the development and healing of DFU, not only participating in local inflammatory responses but also potentially promoting healing by regulating metabolic activities. In particular, in a diabetic environment, metabolic abnormalities may affect endothelial cell function, leading to impaired wound healing [53].

Through CytoTRACE analysis, we found that endothelial cell subsets EC-3 and EC-4 exhibit significant stemness characteristics, while EC-5 and EC-8 have relatively lower stemness. This result indicates that there are significant differences in the differentiation potential of different endothelial cell subsets, which may directly affect their functional properties in the wound healing process. The differences in stemness not only reflect the differentiation status of the cells but may also reveal their potential roles in tissue repair [54]. For example, EC-3 and EC-4 with high stemness may play key roles in the early stages of healing, participating in rapid cell proliferation and migration; while EC-5 and EC-8 with low stemness may function in the later stages of healing, mainly involved in maintaining vascular stability and function. Further Monocle3 analysis revealed the dynamic changes of endothelial cells during development: EC-5 with non-healing characteristics is mainly enriched in the early stages of development, while EC-8 with healing characteristics is primarily distributed in the later stages. This finding provides important clues for understanding the temporal role of endothelial cells in the healing of diabetic foot ulcers (DFU). Specifically, the high abundance of EC-5 in the early stages may be related to its involvement in acute inflammatory responses, and its significant enrichment in DFU patients suggests that it may mediate inflammatory processes under pathological conditions. As the healing process progresses, EC-8 gradually takes the dominant position, indicating its key role in promoting angiogenesis and tissue regeneration. Additionally, MultiNicheNet analysis further revealed the different functions of EC-5 and EC-8 in cell signaling: EC-5 is mainly associated with inflammatory responses and cell death, while EC-8 is closely related to cell adhesion, migration, and cell cycle progression. These findings suggest that by targeting the regulation of these signaling pathways, it may be possible to effectively promote tissue healing or inhibit pathological angiogenesis and inflammatory responses [55]. Intervention strategies targeting different endothelial cell subsets and their signaling pathways may provide new ideas for improving wound healing in DFU patients.

SH3BGRL3, also known as SH3 domain-binding glutamic acid-rich-like protein 3, is a protein that plays a key role in cytoskeletal reorganization and signal transduction. Recent studies have shown that SH3BGRL3 plays a significant role in the development and wound-healing process of diabetic foot ulcers (DFU). Through high-dimensional weighted gene co-expression network analysis (HdWGCNA), we identified SH3BGRL3 as a key gene associated with the progression of DFU. Experimental results further confirm that under high-glucose conditions, SH3BGRL3 can significantly promote the angiogenesis, proliferation, and migration of human umbilical vein endothelial cells (HUVECs), which is crucial for the impaired wound healing process in a diabetic environment. Furthermore, in a diabetic mouse model, overexpression of SH3BGRL3 significantly accelerated wound healing, which may be related to its ability to promote angiogenesis and reduce inflammatory responses. These findings echo studies in tumor biology, where increased SH3BGRL3 expression has been found in various cancers, such as lung adenocarcinoma, bladder cancer, and renal cell carcinoma. In patients with urothelial carcinoma, its protein levels in urine are positively correlated with tumor high grading and invasiveness [56, 57]. Studies have shown that the SH3BGRL3 protein may interact with EGFR and inhibit its degradation, thereby enhancing EGF-induced AKT signaling [58]. Additionally, a recent study found that SH3BGRL3 interacts with the dynamin Myo1c. In SKBR3 cells, confocal microscopy observed the co-localization of SH3BGRL3 with Myo1c and ErbB2 on the plasma membrane, but co-immunoprecipitation experiments and mass spectrometry analysis indicated that SH3BGRL3 does not directly bind to ErbB2 but specifically recognizes the IQ domain-containing neck region of Myo1c, and this interaction is Ca²⁺-dependent. This suggests that the SH3BGRL3-Myo1c interaction may be a regulatory mechanism of cytoskeletal dynamics, playing an important role in cell migration [59]. Therefore, the role of SH3BGRL3 in DFU is not limited to promoting angiogenesis and cell migration but may also involve regulating inflammatory responses and cell proliferation. These functions make SH3BGRL3 a potential therapeutic target, providing new directions for the development of new treatment strategies for DFU. Future research needs to further explore the specific mechanisms of action of SH3BGRL3 and how to modulate its function through drug intervention, thereby providing more effective strategies for the treatment of DFU.

In summary, our study reveals the heterogeneity, differentiation status, and signaling communication of vascular endothelial cells in DFU, as well as genes related to DFU healing. These findings provide new insights into the pathogenesis of DFU and potential molecular targets for the development of new treatment strategies. Future research needs to further explore the specific functions of these endothelial cell subsets and their roles in the development and healing process of DFU.