Diabetic foot ulceration (DFU) is a major problem in diabetic patients as more than 15% of them are expected to develop DFUs within their lifetime. DFUs signifificantly impair quality of life, lead to prolonged hospitalization, and result in more than 70,000 lower extremity amputations per year in the USA alone1. Notably, more than half of the patients undergoing amputation due to DFU are expected to die within 5 years,  mortality rate which is higher than most cancers2. With the expected increase of Diabetes Mellitus (DM), DFUs will represent an even bigger burden for health systems worldwide and may prove to be one of the costliest diabetes complications3.

      Impaired wound healing leading to the development of chronic wounds in diabetic patients manifests exclusively in the foot in the presence of neuropathy and/or vascular disease4,5. Various cell types, including endothelial cells, fifibroblasts, keratinocytes, and immune cells play an important role in the wound healing process, but little is understood about their involvement in impaired wound healing in DFU. Dissecting cell differences within the foot ulcers between DFU patients whose ulcers heal and those who fail to heal and go on to develop a chronic ulcer, the differences between DM patients and non-DM healthy controls, and the differences between foot with DFU and intact forearm skin in both DM and healthy subjects, along with differences in blood immune cells, can considerably increase our understanding of DFU pathogenesis/healing.

      Single-cell RNA-sequencing (scRNASeq) analysis provides deep insight into cell function and disease pathophysiology by allowing the profifiling of the transcriptome landscape of individual cells in heterogeneous tissues. Currently, scRNASeq is widely used in the complex ecosystems of various cancers to map their microenvironment and discover molecular mechanisms and therapeutic targets6, and concerted efforts of the human cell atlas initiative aim to fully profifile all tissues of the human body7. Initial studies in our groups have reported that DM and DFU patients have increased inflflammatory cells and different fifibroblast clusters with a distinctive injury response-associated gene expression profifile, which is believed to be the result of DM related chronic inflflammation8. Spatial transcriptomics (ST) is a more recent method that enables the visualization and quantitation of the transcriptome in individual tissue sections, retaining spatial molecular information unlike scRNASeq9.

       In the present study, we primarily focused on differences between DFU patients whose ulcers heal (DFU-Healers) and those who don’t heal (DFU-Non-healers) within 12 weeks. We hypothesize that diabetic patients with impaired wound healing have aberrant gene and protein expression profifiles that lead to dysregulation of epithelial remodeling and inflflammation pathways. To this end, we investigated the molecular changes via scRNASeq analysis of DFUs and forearm skin biopsies, and peripheral blood mononuclear cells (PBMCs) from patients with healing and non-healing DFUs. As control group we also performed scRNASeq analysis of the foot and forearm biopsies and PBMCs from DM patients with no DFU, and healthy non-DM patients. We also studied different sites of a chronic wound (wound site, wound periphery, and healthy skin) to validate our findings We fifinally employed immunostaining and ST on DFU sections, as well as performed in vitro experiments to confifirm our most striking fifindings associated with DFU wound healing.

Results

