Detailed gene activity patterns show unique molecular signs and disrupted processes in different tissues of human diabetic foot ulcers

A total of 52 subjects were recruited for this study, including 26 patients diagnosed with diabetic foot ulcers (DFU) and 26 control subjects (CTRL). The ethical approval for this study was obtained from the Ethical Review Board at the Third People's Hospital of Yunnan Province (approval number: 2023KY046), and all procedures conformed to the principles outlined in the Declaration of Helsinki. Written informed consent was obtained from each participant before sample collection.

Tissue samples of the experimental group were collected below the knee joints, specifically targeting skin, fat, and muscle tissues surrounding active wounds. For the CTRL group, tissue samples were obtained from patients undergoing surgical treatment for various orthopedic conditions. Biopsies were collected from the following sites: right/left femoral neck fracture, right tibiofibular fracture after surgery, left femur fracture after surgery, right lateral epidermal osteoma and osteoarthritis, comminuted fracture of the proximal left ulna, bilateral femoral head ischemic necrosis, severe osteoarthritis of the left knee, developmental dysplasia of the right hip with total dislocation, left femoral head ischemic necrosis, right femoral head necrosis, and comminuted fracture. Before tissue extraction, rigorous clinical preparation involving disinfection, deiodination, and debridement was performed, preserving the complete structural integrity of each wound. Tissue specimens measuring approximately 0.5 cm² each were surgically excised from defined regions adjacent to the wounds. Collected tissues were immediately placed into labeled, sterile, transparent centrifuge tubes, separately containing skin, adipose, and muscle tissues. All samples were snap-frozen and stored at −80 °C until subsequent RNA extraction and analysis.

Total RNA was extracted from the samples using TRIzol reagent (Invitrogen, Thermo Fisher Scientific), following the manufacturer's protocol. RNA quality and integrity were evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies), and only samples with RNA integrity numbers (RIN) above 7 were selected for library preparation. Poly(A)-enriched RNA sequencing (RNA-seq) libraries were prepared using standard Illumina protocols involving poly(A)-mRNA isolation, fragmentation, first and second-strand cDNA synthesis, adaptor ligation, and PCR amplification.

Sequencing was performed on an Illumina NovaSeq X Plus platform to generate paired-end reads of 150 bp. Raw sequencing data were subjected to quality control analyses, and adaptor sequences were removed using trim_galore (version 0.6.7). High-quality clean reads obtained after trimming were used for subsequent bioinformatics analyses.