DFU healing is signifificantly associated with a subset of fifibroblasts.

     To identify local and systemic factors associated with DFU healing, we examined the cellular landscape of DFUs by scRNASeq analysis of skin specimens from DFU, foot, forearm, and PBMC samples. We analyzed 54 samples from 17 diabetic patients (11 with and 6 without DFU) and 10 healthy non-DM subjects. The study cohort, objectives, and analysis strategy are outlined in Fig.1.a In total, we sequenced 174,962 cells (94,325 from foot, 37,182 from the forearm, and 43,455 from PBMC samples) and created a gene expression matrix for each cell, which we used to perform dimensionality reduction by UMAP and graph-based clustering, thereby identifying 37 orthogonal clusters of cells. The expression of established cell-specifific marker genes assisted in the annotation of these 37 cell clusters into 21 distinct cell types (Fig. 1b, c). We identifified most of canonical cell types observed in the human skin10,11 and PBMCs12, namely: smooth muscle cells, SMCs (TAGLN+, ACTA2+); fifibroblasts, Fibro (DCN+, CFD+); vascular endothelial cells, VasEndo (ACKR1+); T-lymphocytes, T-lympho (CD3D+); CD14+ monocytes, CD14-Mono (CD14+, S100A9+); differentiated keratinocytes, DiffKera (KRT1+, KRT10+); basal keratinocytes, BasalKera (KRT5+, KRT14+); natural killer cells, NK (CCL5+, GZMB+); NK and T cells, NKT (CD3D+, CCL5+); CD16+ monocytes,CD16-Mono (FCGR3A+/CD16+); M1 macrophages, M1-Macro (IL1B+); M2 macrophages, M2-Macro (CD163+); melanocytes and Schwann cells, Melano/Schwann (MLANA+, CDH19+); sweat and sebaceous gland cells, Sweat/Seba (DCD+); lymphatic endothelial cells, LymphEndo (CCL21+), erythrocytes, Erythro(HBB+); dendritic/Langerhans cells, DCs (GZMB+, IRF8+); Blymphocytes, B-lympho (CD79A+, MS4A1+); plasma cells, Plasma (MZB1+), and mast cells, Mast (TPSAB1+) (Fig. 1c).Comparative analysis of cell type abundance revealed substantial variations in the enrichment across clinical groups (Fig. 1d). Statistical analysis on cellular abundance showed signifificant variation for HE-Fibro, NKT, plasma and erythrocytes among the clinical groups (Supplementary Fig. 1). The enrichment of mast cells in non-diabetic subjects is in agreement with our previous work demonstrating excessive degranulated mast cells in diabetic skin, which could affect their ability to survive enzymatic digestion and sequencing13. Erythrocytes were increased in DFU-Nonhealers, most probably due to insuffificient RBC lysis during sample processing. Plasma cells were also enriched in DFU-Nonhealers, reflflecting a possible link of B-lymphocyte differentiation with non-healing wounds. Interestingly, our analysis showed signifificant heterogeneity in the transcriptome profifile of fifibroblasts and identifified a unique population of fifibroblasts that were overrepresented in the samples from DFU-Healers (Fig. 1b, d). We will refer to these as Healing Enriched-Fibroblasts, HE-Fibro, in the rest of the article. Further, the gene signature for each cell cluster was defifined by comparing the expression profifile of the target cluster with the rest of cells based on non-parametric Wilcoxon Rank Sum test (P-value <0.01 and Fold Change >2) (Fig. 1e). In-depth analysis of the HE-Fibro cell cluster revealed high expression of multiple extracellular matrix (ECM) remodeling (MMP1, MMP3) and immune/inflflammation (CHI3L1,TNFAIP6) associated genes (Fig. 1f). MMP1 (Matrix Metalloproteinase-1) interacts with CD49b14,15, an integrin alpha subunit involved in cell adhesion and cell-surface-mediated signaling in T, NK, and NKT cells16, fifibroblasts, and platelets. CHI3L1 (Chitinase-3-Like Protein 1) is a secreted glycoprotein that has been previously associated with pathogenic processes related to inflflammation and ECM remodeling17. TNFAIP6 (tumor necrosis factor alpha induced protein 6) is known to be involved in ECM stability and cell migration, and its expression is correlated with proteoglycan synthesis and aggregation18. This protein has shown anti-inflflammatory effects in various models of inflflammation, which suggest that it is a component of a negative feedback loop capable of downregulating the inflflammatory response19. The distinct and previously undescribed subtype or state of fifibroblasts, HE-Fibro, with overexpression of matrix remodeling, immune and inflflammatory genes, may contribute to successful wound repair in DFU-Healers. The top 10 overexpressed genes from the annotated cell clusters are included in Supplementary Dataset 1.

       Exploring cellular heterogeneity across different anatomical sites. To assess tissue specifific cellular heterogeneity along with gene expression and molecular pathway alterations, we generated the UMAP rendering split based on foot, forearm, or PBMC samples (Fig. 2a). The analysis depicted signifificant variations in the abundance of cell types based on anatomical sites (Fig. 2b). Fibroblasts, smooth muscle cells, melanocytes, sweat gland cells, vascular and lymphatic endothelial cells were enriched in the foot samples. The analysis on keratinocytes revealed a predominance of basal and differentiated keratinocytes in the foot and forearm samples, respectively (Fig. 2a, b). Interestingly, 99.94 ± 1.58% (mean ± stderr) of HE-Fibro cells, were identifified in the foot samples (Fig. 2a, b, Supplementary Table 1), indicating that these are foot specifific cells. The analysis demonstrated that a signifificant fraction of immune cells was contributed by the PBMC samples. CD14+ (98.47 ± 1.06%) and CD16+ (99.63 ± 1.11%) monocytes were observed predominantly in the PBMC samples. Dendritic, NK and NKT cell populations were also predominantly present in PBMC samples. The highest abundance of M2 macrophages (71.28 ± 1.27%) and mast cells (82.51 ± 1.49%) was observed in the foot samples, with lower proportions in the forearm samples, and the lowest in the PBMCs. On the adaptive immune system side, 31.27 ± 0.27% of T-lymphocytes came from the foot samples, 12.42 ± 0.35% from the forearm, and the remaining 56.31 ± 1.21% from the PBMCs (Fig. 2b, Supplementary Table 1). In contrast, more than half of plasma cells were derived from foot (66.01 ± 1.89%), and the remaining were almost equally proportioned between forearm and PBMCs. Most of B-lymphocytes (84.49 ± 1.63%) originated from PBMCs.

       To more closely examine the gene expression landscape of cells with differential abundances between foot and forearm i.e., fifibroblasts and keratinocytes, we performed comparative analysis on their transcriptome profifiles (Fig. 2c–f). The foot fifibroblasts exhibited upregulation of multiple genes associated with ECM remodeling and immune response. This may be attributed to the enrichment of HE-Fibro population in the foot samples (Fig. 2c). Genes that were overexpressed in foot fifibroblasts include the gene for Wnt signaling antagonist, secreted frizzled-related protein 4 (SFRP4) and genes directly related to ECM organization, asporin (ASPN) and tenascin C (TNC). Wnt signaling is crucial for effective wound healing20,21 and its modulation is closely linked with TGFβ expression22, which is in line with the enhanced expression of TGFB1 in these fifibroblasts. Tenascin C is known to upregulate TGFB1 as well as promote expression of type I collagen in fifibroblasts, which is essential for maintaining ECM integrity22. The cellular development (TSPAN8, WIF1) and immune cell traffificking (CCL19) related genes were signifificantly overexpressed in the fifibroblasts from the forearm (Fig. 2c). Pathway analysis on the foot fifibroblasts’ differentially expressed genes (DEGs) revealed signifificant (P value <0.01) activation of ILK, leukocyte extravasation signaling, RhoA signaling, and actin cytoskeleton signaling (Fig. 2d). The comparative analysis between foot and forearm keratinocytes showed signifificant upregulation of basal (KRT6A, KRT16, KRT17) and differentiated (KRT2, KRT10) keratinocyte associated genes in the foot and forearm samples, respectively (Fig. 2e). This discrepancy can be explained by the fact that forearm biopsies represent unwounded tissues with fully stratifified epidermis, as opposed to foot samples that include DFUs with partially formed epithelium, and therefore fewer differentiated keratinocytes. Moreover, the differences between plantar glabrous skin and forearm hairy skin could contribute to the disparity23. In addition to upregulation of alarmins like KRT6A/16/17 in the foot samples24, we also observed upregulation of inflflammatory molecules including S100A8 and S100A9, known to activate the immune system in response to skin injury25. Further pathway analysis on foot keratinocyte DEGs uncovered signifificant activation of immune and inflflammatory pathways including ILK and IL-8 signaling (Fig. 2f).

      Systemic dysregulations revealed by comparative analyses of PBMCs across clinical groups. To better understand the impact of DFU at a systemic level, we performed separate analysis on PBMC samples alone from the four clinical groups, viz Healthy (healthy subjects without DM), DFU-Healer (DM patients with healing DFUs), DFU-Non-healer (DM patients with non-healing DFUs), and Diabetic (DM patients without DFU) (Fig. 3a). The cell annotation was done using well-established marker genes (Supplementary Fig. 2). The DFU-Healers were observed to have higher proportions of naive and early differentiated progenitor Tlymphocytes, T-lympho, expressing CCR7, shown to have a role in activation of various T-cell subsets26 (Fig. 3b). On the other hand, DFU-Non-healers had a higher proportion of cytotoxic NKT cells (IL7R, GZMB, KLRD1), indicating a shift in T-cell subpopulations correlating with DFU healing (Supplementary Fig. 3). We observed statistically signifificant higher CCR7+ T-lympho cells to NKT cells ratio (P value < 0.01) in the DFUHealers as compared to DFU-Non-healers and DM patients without DFU, indicating the association of these T-cells with successful wound healing (Fig. 3b). A signifificantly higher proportion of CCR7+ CD8+ T cells (cluster CD8T2 in Fig. 3a) was also observed in DFU-Healers as compared to DFU-Non-healers (Fig. 3c). Further DEGs analysis on these T-lympho, CD8T2 and NKT cells indicated overexpression of T-cell-specifific genes like IL7R, TCF7, and CCR7 in the DFU-Healers, whereas DFU-Nonhealers overexpressed NKT lineage genes like NKG7, GNLY, CCL5, and KLRD1 (Fig. 3c). Pathway analysis on these T/ NKT cells DEGs demonstrated inhibition of key immune and inflflammation pathways including IL-6, IL-8, CD28 Signaling in T helper cells, and iCOS-iCOSL pathways and activation of RhoGDI and EIF2 signaling in the DFU-Healers, as compared to DFU-Non-healers at the systemic level (Fig. 3d). Further systems biology analysis revealed inhibition of several upstream regulators of immune pathways including CD44, TGFB1, CCL5, and NFKBIA in the T-cells from PBMCs of patients with healing DFUs (Fig. 3e). This was in accordance with the observed reduced gene expression of NFKBIA, CCL5, and TGFB1 in DFU-Healers compared to high expression in DFU-Non-healers (Fig. 3f). In aggregate, these results underscore the enrichment of naive T-cells with a prevalence of immune inhibitory pathways and processes for DFU-Healers, and a state of chronic inflflammation for DFU-Non-healers, at the systemic level.

       T, NK, and NKT cells exhibit distinct cell subpopulations in DFU-Healers and DFU-Non-healers. The focused sub-clustering analysis on the T, NK, and NKT cell populations identifified 17 subclusters (Supplementary Fig. 4a). CD4+ (sub-clusters 0, 4, 10) and CD8+ (sub-cluster 14) naive T-cells (CCR7+, LEF1+), that can self-renew and proliferate readily into other T-cells, were enriched in DFU-Healers (Supplementary Figs. 4b, c and 5c, d). Cluster 6, CD8+ effector T cells (CCL5+, GZMB+, IL32+, GZMK+), enriched in DFU-Healers, also expressed higher levels of CD27, a key molecule in generation and maintenance of T-cell immunity27. NKT (CD8+, CCL5+, GZMB+, IL32+, GZMH+) cells subclusters 5 and 7 were enriched in Diabetic and DFU-Non-healer groups, respectively (Supplementary Figs. 4a, b and 5c, d). The DFU-Non-healer enriched sub-cluster 7 also had high expression of T-cell exhaustion marker, TIGIT (Supplementary Fig. 5d). Sample site-specifific split t-SNE plots revealed separate clustering of T, NKT and NK cells from skin and PBMCs (Supplementary Fig. 5a). Sub-clusters 1, 2, 9, 11, and 13, expressing activation markers CD69 and CD44 were largely made up of skin samples derived from foot (Supplementary Fig. 5a–d). Clusters 1 and 9 were positive for T-cell exhaustion markers (TIGIT+, HAVCR2+, LAG3+) (Supplementary Fig. 5c). In PBMCs, DFU-Healers appeared to have more non-polarized central memory and naive T-cells (Supplementary Fig. 4b, c; Supplementary Fig. 5). CD27 which characterizes central memory T-cells that lack immediate cytotoxicity28, was also more in the DFU-Healers (Supplementary Fig. 5c). In contrast, DFU-Non-healers were enriched with cytotoxic NKT cells (cluster 7), expressing GZMH, GZMA, and GZMB (Supplementary Figs. 4b and 5c). These granzyme molecules have been previously implicated in impaired wound healing development by promoting chronic inflflammation, vascular dysfunction, and reduced cell adhesion29. We also observed a unique CD4+ cluster (cluster 10, Supplementary Fig. 4a) that was predominantly present in DFU-Non-healers and enriched for GIMAP1 and GIMAP4, both shown to be implicated in T helper cell differentiation towards the Th1 lineage30. In the Diabetic group, a GZMH+, GNLY+, and CCL5+ expressing cluster (cluster 5) was prominent, pointing toward the presence of specialized DM associated NKT cells. CCL5/RANTES, a potent chemoattractant of immune cells, has been reported to be strongly downregulated in DFUs compared to acute wounds, and could represent a potential therapeutic target31. In summary, while skin samples derived T/NKT cells did not show signifificant differences between DFU-Healers and DFU-Non-healers, potentially due to the low number of recovered T-cells, a defifinitive enhancement of naive T-cells was seen in PBMCs of DFU-Healers compared to more cytotoxic NKT cells in DFU-Non-healers.

      Analysis of foot ulcer cells reveals the signifificance of localized inflflammation in diabetic wound healing. To map the transcriptome and cellular landscape of the site for DFUs, we performed focused analysis on single-cell profifile of 94,325 cells from 26 foot samples. Split UMAP analysis indicated differential abundance of cell types among the four clinical groups (Fig. 4a, Supplementary Fig. 6). The DFU-Healers had a signifificantly higher number of HE-Fibro cells (P value <0.05) as compared to DFU-Non-healers, Diabetic patients and non-DM healthy controls (Fig. 4b). Additionally, the DFU-Healer group also showed a signifificantly higher proportion of M1 macrophages (classically activated macrophages that promote inflflammation) than M2 macrophages (alternatively activated macrophages with anti-inflflammatory properties), as opposed to DFU-Non-healers (Fig. 4c). Also, a group of SMCs, SMC2, with overexpression of proliferation markers CENPF, PTTG1, MKI67, and TOP2A was signifificantly enriched in DFU-Healers (Fig. 4d, Supplementary Dataset 2), highlighting the presence of a highly proliferative response for promoting wound healing. On the other hand, DFU-Non-healers macrophages overexpressed genes from the complement system like C1QA/B/C, which are associated with M2 macrophage-like anti-inflflammatory responses32 (Fig. 4f). Immunoflfluorescent staining of healing and non-healing DFUs with pan-macrophage marker CD68, M1 markers IL1B and S100A8 and M2 markers DAB2 and CD163 confifirmed more M1 associated macrophages in DFU-Healers with both M1 markers showing increased presence and more M2 macrophages in DFU-Non-healers (Supplementary Fig. 7a–f). SMC cluster in healing DFUs. Other cell types also exhibited variation across clinical groups but did not achieve statistical signifificance due to intragroup variation among patients(Fig. 4e). DEGs analysis on DFU-Healers vs. Non-healers and M1 macrophages vs. M2 macrophages identifified a signature comprising of 195 genes that were differentially expressed in M1 macrophages from DFU-Healers (Fig. 4f).  Immunoflfluorescent staining of healing and non-healing DFUs with pan-macrophage marker CD68, M1 markers IL1B and S100A8 and M2 markers DAB2 and CD163 confifirmed more M1 associated macrophages in DFU-Healers with both M1 markers showing increased presence and more M2 macrophages in DFU-Non-healers (Supplementary Fig. 7a–f).

     Pathway analysis showed activation of the IL17 signaling pathway, a known regulator of inflflammatory response33, in DFUHealers (Fig. 4g). The upstream regulators activated in DFUHealers included HIF1A, TNF, STAT5a/b, TLR7, TLR9, and IL17R/C (Fig. 4h), whereas SOX4, TGFB1, and NANOG, were inhibited (Fig. 4i). Immunoflfluorescent staining with IL17 and HIF1A antibodies showed a trend for more IL17+ cells (p = 0.06) and higher HIF1A expression (p = 0.09) in DFU-Healers (Supplementary Fig. 7g–i).

      Similar analyses were also conducted on the forearm cells (Supplementary Figs. 8a and 9). Differentiated keratinocytes were enriched in DFU-Healers compared to DFU-Non-healers (Supplementary Fig. 8b). We found that LGALS7 or Galectin-7, which has been previously implicated in keratinocyte migration during re-epithelialization of wounded epidermis34, was the top differentially expressed gene in the forearm keratinocytes of DFU-Healers (Supplementary Fig. 8c).

       Healing associated fifibroblasts drive DFU healing by promoting matrix remodeling and inflflammatory response. To further delineate the role of fifibroblasts in wound healing, we performed focused analysis on fifibroblasts that produced 14 sub-clusters, representing different molecular states or subtypes of fifibroblasts (Fig. 5a, Supplementary Dataset 3 for top 10 marker genes). The majority of sub-clusters showed distinct expression profifiles indicating heterogeneity in the fifibroblast population (Fig. 5b). Sub-clusters 0, 1, 2, and 5 comprised most of the cells from unwounded skin. Sub-cluster 0 was characterized by the expression of reticular fifibroblast marker MGP35 and multiple adipocyteassociated genes (APOE, APOD, CFD), consistent with the enhanced adipogenic potential of these cells36. Sub-clusters 2 and 5 contained cells expressing papillary fifibroblast markers PTGDS, APCDD1, and COL23A137,38, while sub-cluster 1 was enriched for WISP2, PI16, SLPI, and SFRP2, which describe fifibroblasts residing both in the papillary and reticular dermis, and are believed to contribute to ECM homeostasis39–41. The evaluation of cellular makeup of clusters unveiled a higher proportion of cells (58–90%) from DFU-Healers in specifific sub-clusters; clusters 3, 4, 6, and 13. These four sub-clusters represent four heterogeneous states or subtypes of the HE-Fibro (Fig. 5a, marked area). Further generation of gene signatures for these sub-clusters revealed that cluster 3 was signifificantly enriched with cells expressing genes related to ECM remodeling including MMP1 and MMP3 (Fig. 5b). Matrix metalloproteinases MMP1 and MMP3 have been well-known early responders to tissue injury, actively regulating the inflflammatory phase of healing by degradation of the ECM, stimulating leukocyte infifiltration for resolution of inflflammation and transition to the proliferative phase42. Cluster 4 exhibited overexpression of POSTN and ASPN (Fig. 5b) that are associated with ECM signaling, adhesion and migration. POSTN (Periostin) is a ligand for alpha-V/beta-3 and alpha-V/beta-5 integrins and supports adhesion and migration of epithelial cells43, and has been shown to play a regulatory role in fifibroblast proliferation and inflflammation44,45. ASPN (Asporin) is an ECM protein that has been found to inhibit TLR2- and TLR4- induced NF-κB activity and pro-inflflammatory cytokine expression in macrophages46. TLR4 mediated inflflammation drives the synergistic effect of hypoxia and hyperglycemia on impairment of diabetic wound healing47, hence overexpression of ASPN might be an important determining factor for healing of DFUs. In a recent study, a distinct ASPN and POSTN enriched cluster of fifibroblasts was described as mesenchymal and shown to have a more reticular dermis localization40. These sub-clusters (3, 4, 6,and 13) were also enriched with genes like IL6, CHI3L1, PLA2G2A, and TIMP1, commonly associated with an inflflammatory signature (Fig. 5b). Based on our analysis, we identifified a healing associated fifibroblast signature consisting of ECM remodeling and inflflammatory response-related genes: MMP1, MMP3, IL6, CHI3L1, ASPN, POSTN, and PLA2G2A.

      Further analysis of these fifibroblasts revealed that IL6/TIMP1/PLA2G2A and CHI3L1 transcripts were detected simultaneously in ~38% of the cells suggestive of a common regulatory mechanism in HE-Fibro. We also noticed that 99.8% of CHIL3L1 expressing cells exhibited signifificant expression of at least one of the ECM remodeling genes including MMP1, MMP3, MMP11, indicating a role of these genes in tissue repair. Based on these preliminary results we posit that CHIL3L1 is one of the key players in driving the healing phenotype of HE-Fibro by expressing pro-inflflammatory and ECM genes together to improve wound repair. Several lines of evidence have previously implicated CHI3L1 in dampening of chronic inflflammation48, promoting M1 macrophage activation49 and stimulating fifibroblast proliferation50 and ECM remodeling51. DEGs analysis on  these DFU-Healers vs. other fifibroblast clusters revealed some ubiquitous markers (HIF1A, TNFAIP6) that are overexpressed in HE-Fibro cells (Fig. 5).

      Pathway analysis indicated the activation of multiple immune and inflflammatory pathways including IL6, HIF1A, and ILK signaling in the fifibroblasts from DFU-Healers (Fig. 5c). Moreover, upstream regulator genes like TNF, HIF1A, and IL6, were activated (Fig. 5d) in the DFU-Healers. HIF1A (Hypoxia inducible factor 1-alpha) is a master-regulator that activates multiple factors to enhance wound healing by promoting cellular motility and proliferation, angiogenesis, re-epithelialization, and cell survival52. HIF1A also upregulates IL6 expression by binding to its promoter region53, thereby promoting inflflammation and cell proliferation.

      IL6 is a pleiotropic cytokine that plays a vital role in wound healing, as studies have shown delayed and impaired wound closure in IL6 knockout mice54–56. It is noteworthy that in our previous work with a large animal model of diabetic wound healing57, IL6 levels post-injury were attenuated compared to acute non-diabetic wounds, suggesting that increased IL6 is advantageous for DFUs. Another notable gene in the molecular interaction network is TNF (Fig. 5d), a potent pro-inflflammatory cytokine that has been previously implicated in the wound healing process58, and is known to be elevated shortly after wounding. TNF has also been shown to upregulate the expression of MMP1 and MMP3 in human dermal fifibroblasts via NFκB/p65 activation59. To corroborate the activation of central regulators TNFA and IL6 at the protein level, we probed healing and nonhealing DFUs and measured the percentage of area stained to fifind signifificantly higher IL6 expression (p = 0.001) and a trend for higher TNFA expression (p = 0.08) (Supplementary Fig. 7j–l) in DFU-Healers.

      Deciphering communication among Healing associated fifibroblasts. Further, to determine possible communication among heterogeneous healing associated fifibroblasts (HE-Fibro subclusters 3, 4, 6, and 13; Fig. 5a), we performed ligand, receptor, and target gene co-expression analysis using the NicheNetR algorithm60. NicheNet predicts which ligands from one or more cell population(s), termed “sender/niche”, will most likely affect gene expression in interacting cell population(s), termed as “receiver/target”. Also, this algorithm can predict which specifific target genes in the “receiver” cell populations are affected by the predicted ligands in the “sender” cell population(s). As subcluster 3 of the HE-Fibro was enriched in both inflflammatory and ECM remodeling genes (IL6, CHI3L1, MMP1), critical to the healing process, it was selected as “sender” cell population to generate ligand candidates, while the remaining HE-Fibro subclusters (4, 6, and 13) were treated as receiver cells (Supplementary Fig. 10). To fifilter out non-specifific ligand and receptors, we also included control fifibroblast sub-clusters (0, 2, and 5), enriched in healthy non-DM and diabetic without DFU patients, as receiver cells. The analysis identifified multiple ligands including IL6, CCL2, and TIMP1 with high correlation between differential expression of ligands in “sender” fifibroblasts and their target genes in the “healer” fifibroblasts (sub-clusters 4, 6, and 13) but not in the “control” fifibroblasts (sub-clusters 0, 2, and 5) (Fig. 5e) This indicates that fifibroblasts from sub-cluster 3 are primarily interacting with the other three sub-clusters of HE-Fibro subset enriched in DFU-Healers rather than the “control” fifibroblasts subset enriched in healthy non-DM and diabetic with no-DFU FN1 was enriched in all the HE-Fibro sub-clusters (3, 4, 6, and 13), while IL6, MMP13, CCL2, PTGS2, and VEGFA were enriched in only HE-Fibro sub-cluster 3 (Supplementary Fig. 11).

       A heatmap displaying the connection between these key ligands that are expressed by HE-Fibro sub-cluster 3 (rows) and marker genes for the HE-Fibro sub-clusters (3, 4, 6, and 13) (columns) is shown in Fig. 5f and Supplementary Fig. 12. For example, CCL2 overexpression in the “receiver” HE-Fibro sub-clusters can be strongly predicted by the presence of IL6 in the “sender” HEFibro sub-cluster, indicating possible molecular interaction between them. The circos plot shows association between ligands from “sender” cells, sub-cluster 3 (lower hemicircle), and DEGs in “receiver” sub-clusters 4, 6, and 13 (upper hemicircle) (Fig. 5g, Supplementary Fig. 13). The analysis identifified IL6, MMP13, CCL2, CXCL12, CTGF, TIMP1, and VEGFA as key regulatory ligands in the HE-Fibro sub-cluster 3, altering the expression of downstream target genes in the HE-Fibro sub-cluster 4, 6, and 13.These identifified ligands and their downstream targets might be responsible for the healing associated phenotype of HE-Fibro. Based on enrichment and specifific regulatory interaction among HE-Fibro sub-clusters in DFU-Healers, we postulate that their role consists of creating a benefificial physiological environment for accelerated DFU healing.

      RNA velocity analysis predicts differentiation of HE-Fibro to other fifibroblasts and SMCs. To gauge the transcriptional dynamics of cell types of interest, we performed in silico trajectory analysis by computing the RNA velocity of Fibro, HE-Fibro, SMC, and M1- and M2-Macro cell populations in DFUs of Healers and Non-healers. We observed that the DFU-Healers cells were more flfluid and interconnected, with distinct thick bridges between SMC and HE-Fibro or Fibro, reflflecting potential transdifferentiation events. The analysis also predicted dual transdifferentiation of HE-Fibro toward SMCs and fifibroblasts, as well as transdifferentiation of SMCs to macrophages.

         Conversely, the DFU-Non-healers displayed more clearly defifined branches of cellular trajectory and weak differentiation bridges among the cell types (Fig. 5h). Further, predicting the individual cells latent times revealed universally lower latent times for Healers cell types as compared to Non-healers cell types (Supplementary Fig. 14), thereby substantiating our hypothesis that Non-healers are stuck and fail to progress in the wound healing cascade. The lower latent times predicted for cell types from Healers indicate the likelihood that these cells are in earlier stages of differentiation and more stem cell-like.

        HE-Fibro display enrichment at the early recovery time point in acute wounds. Our single-cell data identifified signifificant enrichment of HE-Fibro and M1-Macro in the DFUs that eventually healed (Fig. 4b, c). To explore this enrichment in sequential wound samples during the healing process, we performed enrichment analysis in a previously published microarray dataset of acute wound healing61 by calculating a score based on the HEFibro and M1-Macro gene signatures. The analysis revealed signifificant enrichment of HE-Fibro (Adj. P-Value = <0.0001) and M1-Macro (Adj. P-Value = <0.0001) at Days 4–8 post-injury samples in comparison to healthy unwounded skin (Supplementary Fig. 15), which correspond to the inflflammation and early proliferation stages of healing. We did not observe enrichment of HE-Fibro and M1-Macro during the later phases of wound healing. These results provide further evidence that DFUHealers exhibit a more acute-like response and have advanced forward in the wound healing process. Thus, enrichment of HE Fibro and M1-Macro during the initial phases appears crucial for successful wound healing.

       HE-Fibro display enrichment at the early recovery time point in acute wounds. Our single-cell data identifified signifificant enrichment of HE-Fibro and M1-Macro in the DFUs that eventually healed (Fig. 4b, c). To explore this enrichment in sequential wound samples during the healing process, we performed enrichment analysis in a previously published microarray dataset of acute wound healing61 by calculating a score based on the HEFibro and M1-Macro gene signatures. The analysis revealed signifificant enrichment of HE-Fibro (Adj. P-Value = <0.0001) and M1-Macro (Adj. P-Value = <0.0001) at Days 4–8 post-injury samples in comparison to healthy unwounded skin (Supplementary Fig. 15), which correspond to the inflflammation and early proliferation stages of healing. We did not observe enrichment of HE-Fibro and M1-Macro during the later phases of wound healing. These results provide further evidence that DFUHealers exhibit a more acute-like response and have advanced forward in the wound healing process. Thus, enrichment of HEFibro and M1-Macro during the initial phases appears crucial for successful wound healing.

      Spatial transcriptomics and immunohistochemistry further elucidate gene expression patterns in healing and non-healing DFUs. We subsequently selected well-defifined surgically excised DFU sections from healers and non-healers for additional characterization. They both displayed blood vessel proliferation and chronic inflflammatory cell infifiltrates predominantly with perivascular distribution (Fig. 6a, c). We stained for inflflammatory fifibroblast markers, CHI3L1 and TIMP1, together with pan- fifibroblast marker fifibroblast activation protein (FAP) and discovered elevated numbers of triple-positive cells within the ulcer area of healing DFUs, with the cells forming dense aggregates (Fig. 6d). However, in the non-healing ulcers these cells were far fewer and scarcely distributed (Fig. 6b). We also evaluated gene expression using a spatial transcriptomics approach. The GeoMx® platform enables spatial, high-plex quantitation of gene expression in tissue through the use of in situ hybridization (ISH) probes that target mRNA in tissue; attached to the probes are photocleavable and indexed oligonucleotides that can be liberated via UV light and counted with an Illumina® sequencer. Regions of interest (ROIs) were chosen after staining for immune cell marker CD45, vasculature marker αSMA, and epithelial marker pan-Cytokeratin along with nuclear counterstain DAPI, to represent areas within the ulcer, at the edge of ulcer and adjacent non-injured tissue (Fig. 6a, c and Supplementary Fig. 16a, b). Hierarchical clustering analysis of representative healing and non-healing specimens revealed dissimilar gene expression pro- fifiles according to location within the sample: ROIs at similar dermal depth grouped together. (Fig. 6e, f). The non-healing ulcer ROI was particularly distinct from neighboring ROIs (Fig. 6e), while the healing ulcer ROIs appeared more transcriptionally similar (Fig. 6f). Focusing on the ulcer localized ROIs, DE analysis from two Healers (9 ROIs in total) and two Non-healers (4 ROIs in total) showed 148 genes upregulated in Healers and 57 in Nonhealers (Fig. 6g and Supplementary Fig. 16c–f for additional DFUs). Among the most notable ones, HE-Fibro marker PLA2G2A and M1 macrophage marker FOS were overexpressed in Healers (Fig. 6h, i), while M2 macrophage markers TYMP and ANXA1 were upregulated in Non-healers (Fig. 6j, k). Finally, gene ontology (GO) enrichment analysis unveiled cellular response to TNF as top biological function activated in healing ulcers and myeloid leukocyte migration in non-healing ulcers (Fig. 6l). Taken together, these fifindings verify our previous observations at the protein level and specify the location and functional roles of cell types reported in our scRNASeq dataset.

     To further validate the fifinding based on spatial profifiling that HE-Fibro mainly form niches in the wound bed to promote wound healing, we performed scRNASeq analysis on multiple samples from the same patient. ScRNASeq analysis was performed on skin specimens of the same patient from three different sites: wound bed, wound edge, and non-wound excess skin from a pressure sore excision (Supplementary Figs. 17 and 18). The unsupervised analysis and cellular annotation revealed that HE-Fibro were enriched in the wound bed, but not in the wound edge and non-wounded samples (Supplementary Fig. 18a). This unique wound bed enriched cell cluster exhibited signififi- cantly higher expression of HE-Fibro-associated genes like IL6, TNFAIP6, MMP1, and CHI3L1 (Supplementary Fig. 18c). The absence of any other fifibroblast cluster for the wound bed sample suggests that HE-Fibro originate from “normal” fifibroblasts. These results further affifirm an association of HE-Fibro with the wound healing process in an additional type of chronic wound and point toward heterogeneity of fifibroblasts across different regions of ulcers.

         Induced overexpression of CHI3L1 in dermal fifibroblasts inflfluences cell behavior. To explore the effects of inflflammatory marker genes expression in vitro, we selected one of the top enriched genes CHI3L1 and generated dermal fifibroblast cell lines transduced with lentiviral vectors overexpressing CHI3L1 (CHI3L1-OE) or a control sequence (CTRL). Western blotting demonstrated a complete lack of expression in untreated cells and RT-qPCR analyses confifirmed a signifificant upregulation of CHI3L1 with construct 2 (Supplementary Fig. 19a, b), which we selected for further experiments. In adhesion assays, more CHI3L1-OE cells attached to fifibronectin-coated surfaces compared to CTRL (Supplementary Fig. 19c, d), while diminished migration was observed in scratch wound experiments (Supplementary Fig. 19e, f). Altogether, these fifindings shed light on the potential functional roles of the HE-Fibro, indicating that they possess enhanced adherent and decreased migratory capacities and suggest that they are fifirmly anchored on the ECM and mediate healing through secretion of molecules.

Discussion

      In this study, we performed large-scale unbiased scRNASeq to accurately and systematically profifile patients with healing and non-healing DFUs, together with healthy non-DM subjects, and DM patients without DFUs, as controls. For a subset of patients, we also characterized forearm biopsies and PBMCs to evaluate any potential systemic effects of DM in presence of DFUs. To the best of our knowledge, we were the fifirst groups to employ this approach in DFU samples8, and we have now substantially expanded the number of cells sequenced, and incorporated stateof-the-art techniques like spatial transcriptomics, in order to gain novel insights into the transcriptomic landscape of DFU healing.

      We identifified a hitherto unreported fifibroblast cell type associated with healing and expressing multiple immune and ECM remodeling-related genes. We then corroborated the results at the protein level and with the additional sequencing modality of spatial transcriptomics, demonstrating their localization within the ulcer area. It has become increasingly apparent that dermal fifibroblasts are a diverse and highly heterogeneous population with different functional roles in wound healing62–66. Fibroblasts at sites of inflflammation, for instance within tertiary lymphoid structures, have been shown to acquire immune cell features67, while in murine wounds’ granulation tissue, a large proportion of fifibroblasts is of myeloid cell origin68. A number of studies have also emphasized the interactions between fifibroblasts and monocytes or macrophages in the context of inflflammation and wound healing, implying a reciprocal relationship69–71. Our data suggest that specifific fifibroblast subtypes are key players in healing of DFUs and targeting them could be a therapeutic option.

     Mapping the immune landscape of Healers and Non-healers revealed the presence of more M1 macrophages in Healers and M2 in Non-healers, as well as higher numbers of naive and central memory T-cells in Healers, as opposed to more NK and NKT cells in Non-healers. Accumulating evidence suggests that a favorable outcome in wound healing is contingent on a highly regulated balance of macrophage polarized states72. The presence of more M2 macrophages does not necessarily equate better healing, as wound repair studies have shown delayed healing in diabetic73 or wild-type74 mice treated with M2 macrophages. DFUs are most probably populated both by tissue-resident macrophages that have differentiated from bone marrow-derived monocytes with a minimal contribution of yolk sac originating macrophages75 as well as peripherally recruited monocytes responding to inflflammatory cues and similarly differentiating to macrophages at the injury site76. To disentangle the admixture of macrophage subsets in DFUs and fully characterize their origins and repopulation dynamics, defifinitive future studies are required in human skin as has been accomplished in other organs77.

       Impairment in the recruitment of macrophages and neutrophils in DFUs was recently demonstrated78. A dysregulation in the differentiation of peripheral blood-derived T cells and diminished T-cell receptor repertoire diversity has been previously reported for DFU patients79. The majority of T-lymphocytes in our study originated from the blood samples, while macrophages were mostly located at the foot. Interestingly, in PBMCs of healers, inflflammation pathways were mostly inhibited. These fifindings underline fundamental differences between systemic inflflammation and the local wound inflflammatory milieu. Overall, our results provide further evidence to support the claim that localized activated inflflammatory response is required to surmount the chronic inflflammation in DFUs, and progress to the next phases of wound healing80,81, while, conversely, inhibition of inflflammatory processes at the systemic level appears benefificial for healing. Trajectory analyses and comparison with acute wounds’ gene expression demonstrated that DFUHealers had advanced forward in the wound healing process and therefore the observed differences are not a distinct new pattern of healing but rather the difference between a dysregulated chronic inflflammatory environment in DFU-Non-healers and wounds that have progressed to the fifirst (inflflammatory) or early second (proliferative) phases of healing in DFU-Healers.

       Future longitudinal studies interrogating DFU samples collected from the same patient over multiple time points in the course of wound repair82 could help build a map of the diabetic wound healing timeline. However, considering how technically challenging it is to consistently achieve high enough quality in debridement samples for next-generation sequencing83,84, singlenucleus sequencing could be an alternative and complementary approach to offer single-cell resolution without relying on highly viable single-cell suspensions85,86.

       In summary, we present a comprehensive characterization of the DFU ecosystem and report novel cell types and interactions. Our dataset will be a valuable resource for diabetes, dermatology, and wound healing research, and can serve as the baseline for designing in vitro and in vivo experiments for the assessment of therapeutic interventions focusing on one or more cell types. Future studies utilizing pre-enrichment via flflow or magnetic cell sorting could further characterize specifific populations and lead to the discovery of rare cells.