Clean reads obtained after trimming were aligned to the human reference genome (GRCh38) using the STAR aligner (28). Default parameters were employed with minor adjustments to optimize the alignment quality. SAM files generated by STAR were converted into sorted and indexed BAM files using SAMtools (version 1.18; https://github.com/samtools↗).

Gene expression quantification was performed by counting uniquely mapped reads to annotated gene regions with featureCounts (version 2.0.6; https://subread.sourceforge.net/featureCounts.html↗). Read counts for each gene were used as the input for subsequent differential expression analyses.

Raw gene counts were preprocessed and normalized following the integrated Differential Expression and Pathway (iDEP) analysis pipeline. Initial preprocessing involved filtering out low-expression genes based on the criteria of at least 0.5 counts per million (CPM) in a minimum of one sample. Gene identifiers were standardized using the gene ID conversion functionality of the iDEP pipeline (29, 30).

Normalization of gene expression data was carried out using the counts per million (CPM) method implemented by the EdgeR package (31). Subsequently, normalized counts were transformed using a balanced log transformation to reduce heteroscedasticity and improve the data's suitability for downstream analysis. Specifically, gene counts were log-transformed using the formula log2(CPM+c), with the constant (c=4) recommended by the iDEP server. The effectiveness of normalization and transformation was assessed by visual inspection of distribution plots, including box plots and density plots, ensuring comparable distributions across samples within each tissue type.

Differentially expressed genes (DEGs) between DFU and CTRL tissues within each specific tissue type (skin, fat, and muscle) were identified using the DESeq2 R package (32). DESeq2 employs a negative binomial generalized linear model to estimate variance-mean dependence in count data, providing robust statistical inference. Genes were classified as significantly differentially expressed based on adjusted p-values of <0.05 and fold-change criteria of ≥2 (absolute log2 fold change ≥1).

Hierarchical clustering analyses of DEGs were performed using the heatmap.2 function within the R environment to visualize differential expression patterns and sample grouping clearly. Volcano plots depicting the distribution of significant DEGs according to their statistical significance (-log₁₀ adjusted p-value) and expression fold-change (log2) were generated using the ggplot2 R package. Genomic distribution plots of DEGs across chromosomes were also generated to illustrate genome-wide alterations in gene expression profiles. Certain genes with high expression were validated through qPCR. The relative mRNA expression levels were calculated using the 2^-ΔΔCt method, with GAPDH used as the internal control.

To identify significantly enriched functional categories and biological pathways associated with DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted using ShinyGO (33). DEGs identified from each tissue were analyzed separately to elucidate tissue-specific functional implications. Statistical significance for enriched terms was determined using false discovery rate (FDR)-corrected p-values, with a significance threshold of FDR ≤0.05.

Furthermore, significantly enriched KEGG pathways identified from DEGs were visualized, and gene mapping was performed to illustrate the involvement of specific differentially regulated genes within each pathway. Genes significantly enriched in each KEGG pathway were annotated clearly to distinguish upregulated (red) and downregulated (green) genes, facilitating the identification of key regulatory genes implicated in DFU pathogenesis.

Proteins were extracted from human skin, adipose, and muscle tissue specimens by rinsing with phosphate-buffered saline (PBS) and homogenizing in ice-cold RIPA buffer (BP115DG, Boston BioProducts, Ashland, MA, USA) supplemented with EDTA and EGTA using a tissue homogenizer (KZ-III-FP, Servicebio, Wuhan, China). Homogenates were kept on ice for 30 min and centrifuged at 13,000 rpm for 20 min at 4 °C to remove debris. Protein concentrations in the supernatants were determined with a bicinchoninic acid (BCA) protein assay kit (P0010, Beyotime, Shanghai, China). Equal amounts of protein were mixed with 4× Laemmli sample buffer (1610747, Bio-Rad, Hercules, CA, USA) containing 5% β-mercaptoethanol and denatured at 95 °C for 5 min. Proteins were then separated on 10% SDS-PAGE gels and transferred to 0.22-μm polyvinylidene fluoride (PVDF) membranes (Millipore, Bedford, MA, USA) using a Bio-Rad wet transfer system. Membranes were blocked with 5% (w/v) skim milk prepared in Tris-buffered saline containing 0.1% Tween-20 (TBS-T) for 2 h at room temperature, followed by incubation with primary antibodies Col1a1 (HY-P81227, MedChemExpress, China), MMP9 (24317T, Cell Signaling Technology, China), PPARγ (YP-mAb-03325, Uping Bio, China), and IL17 (YP-mAb-12431, Uping Bio, China) overnight at 4 °C. After several washes, membranes were exposed to species-appropriate horseradish peroxidase (HRP)-conjugated secondary antibodies for 2 h at room temperature. Bands were visualized with the GeneGnome XRQ-NPC System (Synoptics, Frederick, MD, USA) or SCG-W3000 Plus (Servicebio, Wuhan, China) and quantified using ImageJ software (NIH, Bethesda, MD, USA).

Computational deconvolution of bulk RNA-seq profiles was performed to estimate the relative abundance of major cell types in skin, adipose, and muscle samples. Transcripts per million (TPM)–normalized data were input into CIBERSORTx (34) using a custom reference matrix derived from publicly available single-cell RNA-seq atlases of human skin, adipose, and skeletal muscle. Default parameters were applied with 100 permutations for significance testing, and cell types with deconvolution p <0.05 were considered reliably estimated. In parallel, cellular enrichment scores were computed using xCell (35) with the xCell64 signature gene set to provide complementary estimates of tissue-resolved cellular composition.

To infer paracrine ligand–receptor interactions across tissue compartments, we employed a two-step strategy. First, lists of differentially upregulated genes were generated for each tissue using our established differential expression pipeline (DESeq2; adjusted p <0.05, log2FC >1). Second, these gene lists were queried against curated ligand–receptor databases CellChat (36) and CellPhoneDB (37) to identify candidate interaction pairs. To prioritize biologically relevant regulatory links, we further applied NicheNet (38), which computes ligand–target regulatory potential scores. For cross-tissue inference, we considered ligands upregulated in one tissue paired with corresponding receptors either upregulated or highly expressed in a second tissue. Only ligand–receptor pairs supported by both databases or those with high NicheNet regulatory scores were retained for downstream analysis. Predicted inter-tissue interactions were visualized using chord diagrams and heatmaps to highlight signaling directionality, functional categories, and tissue-specific paracrine communication patterns.

To assess tissue-specific molecular signatures, we collected skin, fat, and muscle samples from patients diagnosed with DFU. For comparison, the same tissues were obtained from age-matched patients with lower limb fractures, serving as the nDFU control group. The demographic and clinical characteristics of both groups are summarized in Table 1. Quality control and exploratory analyses were conducted on RNA-seq data from DFU samples and CRTL tissues derived from skin, fat, and muscle (Supplementary Figures S1–S3). RNA-seq data were preprocessed by filtering low-expression genes, normalization, and log2 transformation. The boxplots and density plots revealed that the distributions of transformed gene expression values were comparable across the DFU and CTRL groups in skin (Supplementary Figures S1A, B), fat (Supplementary Figures S2A, B), and muscle tissues (Supplementary Figures S3A, B), indicating effective normalization. Furthermore, sequencing reads predominantly mapped to protein-coding genes across all tissue types, with smaller fractions mapping to immunoglobulin (IG), long non-coding RNA (lncRNA), and mitochondrial rRNA (Mt_rRNA) genes (Supplementary Figure S1C, S2C, S3C). Principal component analysis (PCA) plots showed clear clustering of samples by condition in each tissue type, with principal component 1 (PC1) accounting for a substantial proportion of the variance: 39% in skin (Supplementary Figure S1D), 45.5% in fat (Supplementary Figure S2D), and notably, 73.2% in muscle tissue (Supplementary Figure S3D). These distinct clusters reflect substantial differences in transcriptomic profiles between DFU and control samples within each tissue type, especially pronounced in muscle samples. Overall, the results from quality control analyses indicate that the data preprocessing steps were successfully applied, ensuring robust subsequent differential gene expression and enrichment analyses.

The transcriptomic profile of DFU skin displayed marked divergence from that of CTRL tissue. Hierarchical clustering of global gene expression profiles revealed two sharply distinct branches that aligned perfectly with the clinical classifications, highlighting the presence of strong, disease-specific transcriptional signatures (Figure 1A). Differential expression analysis identified 5,210 transcripts that met rigorous statistical thresholds (log2 FC >1, FDR <0.05). Notably, the majority of these genes were downregulated in DFU (n = 3,128), while 2,082 were upregulated, suggesting that transcriptional suppression is the dominant regulatory trend in ulcerated skin (Figure 1B; Supplementary Table S1). The range of deregulation was broad, spanning nearly three orders of magnitude, as visualized in the volcano plot (Figure 1C). Highly induced genes included stress-responsive epidermal markers such as SPRR2A, KRT16, and PI3, while transcripts linked to mitochondrial and metabolic functions were among the most strongly repressed. Chromosomal mapping of Chromosomal mapping of DEGs revealed a genome-wide distribution without discernible clustering, implying that the expression shifts stem from broad regulatory disruptions rather than localized genomic alterations (Figure 1D).

Pathway enrichment analysis further uncovered two overarching biological patterns. Upregulated genes in DFU skin were significantly enriched in pro-inflammatory and proliferative KEGG pathways, including cell cycle, IL17 signaling, cytokine–receptor interaction, p53 signaling, and cellular senescence—consistent with the hyperproliferative and inflamed phenotype observed in chronic, non-healing ulcers (Figure 2A). GO analysis reinforced this observation, revealing strong enrichment for biological processes such as epidermal development, keratinocyte differentiation, keratinization, and tissue remodeling (Figure 2C). Conversely, downregulated genes were predominantly involved in maintaining metabolic stability and structural integrity. KEGG pathways such as PPAR signaling, cell adhesion molecules, and cAMP signaling were notably depleted (Figure 2B). GO terms related to neurovascular development, circulatory system function, and fatty acid metabolism were also significantly under-represented (Figure 2D), indicating compromised neurovascular support and metabolic imbalance within the DFU microenvironment. Together, these findings outline a dual transcriptional program in DFU skin, characterized by widespread repression of genes associated with lipid metabolism, cell adhesion, and neurovascular signaling, alongside strong induction of pathways governing inflammation, proliferation, and epidermal stress responses. This transcriptomic polarity may underlie the pathological features of chronic DFU—persistent inflammation, hyperkeratinized wound edges, and impaired tissue repair—and identifies distinct molecular circuits that may serve as targets for future therapeutic intervention.

Transcriptomic profiling of DFU versus CTRL fat tissues revealed substantial differential gene expression patterns. Hierarchical clustering distinctly differentiated DFU and CTRL samples, reflecting marked transcriptional changes associated with DFU pathology (Figure 3A). Differential expression analysis (log2 fold change >1, adjusted p <0.05) identified 4,630 significantly dysregulated transcripts—1,992 of which were upregulated and 2,638 downregulated in DFU adipose tissue (Figure 3B; Supplementary Table S2). The volcano plot (Figure 3C) showed effect sizes spanning nearly ten-fold in both directions. Immune-related genes, including MMP9, CXCL5, and IL1RN, were among the most strongly upregulated, while key metabolic regulators involved in adipocyte function, such as PPARGC1A, PLIN1, and ADIPOQ, were prominently downregulated. Chromosomal distribution of these DEGs revealed a broad genomic dispersion across all autosomes and the X chromosome, with no evident clustering, suggesting a global reprogramming of transcription rather than localized genomic alterations (Figure 3D).

Pathway enrichment analysis supported this dichotomous expression landscape. KEGG analysis of the upregulated genes highlighted central inflammatory pathways, including cytokine–cytokine receptor interaction, IL17 signaling, chemokine signaling, and neutrophil extracellular trap (NET) formation, indicating heightened immune activity within the ulcer environment (Figure 4A). Consistently, GO terms such as "immune system process", "leukocyte activation", "adaptive immune response", and "cytokine production" were significantly overrepresented (Figure 4C), pointing to robust recruitment and activation of both innate and adaptive immune cells. In stark contrast, downregulated genes were enriched for pathways governing energy metabolism and lipid handling. KEGG terms included overarching metabolic processes as well as specific pathways like PPAR signaling, adipocyte lipolysis regulation, and degradation of branched-chain amino acids such as valine, leucine, and isoleucine (Figure 4B). GO analysis further supported this metabolic repression, with significant depletion of terms associated with fatty acid metabolism, monocarboxylic acid catabolism, and small molecule processing (Figure 4D). These findings describe an adipose environment where metabolic and lipid storage functions are transcriptionally silenced, while inflammatory and immune cell recruitment pathways are strongly activated. Together, these results illustrate a DFU adipose transcriptome defined by the co-existence of metabolic suppression and immunological overdrive. This dual state likely contributes to the chronic wound phenotype: impaired energy metabolism may hinder tissue repair, while sustained leukocyte infiltration and cytokine signaling exacerbate inflammation and tissue damage. These insights highlight both inflammatory pathways and adipocyte metabolic regulators as potential therapeutic targets for restoring homeostasis in DFU-affected tissues.

To characterize molecular alterations in DFU muscle, we performed RNA-seq on CTRL and DFU patients. Hierarchical clustering of the top DEGs revealed clear separation of DFU and CTRL samples (Figure 5A), with reciprocal patterns of up- and down-regulation across groups. In total, 9,923 genes met significance thresholds (log2 FC >1, adjusted p <0.05), of which 5,317 were upregulated and 4,606 downregulated in DFU muscle relative to controls (Figure 5B; Supplementary Table S3). A volcano plot illustrated the magnitude and significance of these changes (Figure 5C). Notably, immune-related mediators such as CXCL5, MMP9 and CCL20 dominated the positive fold-change tail, whereas canonical myofibrillar and oxidative metabolism genes such as ACTN2, MLIP, CAP2, and numerous electron-transport components occupied the negative extreme (Figure 5C). Genome-wide mapping of DEGs further underscored a widespread transcriptional reprogramming in DFU muscle, with up- and down-regulated genes distributed across all autosomes and the mitochondrial chromosome (Figure 5D).

To determine the biological processes underlying these expression changes, we conducted KEGG and GO enrichment analyses separately for up- and down-regulated gene sets (Figure 6). Among upregulated genes, the highest enrichment was seen for cytokine–cytokine receptor interaction, viral protein interaction with cytokine receptors, and chemokine signaling pathways (Figure 6A). Complementary GO analysis highlighted robust activation of immune and inflammatory processes, including "immune response", "leukocyte activation", and "inflammatory response" (Figure 6C), reflecting a potent pro-inflammatory milieu in DFU muscle.

Conversely, downregulated genes were strongly enriched for pathways related to oxidative phosphorylation, Parkinson's and Alzheimer's diseases (which share mitochondrial dysfunction signatures), and thermogenesis (Figure 6B). GO terms such as "oxidative phosphorylation", "ATP synthesis coupled electron transport", and "respiratory electron transport chain" were among the most significantly depleted biological processes (Figure 6D). These findings indicate that the DFU muscle exhibits coordinated downregulation of genes governing mitochondrial bioenergetics, concomitant with upregulation of inflammatory networks. Collectively, these data reveal a dual transcriptional phenotype in DFU muscle characterized by heightened inflammation and immune activation alongside suppression of mitochondrial energy metabolism, which may contribute to impaired muscle function and wound healing in diabetic foot ulcers.

Transcriptomic profiling of DFU tissues revealed both shared and tissue-specific alterations across skin, fat, and muscle compartments (Supplementary Tables S1–S3; Additional Files 1–3), underscoring the complex and multifactorial nature of impaired wound healing in diabetes. Common themes across all three tissues included upregulation of immune and inflammatory pathways and downregulation of genes involved in metabolic homeostasis, mitochondrial function, and structural integrity. Some highly expressed genes were also validated through qPCR (Figure 7), further confirming the findings from the RNA-seq analysis. Western blot analysis (Figure 8) of skin, adipose, and muscle tissues revealed that COL1A1, MMP, and IL17 proteins were consistently upregulated in DFU samples compared with non-DFU controls, whereas PPARγ protein expression was reduced across all tissue types, with the most pronounced decrease observed in adipose tissue and skin. Together, these findings validate the protein-level confirmation of RNA-seq results.

Computational deconvolution of bulk RNA-seq profiles revealed distinct alterations in cellular composition across DFU and control tissues (Figure 9). Major cell type distributions across samples (Figure 9B) indicated that skin (CTRL skin; CS vs. diabetic foot skin;DFS), adipose (CTRL fat; CF vs. diabetic foot fat; DFF), and muscle (CTRL muscle; CM vs. diabetic foot muscle; DFM) each exhibited distinct compositional shifts. Detailed resolution of immune, stromal, and parenchymal subtypes (Figure 9C) highlighted that DFS was enriched in keratinocytes, endothelial cells, and diverse immune subsets compared to CS, reflecting proliferative and inflammatory remodeling of the epidermal barrier. DFF showed expansion of stromal fibroblasts, macrophages, and mesenchymal stromal cells alongside diminished adipocyte fractions relative to CF, consistent with adipocyte dedifferentiation and pro-fibrotic remodeling. DFM samples exhibited loss of mature myofiber signatures and increased representation of progenitor and immune lineages compared to CM, suggesting myofiber attrition with stromal and immune infiltration. Boxplot analyses of representative cell classes (Figure 9D) corroborated these trends: higher epithelial and immune scores in DFS, elevated fibroblast and stromal signatures in DFF, and reduced contractile cell fractions with increased progenitor compartments in DFM. Collectively, these findings demonstrate that DFU tissues undergo compartment-specific shifts characterized by keratinocyte hyperplasia in skin, adipocyte depletion in adipose, and myofiber loss in muscle, each accompanied by expansion of stromal and immune components.

To explore how these compositional imbalances translate into inter-tissue communication, we inferred ligand–receptor interactions across skin, adipose, and muscle. A chord diagram highlighted extensive bidirectional signaling among all three compartments (Figure 9E), while a Sankey plot summarized the functional distribution of paracrine interactions across inflammatory, fibrotic/ECM, angiogenic, metabolic/lysosomal, and immune-modulatory categories (Figure 9F). Detailed analyses of direction-specific ligand–receptor axes are presented in Supplementary Figures S4–S9. To further probe inter-compartmental communication, we constructed tissue-resolved ligand–receptor expression maps by integrating differential expression with curated ligand–receptor databases and NicheNet prioritization. This analysis revealed a dense paracrine signaling network linking skin, adipose, and muscle, in which inflammatory, fibrotic, angiogenic, and metabolic cues were exchanged bidirectionally. Skin-derived ligands prominently included chemokines (CXCL5), cytokines (OSM, TNFSF13B, IL36G), and structural components (COL6A1/3, laminins, POSTN) (Supplementary Figures S4, S6). These signals projected into adipose and muscle receptors linked to immune activation, matrix remodeling, and stress pathways, indicating that inflamed skin actively drives inflammation, fibrosis, and metabolic dysfunction in deeper tissues. In turn, adipose and muscle secreted angiogenic and fibrotic ligands (ANGPT1/4, TGFB3, VEGFC, THBS2), immune modulators (CD276, PDCD1LG2, SIGLEC1), and cytotoxic mediators (GZMA, TNFSF11) (Supplementary Figures S5, S7–S9), which mapped to receptors in skin and adipose. Such reciprocal signaling suggests that stressed adipose and muscle compartments feed back into skin to reinforce keratinocyte dysfunction, aberrant angiogenesis, and chronic immune activation.

Across all tissue pairs, lysosomal and metabolic mediators (GRN, PSAP, DHCR24, RARRES2) emerged as recurrent ligands (), highlighting that disrupted metabolism and degradative stress are communicated between compartments alongside classical inflammatory signals. Collectively, these findings reveal that DFUs are sustained by a bidirectional, multi-tissue paracrine network in which skin, adipose, and muscle amplify one another's inflammatory, fibrotic, and metabolic dysfunction. This inter-tissue dialogue likely underlies the persistence of the non-healing wound environment. 1