Methods

Subjects. Our study includes non-DM patients (n = 10) who underwent foot surgery for various reasons, such as hallux valgus correction, as the healthy controls, and DM patients without foot ulceration (n = 6) who had similar foot surgery. Discarded skin specimens from the dorsum of the foot were collected for analysis. We also enrolled DM patients with plantar foot ulceration (DFU) (n = 11), who underwent surgical resection of the ulcer, providing suffificient wound and peri-wound tissue for analysis. Subjects with any conditions, other than DM, or medications that could affect wound healing were excluded from the study. Four non-DM subjects, two DM patients with no DFU, and fifive DM patients with DFU (Healers; n = 3, Non-healers; n = 2) provided two 3-mm forearm skin biopsies and 20 ml of blood, from which PBMCs were isolated, within 1 week of the foot surgery. DFU patients were followed for 12 weeks post-surgery and were divided into two subgroups: those who healed their ulcers and those who failed to heal them (Healers; n = 7, Non-healers; n = 4). Supplementary Dataset 4 includes clinical details of the subjects included in the study. There were no major differences among the main biological characteristics of the studied groups (Supplementary Table 2). All patients were enrolled and followed at the Joslin-Beth Israel Deaconess Foot Center, Boston, MA, and the study was approved by the Beth Israel Deaconess Medical Center IRB (Reference number 2018P000581). For the scRNASeq analysis of spatially separated samples, multiple samples were collected from an ischial pressure sore of one patient at the Yale Plastic and Reconstructive Surgery – Wound Center, New Haven, CT (collected under IRB approval 1609018360). Informed consent was obtained from all study participants at Beth Israel Deaconess Medical Center and Yale. For comparative analysis of enrichment of various healing associated cell types in sequential wound healing samples, we downloaded microarray burn study data of 5 patients from the ArrayExpress database, accession ID E-MTAB-1323.

       PBMCs isolation. PBMCs were separated using Ficoll-Paque density gradient fractionation, as previously described87, and cryopreserved in freshly prepared freezing media (90% FBS and 10% DMSO).

       Single-cell preparation from skin samples. Skin specimens were kept in sterile PBS on ice until processing, normally within 3 h post-surgery. The skin was cleaned by sequentially immersing in 10% Betadine, 70% ethanol, and PBS for 1 min at a time. Then it was incubated in 5 mg/ml Dispase II (Thermo Fisher Scientifific, 17105041) in HBSS (STEMCELL Technologies, 37150) overnight at 4 °C. The next day, the epidermis was peeled off using forceps, and the tissue was fifinely minced with a No. 10 disposable scalpel. The skin pieces were then placed in an enzyme cocktail consisting of 3.3 mg/ml Collagenase-P (Roche, 11249002001), 3.3 mg/ml Dispase II, and 1.5 mg/ml DNase I (STEMCELL Technologies, 07470) in 0.25% Trypsin-EDTA (Thermo Fisher Scientifific, 25200072) and incubated for 90 min at 37 °C with constant shaking, using glass pipettes for trituration every 20 min. Enzymes were then inactivated with the addition of complete DMEM (+10%FBS, +1% Pen/Strep). The single-cell suspension was passed through 70 and 40 μm cell strainers and centrifuged for 10 min, 500 × g at 4 °C. For red blood cell (RBC) lysis, ACK buffer (Lonza, 10-548E) was added. The process resulted in highly viable, typically >90%, single-cell suspensions. For immediate single-cell capture, the cells were resuspended in 0.04% Ultra-Pure BSA in PBS (Thermo Fisher Scientifific) and concentration was adjusted to 1000 cells/µl. If not processing for scRNASeq immediately, the cells were cryopreserved in freshly prepared freezing media (90% FBS and 10% DMSO).

       Single-cell RNA sequencing. The single-cell preparations of the foot, forearm, and PBMC samples were used fresh or after thawing of viably frozen samples with fifinal resuspension in PBS with 1% BSA. A droplet-based ultra-high throughput scRNASeq system was utilized to capture single cells along with uniquely barcoded primer beads together in tiny droplets, enabling large-scale parallel single-cell gene expression studies. The gene expression (GEX) libraries were prepared using the Chromium 3’V2/3 reagent kits (10x Genomics, 120237 and 1000075). Brieflfly, gel bead-in-emulsions (GEMs) were generated and barcoded by loading single-cell suspensions along with gel beads and reverse transcription (RT) master mix in 10x Genomics Single cell chip (A chip kit, 120236; B chip kit, 1000153) and running on the chromium controller (10x Genomics, 110211). Following RT, the cDNA was amplifified and used to generate GEX libraries. The cDNA and GEX libraries were quantifified using Qubit 3.0 flfluorometer (Life Technologies, 15387293), and quality was assessed using HS DNA chips (Agilent technologies, 5067-4627) with 2100 Bioanalyzer (Agilent Technologies, G2939BA). Sequencing was performed using massively parallel sequencing on the Novaseq S4 platform (Illumina). We produced~40,000–50,000 reads per cell capturing the expression of ~1000–2000 transcripts per cell.

       Data processing and analysis. Raw scRNASeq data was demultiplexed, aligned to the reference human genome (Hg38), and processed for single-cell gene counting using the Cell Ranger Software from 10X Genomics Inc. The single-cell count data was normalized using the SCTransform algorithm in Seurat v3.0 Bioconductor package88 that uses regularized negative binomial models for normalizing sparse single-cell data. The normalized expression profifiles of the samples were merged, and undergone quality control, pre-processing, unsupervised and supervised analysis using various R and Bioconductor packages. The quality fifiltering on scRNASeq data was performed by multiple fifiltering parameters including >50% of mitochondrial genes, cells expressing the lower number of genes (<200 genes), and

       The unsupervised analysis using principal component analysis (PCA) was performed on variable genes to identify principal components, which captured the most variance across the samples. These principal components were used as an input for Uniform Manifold Approximation and Projection (UMAP) analysis89 to determine the overall relationship among the cells. Cells with similar transcriptome profifiles clustered together, and the clusters were subsequently annotated to different cell types based on the expression of specifific well-established cell marker transcripts. Comparative analysis of the single-cell landscape of healing and nonhealing DFUs, along with healthy non-DM subjects and non-DFU DM patients as controls, was performed using split UMAP plots, for determining heterogeneity (based on clusters of cells) and abundance of cell types. The signifificance testing change in abundance of cell types across clinical groups was performed either using one-way ANOVA or Welch’s t-test (p-value < 0.05). Similar analysis was also performed for the 3 different anatomical sites separately from where the samples were collected, i.e., foot, forearm, and peripheral blood. To further characterize cell type-specifific differences among clinical groups, we performed comparative analyses using multiple tests corrected non-parametric Wilcoxon Rank Sum test (P Adjusted value = 0.01, Fold Change = 1.2) on individual cell types like fifibroblasts, keratinocytes, T-lymphocytes, natural killer cells, monocytes, macrophages, mast cells, B-lymphocytes, plasma cells, and dendritic cells. genes only uniquely expressed in <3 cells.

       Pathways and systems biology analysis. To precisely characterize the cell types and understand the molecular mechanism of wound healing, we performed pathways enrichment and systems biology analysis. The analysis was performed on transcripts that were signifificantly dysregulated in the specifific cells by comparing healed vs non-healed samples. Pathways and systems biology analysis was performed using the Ingenuity Pathway Analysis software package (IPA 9.0) (Qiagen). A detailed description of IPA is available at the Ingenuity Systems’ website (http:// www.ingenuity.com). Systems biology analysis was performed by analyzing the upstream transcriptional regulators. The regulatory analysis helps in identifying signifificantly activated or inhibited transcriptional regulators based on upregulation or downregulation of its target genes. The signifificance of transcriptional regulators activation/inhibition was determined using one-tailed Fisher’s exact test. The regulators with a p-value <0.01 and absolute z-score 2 were considered statistically significant.

       Ligand and receptor-based cell interaction analysis. NicheNetR60 was used to identify ligands produced by Healer-specifific fifibroblasts, which could uniquely regulate other healer-specifific fifibroblasts. NicheNetR uses a prior model of ligandtarget interactions derived from a meta-analysis of multiple sources to identify ligands that may explain expression differences in a given gene set. In this workflflow, cells are classifified as either senders or receivers. The expression of sender cells is used to identify possible ligands, while the receiver cells are used to generate a gene set. In this case, Cluster 3 with overexpression of MMP1, MMP3, CHI3L1, CCL20, and TIMP1 from the Healer-specifific fifibroblasts was treated as a sender cluster, while other Healer-specifific clusters (Fig. 5; clusters 4, 6, 13) or non-specifific clusters (Fig. 5; 0, 2, 5) were treated as receivers. The gene sets used were the markers differentially expressed between DFU-Healers and DFU-Non-healer samples within the receiver subsets. The top markers are combined with NicheNetR’s ligand-target weights to compute the Pearson correlation coeffificient between ligands and expression changes in the receiver subset. A high Pearson correlation coeffificient between ligand and target gene set indicates that expression of ligand might be responsible for expression differences. For a ligand to be considered for interaction analysis, it must be expressed in at least 5% of the sender cell population, and its corresponding receptor must be expressed in 5% of the receiver cell population.

       Cellular trajectory and differentiation state analysis. To measure the transcriptional dynamics and characterize differentiation process based on single-cell data from DFU healers and non-healers, we performed RNA velocity analysis using the Velocyto90 and scVelo algorithms91. We performed RNA velocity analysis on HE-Fibro, Fibro, SMCs, M1-Macro, and M2-Macro cells in the DFU-Healer and DFU-Non-healer subsets of the foot cells. We generated spliced and unspliced counts using the Velocyto package and merged the data from all patients. This merged dataset was subsetted to retain only HE-Fibro, Fibro, SMCs, M1-Macro,and M2-Macro cells from foot samples as they showed signifificant association with the wound healing process. Velocity streams and inference of root cells were generated using scVelo version 0.2.391. This python package uses steady state and dynamic models to predict cell trajectory and latent time. The RNA velocities were projected onto a computed UMAP for each subset; the streamline velocity vectors represent directions and flflow of estimated trajectory and differentiation of the cells. The dynamic model also predicts the latent time of the cells, which represents the cell’s position in a biological process.

        Enrichment of HE-Fibro and Macrophages gene signatures in the temporal healing data from acute burn wound study. To assess the enrichment of HEFibro and M1 Macrophages in the sequential wound healing samples, we performed external validation using microarray data from a temporal acute burn wound study. The raw data downloaded from the ArrayExpress database accession number E-MTAB-1323 contain gene expression data of skin from 5 patients on 6 different time points, pre- and post-burn injury. The R packages beadarray92 and limma93 were used to pre-process, normalize, and analyze gene expression data. Gene set enrichment analysis (GSEA) was performed using the GSVA R package94 to compare enrichment of HE-Fibro and M1-Macro gene signatures. These signatures are subsets of differentially expressed genes comparing DFU-Healer and DFU-Non-healer subsets. The HE-Fibro gene signature used in GSEA analysis consists of PLA2GA, MMP1, CHI3L1, TIMP1, SFRP4, FTH1, FN1, MT2A, LUM, CHI3L2, MMP13, HIF1A, CCL20, TPM2, ASPN, MMP3, TNFAIP6, and IL6. The M1-Macro gene signature consists of IL1B, S100A8, VCAN, BCL2A1, LYZ, S100A9, TIMP1, C15orf48, SRGN, NFKBIA, BTG1, NAMPT, PLAUR, SAT1, ID2, TYMP, SLC2A3, SERPINA1, CXCL8, and SOD2.

        After calculating enrichment scores of the signatures at each time point for individual subjects, average enrichment scores and standard deviation of each time point were calculated. Further statistical signifificance of changes in the enrichment score was determined by one-way mixed-effects ANOVA with Bonferroni corrections.

         Immunoflfluorescence staining and imaging. For confifirmation of HE-Fibro presence, 5-μm-thick frozen sections from healing and non-healing DFUs were fifixed in 80% ice-cold acetone for 10 min, blocked with 5% donkey serum in 0.2% PBSTween for 30 min at room temperature, and incubated overnight in a humidifified chamber at 4 °C with primary antibodies: mouse monoclonal anti-FAP (1:50, clone F11-24, sc-65398, Santa Cruz Biotechnology), rabbit polyclonal anti-CHI3L1 (1:100, ab77528, Abcam) and goat polyclonal anti-TIMP1 (1:100, AF970, R&D Systems). Alexa Fluor donkey anti-rabbit 488-, anti-mouse 594- and anti-goat 647- conjugated secondary antibodies (1:1000 ab150061, 1:500 ab150112, and 1:1000 ab150131, respectively, all Abcam) were added the next day for 1 h at room temperature. 2-(4-amidinophenyl)-1H-indole-6-carboxamidine (DAPI) was included for nuclear counterstaining. TrueVIEW Autoflfluorescence Quenching Kit (Vector Labs, SP-8400) treatment was employed to enhance staining. Tissue sections were mounted in ProLong Gold Antifade (Thermo Fisher Scientifific, P36930) and visualized with a Zeiss LSM 880 (Carl Zeiss) inverted confocal microscope and images processed with ZEN 2011 (Carl Zeiss) and ImageJ/FIJI (NIH) software  packages.

         For validation of macrophages and most signifificant pathways, paraffifinembedded immunoflfluorescent staining of healing and non-healing DFUs was performed. 5-μm-thick sections were deparaffifinized, rehydrated and antigen retrieval was achieved with citrate buffer pH 6.0 in a pressure cooker for 20 min. The sections were then blocked with 5% donkey serum in 0.2% PBS-Tween for 1 h at room temperature and incubated overnight in a humidifified chamber at 4 °C with primary antibodies: goat polyclonal anti-IL17 (1:100, AF-317-NA, R&D Systems), goat polyclonal anti-TNFA (1:50, AF-410-NA, R&D Systems), rabbit monoclonal anti-HIF1A (1:100, clone EP1215Y, ab51608, Abcam), rabbit polyclonal anti-IL6 (1:100, ab6672, Abcam), goat polyclonal anti-IL1B (1:50, AF-201-NA, R&D Systems), goat polyclonal anti-S100A8 (1:100, AF3059, R&D Systems), mouse monoclonal anti-DAB2 (1:50, clone E-11, sc-136964, Santa Cruz Biotechnology), rabbit monoclonal anti-CD68 (1:100, clone EPR20545, ab213363, Abcam), mouse monoclonal anti-VIM (1:100, clone V9, MAB3400, Sigma), mouse monoclonal anti-CD163 (1:20, clone GHI/61, sc-20066, Santa Cruz Biotechnology).

        Appropriate admixtures of Alexa Fluor donkey anti-rabbit 488- (1:1000, ab150061) and 594- (1:500, ab150064), anti-mouse 488- (1:1000, ab150109) and 647- (1:1000, ab150107), and anti-goat 594- (1:500, ab150132) and 647- (1:1000, ab150131) conjugated secondary antibodies, all from Abcam, were added the next day for 1 h at room temperature. DAPI was included for nuclear staining. Sections were quenched for 5 min using the TrueView Autoflfluorescence Quenching kit to decrease background (Vector Laboratories) and covered with anti-fade mounting medium. Images were obtained at ×20 magnifification with an Axio Imager A2 upright microscope using Zen Blue edition software (Zeiss).

        CD68/Vimentin (VIM) were used as guide stains to capture images within the ulcers with similar abundance of positively stained cells. Quantifification was performed on ImageJ/FIJI by counting the number of double-positive CD68 cells with respective macrophage polarization markers (DAB2, S100A8, IL1B, and CD163) and dividing by the area of the tissue for normalization; counting the number of IL17 positive cells per area; computing the percentage of stained area for TNFA, IL6, and HIF1A. Two measurements were averaged per sample.

         Spatial transcriptomics. The spatial transcriptome profifiling was performed using NanoString’s GeoMx Digital Spatial profifiling platform on unfifixed frozen 5-μm tissue sections. Samples were processed as follows: (1) 10% neutral buffered formalin (NBF) fifixation overnight, (2) target retrieval (1X Tris EDTA, pH 9.0 for 20 min), (3) proteinase K digestion (1 µg/mL for 15 min), (4) post-fifixation (10% NBF, Tris-glycine stop buffer), (5) in situ hybridization overnight with the GeoMx Cancer Transcriptome Atlas probe panel (1800-plex), (6) stringent washes (50:50 formamide/4X SSC), and (7) flfluorescent antibody/marker (aSMA, 1:100, Clone: 1A4, Abcam; CD45, 1:100, Clone: 2B11 + PD7/26, Novus; PanCK, 1:50, Clone: AE1/AE3, Novus) incubation, 1 h at room temperature. Sections were then loaded onto the GeoMx® Digital Spatial Profifiler (Nanostring, GMX-DSP). For profifiling, circular regions of interest (ROIs), ~500 μm in diameter, located within the ulcers or in neighboring non-ulcerated tissue were selected to include high concentrations of CD45+ immune cells in close proximity to vessels (αSMA + structures). After ROI selection, the GeoMx instrument illuminated each ROI separately with UV light to cleave, aspirate, and deposit the oligonucleotides from the hybridized ISH probes for downstream sequencing into a 96-well plate. Library preparation (PCR, AMPure bead purifification) was performed, followed by paired-end sequencing with an Illumina NextSeq 550. Sequencing data (FASTQs) was then processed with a custom GeoMx NGS pipeline (DCCs) to be analyzed in part with the GeoMx Data Analysis Suite. Raw reads were processed for high quality with TrimGalore and FLASH95. Reads were then aligned to analyte barcode with Bowtie296. PCR duplicates were discarded using UMI-tools with the Hamming distance set at three. Poorly performing probes were removed from analysis if they were outliers (Grubbs test) or had low counts relative to other probes targeting the same gene. Raw probe count data (up to 5 unique probes per gene) were condensed into gene level count data and normalized with the quartile 3 gene count value per ROI individually. Complete-linkage hierarchical clustering was performed on normalized counts and represented by heatmap using the R function pheatmap. Unpaired t-test with Benjamini–Hochberg procedure for adjusted p-values was used to calculate differentially expressed genes with a threshold p < 0.05. Signifificantly expressed genes were entered on Metascape (Version 3.5, http://metascape.org) for enrichment analysis with Gene Ontology (GO) Biological Processes (Version 2020-09-16). All genes in the human genome were used as the enrichment background. P-values were calculated based on cumulative hypergeometric distribution and Q-values were calculated using the Benjamini–Hochberg procedure for multiple testing. A term was considered overrepresented when p < 0.01, had a minimum count of 3 and an enrichment factor >1.5, which is the ratio between the observed counts and the counts expected by chance. Volcano plot and GO bar graphs were designed with Prism 8.4.2 (GraphPad).

        Cell culture. Normal human dermal fifibroblast cells (BJ CRL-2522) were obtained from ATCC and maintained in Eagle’s minimum essential medium (EMEM) (ATCC, 30-2003), supplemented with 1% (v/v) penicillin/streptomycin (P/S) and 10% (v/v) fetal bovine serum (FBS) (Sigma-Aldrich, F1435). For passaging, cells at ~80% conflfluence were detached through a 5- to 10-min incubation with 0.05% Trypsin/EDTA and further resuspended in complete EMEM. Cells were then centrifuged at 1200 rpm for 5 min. The cells were replated at a concentration of 6000 cells/cm2 and/or cryopreserved with 90% FBS and 10% DMSO freezing media. Cells were maintained in 95% O2, 5% CO2 at 37 °C and routinely tested for mycoplasma contamination (PromoKine, PK-CA91-1096).

        Transduction of fifibroblasts with Precision LentiORF viral vectors. Cells were seeded in 6-well culture plate at 150,000 cells per well and pre-incubated with 5 µg/ml polybrene for 10 min at 37 °C. Afterward, cells were incubated overnight with  culture medium containing 5 µg/ml polybrene and the viral particles carrying the CHI3L1 gene (OHS5899-202624268, Horizon Discoveries) or the positive control viral particles (OHS5833) at a multiplicity of infection of 10. After removal of the particles containing medium, cells were incubated in culture medium with 10 µg/ml Blasticidin to positively select transduced cells. Transduction effificiency was evaluated with assessment of GFP expression, for both target and control constructs and RFP expression, for the control construct, on a K2 Cellometer (Nexcelom Bioscience) and with live-cell imaging on a Zeiss LSM 880 microscope.

        Real-time qPCR. RNA was extracted from 100,000 cells using the miRNeasy Mini Kit (Qiagen, 217004). RNA quantifification was done by using the Qubit RNA BR Assay kit (Cat. No. Q10210) and the Qubit 3 Fluorometer. cDNA was used at a concentration of 15 ng/ml from 1 μg of RNA and reverse transcribed with the miScript II RT kit (Cat No. 218161). RT-qPCR analysis was run for the samples using a miScript SYBR Green PCR Kit (Cat. No. 218073) on a Stratagene Mx3005P (Agilent Technologies). Housekeeping gene GAPDH primers were purchased from Qiagen (Cat. No. QT00079247) and CHI3L1 primers, were obtained from MGHHMS primer bank with the following sequences: FW: 5′-GAA GAC TCT CTT GTC TGT CGG A-3′ and RV: 5′-AAT GGC GGT ACT GAC TTG ATG-3′. Data were normalized to the expression of GAPDH and were analyzed using the 2−ΔΔCT  method.

        Western blotting. 500,000 cells were centrifuged at 130 × g for 5 min at 4 °C, washed with ice-cold PBS, and centrifuged again at 2400 × g for 5 min at 4 °C.

        The pellet was then resuspended in ice-cold RIPA buffer (Prod# 89901) supplemented with 10 µl/ml protease and phosphatase inhibitors (Prod# 78430 and 78420) and incubated for 15 min on ice with periodical pipetting and vortexing. Samples were then centrifuged at 14,000 × g for 15 min at 4 °C and supernatants were collected. Protein concentration was measured using the Pierce™ BCA Protein Assay Kit (Cat. No. 23225). The protein samples were reduced by using a 6x Laemmli buffer, and boiled at 95 °C for 5 min. 30 µg of protein per sample was loaded into 12% SDS-PAGE gels and run at a constant 200 V for 40 min. The gel was washed with Tris-buffered saline-Tween 20 (TBST), and incubated in blotting buffer for 10 min. The transfer ran overnight in a cold room at a constant 90 mA. Once transfer was complete, the blot was washed in TBST, and blocked with 5% BSA for 1 h at room temperature. Blot was then incubated in 5% BSA with CHI3L1 (Abcam, ab77528 1:1000) or GAPDH (ab9485, 1:5000) antibodies for 1 h at room temperature. Afterward, the blot was washed in TBST, and incubated with a secondary antibody (ab205718 1:10,000) for 1 h at room temperature. Finally, a chemiluminescent substrate (Cat. #1705062) was added and the blot was visualized using the ChemiDocTM Touch Imaging System (Bio-Rad). For stripping, the blot was washed in TBST and incubated in stripping buffer (Prod# 46430) for 45 min at room temperature. After stripping, the blot was washed with TBST, blocked with 5% BSA for 1 h at room temperature, and reprobed as previously described.

      Adhesion assay. Transduced BJ cells were plated at 50,000 cells per well in 12-well plates pre-coated with 10 µg/ml human fifibronectin (Prod# 33016-015) and incubated for 1 h at 37 °C. Afterward, cells were washed with PBS, fifixed with 4% paraformaldehyde for 15 min, and stained with 0.05% crystal violet for 30 min at room temperature. Pictures of the adherent cells were taken on a Primo Vert inverted microscope (Carl Zeiss) with an Axiocam 105 camera. Pictures of two random fifields were taken per well, for at least three wells per condition, and cells were counted using ImageJ/FIJI software. Three independent experiments were  performed.

      Scratch assay. Transduced BJ cells were plated at 50,000 cells per well in 24-well plates. The cell monolayer was scratched in a straight line using a 200 µl pipette tip. Debris were removed by washing once with media, then cells were incubated in medium supplemented with 5% FBS throughout the experiment. Images were taken immediately after the scratch, in 6 h, and in 12 h for four wells per condition. Scratch areas were analyzed using ImageJ/FIJI. Three independent experiments were performed.

     Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability

     Spatial transcriptomics and scRNASeq data have been submitted to NCBI’s Gene Expression Omnibus (GEO) and are accessible through GEO accession numbers GSE166120 and GSE165816. Burn wound gene expression data were downloaded from ArrayExpress accession number E-MTAB-1323. An interactive data resource and analytical tool developed based on this DFU single-cell data are available online at https://bhasinlab.bmi.emory.edu/Diacomp. Source data are provided with this paper.

Code availability

    Code for data analyses is described in the “Methods” and is available from the corresponding authors upon request. Received: 18 March 2021; Accepted: 29 November 2021;published online :10 January 2022

References

1.Jeffcoate, W. J. et al. Current challenges and opportunities in the prevention and management of diabetic foot ulcers. Diabetes Care 41, 645–652 (2018).

2.Armstrong, D. G. et al. Five year mortality and direct costs of care for people with diabetic foot complications are comparable to cancer. J. Foot Ankle Res. 13, 16 (2020).

3.Sen, C. K. Human wounds and its burden: an updated compendium of estimates. Adv. Wound Care 8, 39–48 (2019).

4.Armstrong, D. G., Boulton, A. J. M. & Bus, S. A. Diabetic foot ulcers and their recurrence. N. Engl. J. Med. 376, 2367–2375 (2017).

5.Baltzis, D., Eleftheriadou, I. & Veves, A. Pathogenesis and treatment of impaired wound healing in diabetes mellitus: new insights. Adv. Ther. 31, 817–836 (2014).

6.Stubbington, M. J. T. et al. Single-cell transcriptomics to explore the immune system in health and disease. Science 358, 58–63 (2017).

7.Regev, A. et al. The human cell atlas. Elife 6, e27041 (2017).

8.Theocharidis, G. et al. Integrated skin transcriptomics and serum multiplex assays reveal novel mechanisms of wound healing in diabetic foot ulcers. Diabetes 69, 2157–2169 (2020).

9.Merritt, C. R. et al. Multiplex digital spatial profifiling of proteins and RNA in fifixed tissue. Nat. Biotechnol. 38, 586–599 (2020).

10.He, H. et al. Single-cell transcriptome analysis of human skin identififies novel fifibroblast subpopulation and enrichment of immune subsets in atopic dermatitis. J. Allergy Clin. Immunol. 145, 1615–1628 (2020).

11.Reynolds, G. et al. Developmental cell programs are co-opted in inflflammatory skin disease. Science 371, eaba6500 (2021).

12.Stoeckius, M. et al. Simultaneous epitope and transcriptome measurement in single cells. Nat. Methods 14, 865–868 (2017).

13.Tellechea, A. et al. Mast cells regulate wound healing in diabetes. Diabetes 65, 2006–2019 (2016).

14.Arase, H. et al. Cutting edge: the mouse NK cell-associated antigen recognized by DX5 monoclonal antibody is CD49b (alpha 2 integrin, very late antigen-2). J. Immunol. 167, 1141–1144 (2001).

15.Dumin, J. A. et al. Pro-collagenase-1 (matrix metalloproteinase-1) binds the alpha(2)beta(1) integrin upon release from keratinocytes migrating on type I collagen. J. Biol. Chem. 276, 29368–29374 (2001).

16.Gagliani, N. et al. Coexpression of CD49b and LAG-3 identififies human and mouse T regulatory type 1 cells. Nat. Med. 19, 739–746 (2013).

17.Kazakova, M. H. & Sarafifian, V. S. YKL-40-a novel biomarker in clinical practice? Folia Med. 51, 5–14 (2009).

18.Sanggaard, K. W. et al. The TSG-6 and I alpha I interaction promotes a transesterifification cleaving the protein-glycosaminoglycan-protein (PGP) cross-link. J. Biol. Chem. 280, 11936–11942 (2005).

19.Milner, C. M. & Day, A. J. TSG-6: a multifunctional protein associated with inflfammation. J. Cell Sci. 116, 1863–1873 (2003).

20.Burgy, O. & Konigshoff, M. The WNT signaling pathways in wound healing and fifibrosis. Matrix Biol. 68-69, 67–80 (2018).

21.Whyte, J. L. et al. Augmenting endogenous Wnt signaling improves skin wound healing. PLoS ONE 8, e76883 (2013).

22.Vallee, A. et al. Interactions between TGF-beta1, canonical WNT/beta-catenin pathway and PPAR gamma in radiation-induced fifibrosis. Oncotarget 8, 90579–90604 (2017).

23.Vela-Romera, A. et al. Characterization of the human ridged and non-ridged skin: a comprehensive histological, histochemical and immunohistochemical analysis. Histochem Cell Biol. 151, 57–73 (2019).

24.Zhang, X., Yin, M. & Zhang, L. J. Keratin 6, 16 and 17-critical barrier alarmin molecules in skin wounds and psoriasis. Cells 8, 807 (2019).

25.Kerkhoff, C. et al. Novel insights into the role of S100A8/A9 in skin biology. Exp. Dermatol. 21, 822–826 (2012).

26.Kursar, M. et al. Differential requirements for the chemokine receptor CCR7 in T cell activation during Listeria monocytogenes infection. J. Exp. Med. 201, 1447–1457 (2005).

27.Hendriks, J. et al. CD27 is required for generation and long-term maintenance of T cell immunity. Nat. Immunol. 1, 433–440 (2000).

28.Mahnke, Y. D. et al. The who’s who of T-cell differentiation: human memory T-cell subsets. Eur. J. Immunol. 43, 2797–2809 (2013).

29.Zeglinski, M. R. & Granville, D. J. Granzymes in cardiovascular injury and disease. Cell Signal 76, 109804 (2020).

30.Filen, J. J. et al. Quantitative proteomics reveals GIMAP family proteins 1 and 4 to be differentially regulated during human T helper cell differentiation. Mol. Cell Proteom. 8, 32–44 (2009).

31.Bekeschus, S. et al. Distinct cytokine and chemokine patterns in chronic diabetic ulcers and acute wounds. Exp. Dermatol. 26, 145–147 (2017).

32.Spivia, W. et al. Complement protein C1q promotes macrophage antiinflflammatory M2-like polarization during the clearance of atherogenic lipoproteins. Inflflamm. Res. 63, 885–893 (2014).

33.Mendes, K. L., Lelis, D. F. & Santos, S. H. S. Nuclear sirtuins and inflflammatory signaling pathways. Cytokine Growth Factor Rev. 38, 98–105 (2017).

34.Advedissian, T. et al. E-cadherin dynamics is regulated by galectin-7 at epithelial cell surface. Sci. Rep. 7, 17086 (2017).

35.Janson, D. G. et al. Different gene expression patterns in human papillary and reticular fifibroblasts. J. Invest. Dermatol. 132, 2565–2572 (2012).

36.Korosec, A. et al. Lineage identity and location within the dermis determine the function of papillary and reticular fifibroblasts in human skin. J. Invest. Dermatol. 139, 342–351 (2019).

37.Haydont, V. et al. Genome-wide profifiling of adult human papillary and reticular fifibroblasts identififies ACAN, Col XI alpha1, and PSG1 as general biomarkers of dermis ageing, and KANK4 as an exemplary effector of papillary fifibroblast ageing, related to contractility. Mech. Ageing Dev. 177, 157–181 (2019).

38.Nauroy, P. et al. Human dermal fifibroblast subpopulations display distinct gene signatures related to cell behaviors and matrisome. J. Invest. Dermatol. 137, 1787–1789 (2017).

39.Ascension, A. M. et al. Human dermal fifibroblast subpopulations are conserved across single-cell RNA sequencing studies. J. Invest. Dermatol. 141, 1735–1744 (2021).

40.Sole-Boldo, L. et al. Single-cell transcriptomes of the human skin reveal age related loss of fifibroblast priming. Commun. Biol. 3, 188 (2020).

41.Tabib, T. et al. SFRP2/DPP4 and FMO1/LSP1 defifine major fifibroblast populations in human skin. J. Invest. Dermatol. 138, 802–810 (2018).

42.Dobaczewski, M., Gonzalez-Quesada, C. & Frangogiannis, N. G. The extracellular matrix as a modulator of the inflflammatory and reparative response following myocardial infarction. J. Mol. Cell Cardiol. 48, 504–511 (2010).

43.Gillan, L. et al. Periostin secreted by epithelial ovarian carcinoma is a ligand for alpha(V)beta(3) and alpha(V)beta(5) integrins and promotes cell motility. Cancer Res. 62, 5358–5364 (2002).

44.Yoshihara, T. et al. Periostin plays a critical role in the cell cycle in lung fibroblasts. Respir. Res. 21, 38 (2020).

45.Izuhara, K. et al. Periostin in inflflammation and allergy. Cell Mol. Life Sci. 74, 4293–4303 (2017).

46.Yamaba, S. et al. PLAP-1/Asporin regulates TLR2- and TLR4-induced inflflammatory responses. J. Dent. Res. 94, 1706–1714 (2015).

47.Portou, M. J. et al. Hyperglycaemia and ischaemia impair wound healing via toll-like receptor 4 pathway activation in vitro and in an experimental murine model Eur. J. Vasc. Endovasc. Surg. 59, 117–127 (2020).

48.Wang, R. et al. Inflflammatory-sensitive CHI3L1 protects nucleus pulposus via AKT3 signaling during intervertebral disc degeneration. FASEB J. 34, 3554–3569 (2020).

49.Zhao, T. et al. Chitinase-3 like-protein-1 function and its role in diseases.Signal Transduct. Target Ther. 5, 201 (2020).

50.Zhou, Y. et al. Chitinase 3-like 1 suppresses injury and promotes fifibroproliferative responses in Mammalian lung fifibrosis. Sci. Transl. Med. 6, 240ra76 (2014).

51.Deutschmann, C., Roggenbuck, D. & Schierack The loss of tolerance to CHI3L1—a putative role in inflflammatory bowel disease? Clin. Immunol. 199, 12–17 (2019).

52.Dallas, A. et al. Acceleration of diabetic wound healing with PHD2- and miR-210-targeting oligonucleotides. Tissue Eng. Part A 25, 44–54 (2019).

53.Niu, X. et al. Hypoxia regulates angeogenic-osteogenic coupling process via up-regulating IL-6 and IL-8 in human osteoblastic cells through hypoxiainducible factor-1alpha pathway. Cytokine 113, 117–127(2019).

54.Gallucci, R. M. et al. Interleukin-6 treatment augments cutaneous wound healing in immunosuppressed mice. J. Interferon Cytokine Res. 21, 603–609 (2001).

55.Lin, Z. Q. et al. Essential involvement of IL-6 in the skin wound-healing process as evidenced by delayed wound healing in IL-6-defificient mice. J. Leukoc. Biol. 73, 713–721 (2003).

56.Wang, X. P. et al. The interleukin-6 cytokine system regulates epidermal permeability barrier homeostasis. J. Invest. Dermatol. 123, 124–131 (2004).

57.Pradhan, L. et al. Gene expression of pro-inflflammatory cytokines and neuropeptides in diabetic wound healing. J. Surg. Res. 167, 336–342 (2011).

58.Barrientos, S. et al. Growth factors and cytokines in wound healing. Wound Repair Regen. 16, 585–601 (2008).

59.Lu, Y. E. & Chen, Y. J. Resveratrol inhibits matrix metalloproteinase-1 and -3 expression by suppressing of p300/NFkappaB acetylation in TNF-alphatreated human dermal fifibroblasts. Chem. Biol. Interact. 337, 109395 (2021).

60.Browaeys, R., Saelens, W. & Saeys, Y. NicheNet: modeling intercellular communication by linking ligands to target genes. Nat. Methods 17, 159–162 (2020).

61.Greco, J. A. 3rd et al. A microarray analysis of temporal gene expression profifiles in thermally injured human skin. Burns 36, 192–204 (2010).

62.Driskell, R. R. et al. Distinct fifibroblast lineages determine dermal architecture in skin development and repair. Nature 504, 277–281 (2013).

63.Guerrero-Juarez, C. F. et al. Single-cell analysis reveals fifibroblast heterogeneity and myeloid-derived adipocyte progenitors in murine skin wounds. Nat. commun. 10, 650 (2019).

64.Lynch, M. D. & Watt, F. M. Fibroblast heterogeneity: implications for human diseace. J. Clin. Invest 128, 26–35 (2018).

65.Mascharak, S., desJardins-Park, H. E. & Longaker, M. T. Fibroblast heterogeneity in wound healing: hurdles to clinical translation. Trends Mol. Med. 26, 1101–1106 (2020).

66.Philippeos, C. et al. Spatial and Single-Cell Transcriptional Profifiling Identififies Functionally Distinct Human Dermal Fibroblast Subpopulations. J. Invest. Dermatol. 138, 811–825 (2018).

67.Nayar, S. et al. Immunofifibroblasts are pivotal drivers of tertiary lymphoid structure formation and local pathology. Proc. Natl Acad. Sci. USA 116, 13490–13497 (2019).

68.Sinha, M. et al. Direct conversion of injury-site myeloid cells to fifibroblast-like cells of granulation tissue. Nat. Commun. 9, 936 (2018).

69.Donlin, L. T. et al. Modulation of TNF-induced macrophage polarization by synovial fifibroblasts. J. Immunol. 193, 2373–2383 (2014).

70.Humeres, C. et al. Cardiac fifibroblast cytokine profifiles induced by proinflflammatory or profifibrotic stimuli promote monocyte recruitment and modulate macrophage M1/M2 balance in vitro. J. Mol. Cell Cardiol. 101, 69–80 (2016).

71.Ploeger, D. T. et al. Cell plasticity in wound healing: paracrine factors of M1/ M2 polarized macrophages inflfluence the phenotypical state of dermal  fibroblasts. Cell Commun. Signal 11, 29 (2013).

72.Spiller, K. L. & Koh, T. J. Macrophage-based therapeutic strategies in regenerative medicine. Adv. Drug Deliv. Rev. 122, 74–83 (2017).

73.Jetten, N. et al. Wound administration of M2-polarized macrophages does not improve murine cutaneous healing responses. PLoS ONE 9, e102994 (2014).

74.Dreymueller, D. et al. Embryonic stem cell-derived M2-like macrophages delay cutaneous wound healing. Wound Repair Regen. 21, 44–54 (2013).

75.Tamoutounour, S. et al. Origins and functional specialization of macrophages and of conventional and monocyte-derived dendritic cells in mouse skin. Immunity 39, 925–938 (2013).

76.Wynn, T. A. & Vannella, K. M. Macrophages in tissue repair, regeneration, and fifibrosis. Immunity 44, 450–462 (2016).

77.Bajpai, G. et al. The human heart contains distinct macrophage subsets with divergent origins and functions. Nat. Med. 24, 1234–1245 (2018).

78.Sawaya, A. P. et al. Deregulated immune cell recruitment orchestrated by FOXM1 impairs human diabetic wound healing. Nat. Commun. 11, 4678 (2020).

79.Moura, J. et al. Impaired T-cell differentiation in diabetic foot ulceration. Cell Mol. Immunol. 14, 758–769 (2017).

80.Stone, R. C. et al. A bioengineered living cell construct activates an acute wound healing response in venous leg ulcers. Sci. Transl. Med. 9, eaaf8611(2017).

81.Tellechea, A. et al. Topical application of a mast cell stabilizer improves impaired diabetic wound healing. J. Invest. Dermatol. 140, 901–911 e11 (2020).

82.Nassiri, S. et al. Relative expression of proinflflammatory and antiinflflammatory genes reveals differences between healing and nonhealing human chronic diabetic foot ulcers. J. Invest. Dermatol. 135, 1700–1703 (2015).

83.Januszyk, M. et al. Characterization of diabetic and non-diabetic foot ulcers using single-cell RNA-sequencing. Micromachines 11, 815 (2020).

84.Stojadinovic, O. et al. Quality assessment of tissue specimens for studies of diabetic foot ulcers. Exp. Dermatol. 22, 216–218 (2013).

85.Lake, B. B. et al. A comparative strategy for single-nucleus and single-cell transcriptomes confifirms accuracy in predicted cell-type expression from nuclear RNA. Sci. Rep. 7, 6031 (2017).

86.Ding, J. et al. Systematic comparison of single-cell and single-nucleus RNA sequencing methods. Nat. Biotechnol. 38, 737–746 (2020).

87.Tan, Y. S. & Lei, Y. L. Isolation of tumor-infifiltrating lymphocytes by fificollpaque density gradient centrifugation. Methods Mol. Biol. 1960, 93–99 (2019).

88.Butler, A. et al. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol. 36, 411–420 (2018).

89.Maaten, L.J.P.v.d. &. Hinton, G. E. Visualizing high-dimensional data using t-SNE. J. Mach. Learn. Res. 9, 2579–2605 (2008).

90.La Manno, G. et al. RNA velocity of single cells. Nature 560, 494–498 (2018).

91.Bergen, V. et al. Generalizing RNA velocity to transient cell states through dynamical modeling. Nat. Biotechnol. 38, 1408–1414 (2020).

92.Dunning, M. J. et al. beadarray: R classes and methods for Illumina beadbased data. Bioinformatics 23, 2183–2184 (2007).

93.Ritchie, M. E. et al. Limma powers differential expression analyses for RNAsequencing and microarray studies. Nucleic Acids Res. 43, e47 (2015).

94.Hanzelmann, S., Castelo, R. & Guinney, J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinforma. 14, 7 (2013).

95.Magoc, T. & Salzberg, S. L. FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27, 2957–2963 (2011).

96.Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359 (2012).

Acknowledgements

      This work was supported by the NIDDK-sponsored Diabetic Complications Consortium grant 5U24DK115255-04 (A.V. and M.B.). A.V. received funding from the National Rongxiang Xu Foundation. G.T. received a George and Marie Vergottis Foundation Postdoctoral Fellowship award.

Author contributions

      V.H., A.V., and M.B. obtained funding. H.C.H., V.H., A.V., and M.B. supervised the project. G.T., A.V., and M.B. designed the experiments. G.T., B.E.T., D.S., S.S.B., T.S.S., R.F.S., I.M., P.W., and A.L. conducted the experiments. G.T., B.E.T., D.S., H.L.M., W.J.R.P., B.D., W.P., A.K., I.V., S.S.B., A.V., and M.B. analyzed data. G.T., B.E.T., D.S., S.S.B., H.L.M., W.J.R.P., A.V., and M.B. wrote the manuscript.