Integrated transcriptomic, proteomic, and metabolomic analysis unveils key roles of protein and nucleic acid interactions in diabetic ulcer pathogenesis

Background: Diabetes mellitus significantly increases the risk of complications, particularly diabetic foot ulcers (DFUs). However, the underlying mechanism remains unclear. This study aimed to assess the overall therapeutic approach in diabetic ulcers.

Methods: Using integrated high-throughput multi-omics approaches, including transcriptomics, proteomics, and metabolomics, we constructed a compound-reaction-enzyme-gene network to identify the key molecular mechanisms involved in the pathogenesis of DFUs. Major findings were further validated in mouse models of diabetic and control ulcers.

Results: Transcriptomics identified 653 differentially expressed genes (DEGs) between diabetic ulcers and control groups. Pathway analysis indicated that these genes were mostly related to inflammation, including the cytokine–cytokine receptor interaction, TNF signaling pathway, and NF-κB signaling pathway. Proteomics revealed 464 upregulated and 419 downregulated proteins, indicating many differentially expressed proteins (DEPs). The pathways with the highest representation of DEPs included diabetic cardiomyopathy, PPAR signaling pathway, and HIF-1 signaling pathway. Metabolomics identified 1,304 metabolites, predominantly lipids (32.1%) and organic acids (20.2%). Principal component analysis and partial least squares discriminant analysis confirmed the model’s effectiveness in distinguishing sample groups, whereas bioinformatics analysis revealed significant metabolic pathways, particularly amino acid biosynthesis.

Conclusion: Our findings identified critical molecular signatures associated with DFUs and lay the groundwork for developing innovative therapeutic strategies to improve clinical outcomes in patients with this challenging condition.

1 Introduction

Diabetes mellitus (DM) arises from a complex interplay of genetic predisposition, immune dysregulation, infections, lifestyle, and psychosocial factors, leading to insulin resistance or pancreatic dysfunction. Currently, the global prevalence of diabetes among adults aged 20–79 years is estimated at 8.8%, with projections suggesting a potential increase to 693 million cases by 2045 if current trends persist (1). DFUs are significant complications associated with two chronic conditions of diabetes: peripheral neuropathy (PN) and peripheral arterial disease (PAD). Foot ulcers are a prominent manifestation of DFUs (2). DFUs affect up to 15% of individuals with diabetes and are the leading causes of hospitalization and both minor and major amputations in this population (2, 3).

Patients with DFUs often encounter precarious situations because foot ulcers are merely one aspect of a complex clinical picture. This intricacy includes persistent complications such as PN and PAD, compounded by comorbidities that collectively undermine the patient’s overall well-being. Notably, the 5-year mortality rate for individuals who develop new DFUs has been reported to surpass that of several cancers by approximately 25–60% (4, 5). Cardiovascular and renal diseases are leading contributors to mortality in this population. Patients with diabetic feet exhibit a complex interplay of inflammatory markers that adversely affects the cardiovascular system, thereby exacerbating cardiovascular damage and increasing morbidity. Consequently, effective management of DFU requires timely treatment of foot ulcers and a thorough assessment of comorbidities that may influence clinical outcomes (2, 6). However, the underlying mechanisms and potential targets of DFUs require further investigation.

High-throughput multi-omics techniques, including transcriptomics, proteomics, and metabolomics, hold promise for uncovering novel pathological mechanisms and identifying potential therapeutic targets. Transcriptomics enables qualitative and quantitative assessment of mRNA levels across the genome (7). In contrast, proteomics centers on elucidating alterations in protein expression under specific conditions and monitoring dynamic fluctuations within cellular environments (8). Metabolomics entails the comprehensive analysis of endogenous small biomolecules, primarily reflecting the end products of physiological processes mediated by proteins (9). However, a singular omics approach may fail to capture the full spectrum of changes associated with DFUs, given the dynamic and multifaceted nature of the condition. Independent analyses of extensive high-quality data across various omics levels often overlook the intricate interactions between different molecular entities, potentially missing critical biological insights (10). Therefore, a holistic strategy integrating data from multi-omics platforms is imperative for a comprehensive understanding of key pathological processes. Such integration can reveal new perspectives on complex biological systems and clarify networks of interactions at the molecular level.

In this study, we developed a comprehensive strategy to investigate diabetic ulcers, culminating in establishing a compound-reaction-enzyme-gene network. This network was constructed by integrating transcriptomic, proteomic, and metabolomic data to identify key targets and mechanisms for diabetic ulcer treatment. This approach may facilitate a deeper understanding of the molecular mechanisms underlying diabetic ulcers. Our findings present a detailed molecular map of these ulcers, advancing knowledge of their pathogenesis, and paving the way for exploring novel therapeutic interventions for affected patients.

2 Experimental methodologies

2.1 Clinical research

2.1.1 Data sources and study population

The National Health and Nutrition Examination Survey (NHANES) comprises a series of cross-sectional, population-based studies aimed at evaluating the health and nutritional status of adults and children in the United States. The National Center for Health Statistics Research Ethics Review Board approved the NHANES study protocol, and all participants provided written informed consent. In accordance with the guidelines issued by the National Center for Health Statistics (NCHS), secondary use of NHANES data does not require additional IRB approval. This cross-sectional analysis utilized NHANES data (www.cdc.gov/nchs/nhanes/) collected between 1999 and 2004. A flowchart of the study design is shown in Figure 1A. Exclusion criteria included: (1) age under 40 years; (2) missing data on triglycerides, fasting glucose, body mass index (BMI), and waist circumference (WC); and (3) no DM diagnosis. Ultimately, 31,126 participants were enrolled, with 1,275 meeting the inclusion criteria. Among these, 100 participants (7.8%) were diagnosed with non-healing lower extremity ulcers (NHLU).

Figure 1

Figure 1. (A) Flowchart illustrating the methodology of this study. (B) Receiver Operating Characteristic (ROC) curves depicting the predictive utility of the TyG index, TyG-WC, and TyG-BMI for all-cause mortality in NHLU. (C) Smooth curve fitting analysis, with the solid red line representing the fitted curve between variables, and the blue bands indicating the 95% confidence interval surrounding this fit. (D) Development of a mouse model to study diabetic ulcers (this figure is created by biorender) and the representative wound healing photographs. (E) Representative images of ulcerated wounds across three groups of mice. (F) Assessment of dermal thickness among the three groups of mice using hematoxylin and eosin (HE) staining, scale bar = 50 μm.

2.1.2 Predictor and outcome variables

The primary outcome variable for this study was the presence of NHLU, defined by an affirmative response to the following question: “Have you had an ulcer or sore on your leg or foot that took more than four weeks to heal?” In recent years, studies have shown that triglyceride-glucose (TyG) index could be used as a reliable marker of insulin resistance (11, 12). Meanwhile, some reports have shown that TyG is associated with stroke and early diabetic nephropathy and is a reliable predictor of cardiovascular and all-cause mortality in prediabetes (13–15). The TyG index and its variants, TyG-WC and TyG-BMI, were calculated as follows: (1) TyG = ln [triglycerides (mg/dL) × glucose (mg/dL)/2]; (2) TyG-BMI = TyG × BMI; and (3) TyG-WC = TyG × WC. Additional covariates examined included age, sex, glycosylated hemoglobin, total cholesterol (TC), HDL, LDL, C-reactive protein (CRP), BMI, WC, systolic blood pressure (SBP), diastolic blood pressure, smoking history, and alcohol consumption.

2.2 Animals and regents

Twenty healthy male C57BL/6J mice (age: 8–10 weeks; body weight: 30–35 g) were purchased from the Beijing Viton Lever Company and housed in the Medical Research Center animal facility at the First Affiliated Hospital of Shandong First Medical University. The mice were divided into three cages, with five mice per cage. Two cages received a high-sugar, high-fat diet comprising 66.5% rat and mouse maintenance chow, 10% lard, 20% sucrose, 25% cholesterol, and 1% sodium cholate, while one cage received standard chow. Experiments commenced after a 4-week acclimatization period. All procedures involving mice followed the ARRIVE guidelines established by the Animal Protection Society of the First Affiliated Hospital of Shandong First Medical University (Approval no. SYDWLS【2021】002).

2.3 Establishment of the diabetic ulcer mouse model

Ten C57BL/6J mice were fed a high-fat, high-sugar diet for 4 weeks, followed by a 12-h fasting period. They were then intraperitoneally injected with 2% streptozotocin (STZ) solution at 35 mg/kg for five consecutive days. Three days post-injection, blood glucose levels were measured in the tail vein to confirm the establishment of the diabetic mouse model. Subsequently, an oral glucose tolerance test was conducted to verify successful diabetes induction. Following the establishment of diabetes, hyperglycemia was maintained for an additional 2 weeks.

A skin ulcer model was created using ten diabetic and five normal mice. Anesthesia was induced through an intraperitoneal injection of 1% sodium pentobarbital, after which the mice were depilated and disinfected. Two points were marked on each mouse’s back, approximately 7 mm from the midline and 4 cm from the base of the neck, serving as the centers for wound creation. A sterile, disposable 5 mm biopsy punch was used to delineate the circular wound, and iris scissors (elbow type) were used to excise the tissue, creating an ulcer. The duration of wound induction was designated as 0 h. Afterward, the mice were housed individually, and photographs were taken and recorded on days 0, 3, and 7 (Figures 1E, F). Ulcer size in each group was observed, and the wound healing rate was calculated using ImagePro Plus 6.0, with the formula: wound healing rate (%) = [(initial wound area - wound area on the observation day)/initial wound area] × 100%. The diabetic ulcer mouse model and the representative wound healing photographs are illustrated in Figure 1D.

2.4 Multi-omics experimental phase investigation

2.4.1 Transcriptomic analysis through RNA sequencing

To investigate transcriptional changes, high-throughput RNA sequencing was performed using the Illumina NovaSeq 6000 platform with 150 bp paired-end reads. Total RNA was extracted and subjected to quality assessment for purity and integrity using the NanoDrop 2000 (Thermo Fisher Scientific, USA) and LabChip GX Touch (PerkinElmer, USA). Ribosomal RNA (rRNA) was depleted using the Ribo-Zero Gold Kit (Illumina, USA) to enrich for non-rRNA transcripts. Strand-specific cDNA libraries were then constructed using the NEBNext Ultra Directional RNA Library Prep Kit (NEB, USA) following the manufacturer’s protocol. The library preparation workflow included RNA fragmentation, end repair, adapter ligation, size selection, and PCR enrichment.

The quality and concentration of the libraries were further evaluated using Qubit 3.0 (Invitrogen, USA), Agilent 2100 Bioanalyzer (Agilent Technologies, USA), and Bio-RAD CFX96 system (Bio-Rad, USA). Libraries were pooled based on effective concentration and sequencing depth requirements before being subjected to sequencing, which employed the sequencing-by-synthesis principle to generate high-quality reads. Raw sequencing data were initially processed using FastQC and Trim Galore to remove adapters, ambiguous reads (N bases), and low-quality sequences (Qphred ≤ 20 over more than 50% of the read length). Clean reads were aligned to the human reference genome (GRCh38) using HISAT2 (v2.2.1). Transcript assembly and quantification were performed with StringTie, and gene expression levels were calculated using both FPKM (Fragments Per Kilobase of transcript per Million mapped reads) and TPM (Transcripts Per Million).

Differential expression analysis was performed using DESeq2 for samples with biological replicates, applying negative binomial distribution models and adjusting P-values through the Benjamini–Hochberg method. Genes with adjusted P-values ≤ 0.05 were classified as differentially expressed. DESeq was used for samples without replicates. Finally, differential gene enrichment analysis was conducted using ClusterProfiler to identify significant Gene Ontology (GO) terms and KEGG pathways, applying a threshold of P-values < 0.05. The transcriptomic data analysis process is illustrated in Supplementary Figure S1A.

2.4.2 Proteomic analysis

Protein extraction and digestion began with sample lysis using SDT buffer (4% (w/v) SDS, 100 mM Tris/HCl, pH 7.6, 0.1 M DTT), followed by protein quantification using a BCA Protein Assay Kit (Bio-Rad, USA). Proteins were digested according to the filter-aided sample preparation method described by Mann (16). Specifically, 200 μg of protein was treated with UA buffer (8 M Urea, 150 mM Tris-HCl, pH 8.0) to eliminate low-molecular-weight components using Microcon units (10 kD). This was succeeded by treatment with 100 μL of iodoacetamide (100 mM IAA in UA buffer) to block cysteine residues. An overnight trypsin digestion was performed, and the resulting peptides were desalted using Empore™ SPE Cartridges C18 (standard density, bed I.D. 7 mm, volume 3 mL, Sigma), concentrated, and reconstituted in 40 µL of 0.1% (v/v) formic acid. For SDS-PAGE analysis, 20 μg of each sample was combined with 5X loading buffer and boiled for 5 min.

Proteins were resolved on a 12.5% SDS-PAGE gel, with a constant current of 14 mA for 90 min, and visualized using Coomassie Blue R-250 staining. Peptides were then labeled with tandem mass tag (TMT) reagents (Thermo Scientific) and subjected to high-pH reversed-phase fractionation using a kit from Thermo Scientific, resulting in 10 separate peptide fractions. If there are too many samples, a bridge will be used for correction to reduce the batch effect. These fractions were desalted on C18 cartridges (Empore SPE Cartridges C18, standard density, bed I.D. 7 mm, volume 3 mL, Sigma) and concentrated by vacuum centrifugation. Liquid chromatography-tandem mass spectrometry (LC-MS/MS) analysis was performed using a Q Exactive mass spectrometer (Thermo Scientific), with a linear gradient of acetonitrile for peptide separation and data acquired in a data-dependent manner for identification and quantification. Raw MS data were processed using the MASCOT engine (version 2.2; Matrix Science, London, UK) in Proteome Discoverer 1.4 for protein identification and quantification.

After completing these steps, we conducted bioinformatics analysis to elucidate protein characteristics and interactions. Hierarchical clustering was performed using Cluster 3.0 (http://bonsai.hgc.jp/~mdehoon/software/cluster/software.htm) and Java Treeview software (http://jtreeview.sourceforge.net), with the Euclidean distance algorithm and average linkage clustering. The results were presented as heatmaps and dendrograms. Protein subcellular localization was predicted using the CELLO multiclass SVM classification system (http://cello.life.nctu.edu.tw/). Domain annotation was performed using InterProScan to identify the protein domain signatures from the Pfam database. GO annotations were derived using NCBI BLAST+ (ncbi-blast-2.2.28+-win32.exe) and InterProScan, with mapping and visualization achieved using Blast2GO and R scripts. KEGG annotations were obtained by blasting proteins against the KEGG database (http://geneontology.org/) to identify orthologs and associated pathways. Enrichment analyses were conducted using Fisher’s exact test and the Benjamini–Hochberg correction for multiple testing, highlighting functional categories and pathways with P-values < 0.05. Finally, protein–protein interactions were explored using the IntAct database (17) (http://www.ebi.ac.uk/intact/) and STRING (http://string-db.org/), with networks visualized in Cytoscape (http://www.cytoscape.org/, version 3.2.1) to assess the significance of each protein within the interaction network based on degree. The proteomic analysis is illustrated in Supplementary Figure S1B.

2.4.3 Metabolomic analysis

Chemicals: Ammonium acetate was procured from Sigma Aldrich, acetonitrile from Merck, and ammonium hydroxide and methanol from Fisher Scientific. After dissection, mouse tissues were rapidly frozen in liquid nitrogen. Approximately 80 mg of tissue was diced on dry ice, placed in a 2 mL Eppendorf tube, and homogenized with a mixture of 200 μL H2O and five ceramic beads using a homogenizer. For metabolite extraction, 800 μL of methanol/acetonitrile (1:1, v/v) solution was added to the homogenate. The mixture was centrifuged at 14,000 g for 20 min at 4°C, and the supernatant was dried using a vacuum centrifuge. For LC-MS analysis, the dried samples were re-dissolved in 100 μL of a 1:1 (v/v) acetonitrile/water solution and centrifuged again at 14,000 g for 15 min at 4°C, and the supernatant was injected for analysis.

To ensure the stability and reproducibility of the instrument analysis, quality control samples were generated by pooling 10 μL from each sample. In this experiment, the number of peaks with RSD ≤ 30% in the QC samples accounted for more than 80% of the total number of peaks in the QC samples, indicating that the instrument analysis system has good stability and the data can be used for subsequent analysis. These quality control samples were incorporated into the analysis at regular intervals and evaluated after every five samples. Chromatography-mass spectrometry analysis was performed using a Vanquish ultra-high-performance liquid chromatography (UHPLC) system with a HILIC column (Vanquish UHPLC, Thermo Scientific), with samples processed randomly. A Q Exactive series mass spectrometer was used to collect both primary and secondary spectra. The raw data were then converted to mzXML format using ProteoWizard MSConvert. Peak alignment, retention time correction, and peak area extraction were conducted using XCMS. The data extracted by XCMS were subjected to metabolite structure identification and pre-processing, followed by quality evaluation and subsequent analyses.

After sum normalization, data were analyzed using the R package (ropls) for multivariate analysis, including principal component analysis (PCA) and orthogonal partial least-squares discriminant analysis (OPLS-DA). Model robustness was evaluated through 7-fold cross-validation and response permutation testing. Variable importance in projection (VIP) values were calculated to determine each variable’s contribution to classification, with significant metabolites identified using VIP > 1 and P < 0.05. Student’s t-test was used to determine significance between independent sample groups, and Pearson’s correlation analysis was performed to explore the relationships between variables. The metabolomic analysis is illustrated in Supplementary Figure S1C.

2.5 Statistical analysis

In this study, continuous variables were presented as means with standard deviations and categorical variables as counts with their corresponding proportions. Differences between the two groups were analyzed using weighted linear regression for quantitative variables and weighted chi-squared tests for qualitative variables. We utilized multiple logistic regression models to estimate the association of TyG, TyG-BMI, and TyG-WC with NHLU, with the first group as a reference. We developed three multivariate models to further elucidate the clinical implications of the logistic regression findings. Model 1 incorporated TyG, TyG-BMI, and TyG-WC as predictors of NHLU. Model 2 was adjusted for age and sex, while Model 3 included additional covariates, including smoking status, alcohol consumption, fasting blood glucose (FBG), hemoglobin A1c (HbA1c), triglycerides (TG), TC, CRP, and SBP. Analyses were conducted using EmpowerStats (http://www.empowerstats.com/cn/) and R. All probability values were two-sided, with interaction P-values < 0.05 considered statistically significant.

3 Results

3.1 Clinical characteristics

This study enrolled a total of 1,275 patients with type 2 diabetes mellitus (T2DM), of whom 100 presented with NHLU. Table 1 shows the clinical characteristics of the participants categorized by NHLU status. Significant differences were identified in WC, FBG, TyG index, TyG-WC, and TyG-BMI (all with P < 0.05). Table 2 presents the prevalence of NHLU among patients with T2DM, stratified by quartiles of TyG index, TyG-WC, and TyG-BMI. The Q4 group demonstrated the highest prevalence of NHLU across all three indicators, with significant differences observed in each case (P < 0.05). Specifically, within the TyG-WC quartiles, NHLU prevalence was the highest in group Q4, followed by Q3, Q2, and Q1, all with significant differences (P < 0.05). For the TyG quartiles, NHLU prevalence was significantly higher in Q4 and Q3 than in Q2 and Q1 (P < 0.05). Additionally, among those grouped by TyG-BMI quartiles, Q4 exhibited a greater NHLU prevalence than Q1, with significant differences observed between the groups (P < 0.05).

Table 1

Table 1. Clinical characteristics of participants with and without non-healing lower extremity ulcers (NHLU).

Table 2

Table 2. Clinical characteristics of participants with and without NHLU.

3.2 Multivariate analysis of determinants of NHLU in study participants

A multivariate logistic regression model was constructed to investigate the associations between the TyG index, TyG-WC, TyG-BMI, and NHLU in patients with DM, as outlined in Table 3. In Model 1, no adjustments were made for confounding variables. Elevated levels of the TyG index, TyG-WC, and TyG-BMI (group Q4) exhibited significant positive correlations with NHLU when compared to group Q1, yielding odds ratios of 0.59 (95% CI: 0.02–1.17; P = 0.0440), 0.91 (95% CI: 0.33–1.49; P = 0.0022), and 0.94 (95% CI: 0.33–1.55; P = 0.0023), respectively. In Model 2, adjusted for age and sex, the associations remained significant. The highest TyG index, TyG-WC, and TyG-BMI levels in group Q4 remained positively associated with NHLU compared to the reference group, with odds ratios of 0.60 (95% CI: 0.02–1.19; P = 0.0441), 1.03 (95% CI: 0.40–1.65; P = 0.0013), and 1.19 (95% CI: 0.53–1.86; P = 0.0004), respectively. In Model 3, adjusted for age, sex, smoking status, alcohol consumption, FBG, HbA1c, TG, TC, CRP, and SBP, the highest TyG-WC and TyG-BMI levels remained positively associated with NHLU, with odds ratios of 0.93 (95% CI: 0.10–1.77; P = 0.0289) and 1.25 (95% CI: 0.34–2.15), respectively. However, no significant relationship was observed between the TyG index and NHLU.

Table 3

Table 3. Logistic regression analysis of different indexes on the risk of NHLU in T2DM patients.

3.3 Predictive value of TyG-related indices for NHLU in patients with T2DM

To evaluate the predictive capacity of the TyG index, TyG-WC, and TyG-BMI for NHLU in patients with DM, we conducted a receiver operating characteristic (ROC) curve analysis (Figure 1B; Table 4). The identified cutoff values were as follows: TyG index at 9.43 (sensitivity, 55%; specificity, 61%; area under the curve [AUC], 0.58), TyG-WC at 1150.80 (sensitivity, 34%; specificity, 85%; AUC, 0.61), and TyG-BMI at 307.22 (sensitivity, 46%; specificity, 70%; AUC, 0.59). The predictive model showed an overall sensitivity of 44.6% and an AUC of 0.607.

Table 4

Table 4. Ability of TyG, TyG-BMI, TyG-WC to predict NHLU in DM patients.

3.4 Linear relationship between TyG-WC and TyG-BMI in DM and NHLU

We conducted a smoothed curve-fitting analysis to assess the relationship between TyG indices in DM and NHLU (Figure 1C). After adjusting for all covariates in Model 3, we identified a significant linear relationship between TyG-WC and TyG-BMI in the context of DM and NHLU (P < 0.05).

3.5 Transcriptomic analysis of the diabetic ulcers mice model

We assessed the expression levels of the samples using the formula:

$$FPKM=\frac{{10}^3\ast F}{NL/{10}^6}.$$

The resulting expression density map illustrated gene expression concentrations within the peak area (Figure 2). A box plot was generated to illustrate the expression levels across samples and revealed a consistent trend (Figure 2). Additionally, a saturation analysis confirmed that the data volume in this study was adequate for further analysis (Figure 2).

Figure 2

Figure 2. (A) Expression distribution diagram. The horizontal axis is log2(FPKM), which represents the logarithmic value of the gene expression level of each sample. The higher the value, the higher the gene expression level. The vertical axis is the ratio of the number of genes at the corresponding expression level to the total number of genes. (B) Box plot of expression level. The horizontal axis is the sample name, and the vertical axis is the logarithmic value of the gene expression of each sample. From top to bottom, they represent the maximum value, upper quartile, median, lower quartile and minimum value. (C) Saturation analysis chart. Different colors represent genes classified based on expression levels. The horizontal axis is the ratio of sampled reads to the total number of reads during random sampling, and the vertical axis represents the percentage of genes whose expression levels could be accurately estimated through sampled reads. (D) Statistical graph of differentially expressed genes. The horizontal axis represents different sets of differentially expressed genes, green represents all differentially expressed genes, blue represents upregulated genes, yellow represents downregulated genes, and the vertical axis represents the number of differentially expressed genes. (E) Volcano plot of DEGs. The horizontal axis represents the expression fold change of genes in different samples, and the vertical axis represents the statistical significance of the expression change. Blue represents significantly downregulated genes, yellow represents significantly upregulated genes, and gray represents genes with no significant change in expression. (F) Differential gene clustering diagram. Behavioral genes are listed as samples. Red indicates high gene expression, and blue indicates low gene expression. The horizontal axis represents the clustering of samples, and the vertical axis represents the clustering of genes. (G) GO statistical bar chart of DEGs. The horizontal axis is GO Term, the left vertical axis is the percentage of the number of genes, and the right vertical axis is the number of genes. (H) GO enrichment bar chart. The vertical axis represents the GO entry, the horizontal axis represents the number of genes enriched in the entry, and the color represents p.adjust. The redder the color, the more significant it is. (I) GO enrichment bubble chart. The bubble size indicates the number of genes enriched in this entry. The larger the bubble, the more genes there are. (J) Distribution of p.adjust values of enriched GO terms (including BP, CC and MF). (K) KEGG enrichment bar chart. (L) KEGG enrichment bubble chart. (M) Distribution of p.adjust values of enriched pathways. Different colors represent different degrees of enrichment; the redder the color, the more significant the enrichment. (N) Statistics of variable splicing results. The horizontal axis represents the sample name, and the vertical axis represents the percentage of alternative splicing event types. (O) Statistics of differential alternative splicing results. (P) Distribution diagram of mutation statistics. The horizontal axis represents the sample name, and the vertical axis represents the number of mutations. (Q) SNP mutation frequency distribution and InDel length distribution. A>C represents the number of SNP sites that mutated from A to C. (R) Distribution diagram of differential transcription factor protein families. The vertical axis represents different transcription factor families, and the horizontal axis represents the number of differentially expressed genes annotated to the transcription factor. (S) Cluster diagram of differential transcription factors. Red indicates high gene expression, and blue indicates low gene expression. The horizontal axis represents the clustering of samples, and the vertical axis represents the clustering of genes.

Differential gene screening identified 553 upregulated and 100 downregulated genes in the TC group compared with the TST group, totaling 653 differentially expressed genes (DEGs) (Figure 2). DEGs were defined as genes showing statistically significant differences in expression levels between two comparison groups, identified using the DESeq2 package with a threshold of |log₂ fold change| > 1 and adjusted p-value < 0.05. A volcano plot was generated to visually represent the distribution of these DEGs across groups (Figure 2). To further elucidate gene expression differences and identify novel functional genes, hierarchical clustering analysis was performed on all screened DEGs. Genes exhibiting similar expression patterns across samples were clustered and visualized using heatmaps (Figure 2; Data S12). Following DEG identification, we conducted GO functional enrichment and KEGG pathway enrichment analyses to elucidate the fundamental molecular mechanisms underlying the biological processes involved. The GO statistical results, summarized in a bar graph (Figure 2), indicate that molecular functions were predominantly enriched in binding, cellular components were primarily associated with cellular, anatomical entities, and biological processes largely involved cellular processes and biological regulation. The 30 most significant GO terms were visualized in bar and bubble charts (Figures 2H; Data S13). Additionally, a distribution diagram illustrates the adjusted enrichment significance values across samples (Figure 2).

Subsequently, we identified the 30 most significant pathways (or all if fewer than 30) through KEGG enrichment analysis, with results indicating significant enrichment in pathways such as cytokine–cytokine receptor interaction, NF-κB signaling, and TNF signaling (Figures 2K; Data S14). A distribution map illustrates the adjusted P-values of the enriched pathways across all comparison groups (Figure 2).

Alternative splicing is a critical regulatory mechanism that contributes to gene expression diversity and is essential for growth and development. Using ASProfile software, we analyzed and quantified the variable splicing events per sample based on known gene models (Figure 2). Additionally, rMATS was used to classify and count alternative splicing events for each comparison group (Figure 2). We also investigated single nucleotide polymorphisms (SNPs), which represent allele mutations across the genome. Using SAMtools, we compared the alignment files, post-sorting and PCR duplication removal against the reference sequence to detect variations, revealing approximately 3,000 SNP mutations (Figure 2). We subsequently recorded the frequency and length of each mutation based on the identified SNP sites (Figure 2; Supplementary Figure S2).

Transcription initiation in eukaryotes is complex and often involves multiple protein factors. We predicted transcription factors and compared our findings with the Animal Transcription Factor Database (AnimalTFDB) to categorize these factors by family. The distribution of the predicted transcription factors mapped to their respective protein families is illustrated in Figure 2, with notable prevalence in the zf-C2H2 and TF_bZIP families. To reflect the expression patterns of the differentially expressed transcription factors under varying experimental conditions, we conducted a hierarchical cluster analysis using R. This analysis revealed that transcription factors within the same cluster exhibited similar expression trends under the same experimental conditions (Figure 2).

3.6 Proteomic analysis of diabetic ulcers mice model

Through proteomic analysis, we generated a total of 811,606 secondary spectra, which included 80,869 spectra matched to the database, 31,709 peptides (of which 28,252 were unique), and 4,985 identified proteins, of which 4,983 were quantifiable (Supplementary Figure S3; Data S1 and S2). To investigate differential protein expression across groups, we applied a 1.2-fold change (FC) threshold at P < 0.05, identifying 883 differentially expressed proteins (DEPs) between the PC and PST groups. DEPs refer to proteins whose abundance significantly varies between experimental conditions, as determined from proteomic datasets using Student’s t-test or ANOVA, with FDR correction applied where appropriate. These DEPs included 464 upregulated and 419 downregulated proteins (Figure 3; Data S3).

Figure 3

Figure 3. (A) Bar graph of protein quantification difference results. (B) Volcano plot of the PSTvsPC group. The horizontal axis is the difference fold (logarithmic transformation with base 2), and the vertical axis is the significant P-value of the difference (logarithmic transformation with base 10). The red points in the figure are upregulated proteins with significant differential expression, the blue points are downregulated proteins with significant differential expression, and the gray points are proteins with no difference. (C) PSTvsPC group differentially expressed protein clustering analysis results. Hierarchical clustering results are presented in a tree-type heatmap, where each column represents a group of samples (the horizontal axis is sample information), and each row represents a protein (i.e., the vertical axis is the protein with significant differential expression). (D) Pie chart of subcellular localization of differentially expressed proteins in the PSTvs PC group. (E) PSTvsPC group differentially expressed protein domain analysis diagram. (F) PSTvsPC group domain enrichment analysis. The horizontal axis in the figure is the enrichment factor (Rich Fator ≤ 1), which indicates the ratio of differentially expressed proteins annotated to the GO category to the number of all identified proteins annotated to the category. The vertical axis indicates the statistical results of differential proteins under each domain classification; the color of the bubble indicates the significance of the enriched domain classification. (G) GO annotation statistics of differentially expressed proteins in the PSTvsPC group (level 2). (H) GO annotation statistics of the top 20 differentially expressed proteins in the PSTvsPC group. (I) A diagram of biological metabolic processes showing the relationship between different metabolic processes and the relevant values. (J) KEGG pathway enrichment bubble diagram of PSTvsPC group. (K) PSTvsPC group differentially expressed protein interaction network diagram. (L) Interaction network diagram for functional classification.

To identify significant protein differences, we generated a volcano plot, displaying FC and P-values from t-tests. Proteins exhibiting substantial downregulation are indicated in blue (FC < 0.83, P < 0.05), whereas significantly upregulated proteins are indicated in red (FC > 1.2, P < 0.05). Proteins without significant differences are represented in gray (Figure 3). Subsequent cluster analysis revealed distinct differences in protein expression between the PC and PST groups (Figure 3).

Subcellular localization of all DEPs was analyzed using the CELLO subcellular structure prediction software (18). The results indicated that most DEPs were localized in the cytoplasm (356 proteins, 46%), followed by the nucleus (171 proteins, 22%) (Figure 3; Data S4). To predict DEP domains, we utilized InterProScan to identify the top 20 protein domains, with the most prominent being intermediate filament proteins, protein kinase domains, and trypsins (Figure 3). To elucidate the domain enrichment characteristics of the DEPs, we conducted domain enrichment analysis using Fisher’s exact test, revealing significant enrichment in sushi repeat (SCR repeat) and trypsin domains (Figure 3).

To comprehensively understand the function, localization, and biological pathways of the proteins, we performed GO functional annotation using Blast2Go (https://www.blast2go.com/) and assessed secondary functional annotation levels (Figure 3). Fisher’s exact test was applied to identify enriched functional categories among DEPs, revealing significant alterations in key biological processes, including generating precursor metabolites and energy, defense responses, and purine nucleotide metabolism. Additionally, significant changes were observed in molecular functions, including peptidase and endopeptidase inhibitor activities, actin binding, and localization within components like the oxidoreductase complex and extracellular space (Figure 3; Data S5).

To illustrate the hierarchical relationships among GO terms associated with DEPs, we employed a topGO-directed acyclic graph. This top-down visualization presents the functional scope, with branches reflecting inclusion relationships and lower branches corresponding to more specific functional categories (Supplementary Figure S4). The proteins were systematically annotated using the KEGG pathway database (Figure 3), tallying the DEPs associated with each pathway, and revealing diabetic cardiomyopathy (DCM), the PPAR signaling pathway, and the HIF-1 signaling pathway as the most prominent pathways (Figure 3; Data S6). Visual representations of the PPAR signaling pathway and HIF-1 signaling pathway are shown in Supplementary Figure S5.

Given that highly aggregated proteins may share similar functions and that highly connected proteins can act as pivotal nodes influencing metabolic or signaling pathways, we further investigated protein–protein interactions. The resulting interaction network diagram of DEPs in the PST vs. PC group, along with specific clusters, is presented in Figures 3, L.

3.7 Metabolomic analysis of diabetic ulcers mice model

For metabolite identification, we used an in-house database (19, 20) to match retention times, molecular masses (error margin < 10 ppm), secondary fragmentation spectra, and collision energy data from biological samples. A rigorous manual validation confirmed that all identifications achieved a minimum Level 2 classification. After integrating positive and negative ion modes, 1,304 metabolites were identified, with the breakdown of metabolites detected in each mode presented in Table 5 and detailed qualitative and quantitative results in Data S7. Metabolites were categorized by their chemical taxonomy, revealing lipids and lipid-like molecules as the largest proportion (32.055%), followed by organic acids and their derivatives (20.245%) (Figure 4).

Table 5

Table 5. Statistics of metabolites identified in positive and negative ion modes.

Figure 4

Figure 4. (A) The proportion of identified metabolites in each chemical category. The blocks of different colors in the figure represent different chemical classification entries, and the percentage represents the percentage of metabolites in the chemical classification entry to the total number of identified metabolites. Metabolites without chemical classification are defined as undefined. (B) Volcano plots of positive and negative ion modes. Upregulated metabolites are represented by rose red, downregulated metabolites are represented by blue, and non-significantly different metabolites are represented by black. (C) Negative ion mode volcano plot (colors are related to the chemical classification of differential metabolites). (D) PCA score plots for positive and negative ion modes. t [1] represents principal component 1, t [2] represents principal component 2, and the ellipse represents the 95% confidence interval. Points of the same color represent the biological replicates within the group, and point distribution reflects the differences between and within the groups. (E) PLS-DA score plots for positive and negative ion modes. (F) PLS-DA permutation test in positive and negative ion modes. (G) OPLS-DA score plots for positive and negative ion modes. (H) OPLS-DA permutation test in positive and negative ion modes. (I) Analysis of fold difference in expressing significant differential metabolites in positive and negative ion modes. (J) Hierarchical clustering heat map of significantly different metabolites in positive and negative ion modes. Each row in the figure represents a differential metabolite (i.e., the vertical axis is the metabolite with significant differential expression), and each column represents a group of samples (i.e., the horizontal axis is the sample information). (K) Positive and negative ion mode correlation heatmap. Red indicates a positive correlation, blue indicates a negative correlation, and white indicates a non-significant correlation. (L) Positive and negative ion mode chord diagram. The starting point of the inner circle link in the figure represents each significantly different metabolite, and the arc on the outer circle represents the category to which the significantly different metabolites belong. The colored lines represent the correlation within each metabolite, and the lines are the same color as the subclasses. The dark gray lines represent the correlation between metabolites of different categories. (M) Positive and negative ion mode network diagram. The dots in the figure represent significantly different metabolites. The size of the dots is related to the degree of connectivity. The larger the degree, the larger the dots. The color of the line represents the correlation. Red indicates positive correlation, and blue indicates negative correlation. The thickness of the line represents the absolute value of the correlation coefficient. The thicker the line, the greater the correlation. (N) KEGG enrichment pathway diagram (bar graph). The vertical axis in the bar graph represents each KEGG metabolic pathway, and the horizontal axis represents the number of differentially expressed metabolites contained in each KEGG metabolic pathway. (O) The differential abundance score map of all enriched metabolic pathways and the score map after classification according to Pathway_Hierarchy.

Metabolomic data, characterized by high dimensionality and strong intervariable correlations, were analyzed using multivariate statistical methods to integrate univariate analyses with multidimensional assessments, identifying group differences and potential biomarkers. Differential analysis included all metabolites detected in both ion modes and those that remained unidentified, visualized as volcano plots with significant changes indicated by FC > 1.5 or < 0.67 and P < 0.05 (Figure 4). We used distinct colors to visually represent the classification of these differential metabolites, as depicted in Figure 4.

PCA, an unsupervised method, was employed to linearly combine all identified metabolites into a new set of variables. This method aims to reflect as much information as possible from the original variables while achieving dimensionality reduction. PCA results provided insights into overall distribution trends and inter-group differences (Figure 4), with model parameters validated through a seven-fold cross-validation (Data S8a; b), indicating an R² value approaching 1, confirming the reliability of our findings. Subsequently, we established a discriminant model using Partial Least Squares Discriminant Analysis (PLS-DA) to identify the differential lipid substances associated with the groups (Figure 4). The PLS-DA model effectively separated the two sample groups, with robust stability confirmed by evaluation parameters (R²Y, Q²) (Data S8c, d). To mitigate the risk of overfitting in the supervised model, we performed a permutation test (Figure 4). The decreasing R² and Q² values of the random models confirmed the absence of overfitting and the model’s robustness.

To enhance the rigor of our analysis, we applied OPLS-DA, which further distinguished the two sample groups. Circular cross-validation confirmed model stability and reliability (Figure 4; Data S8e, f), and a permutation test validated the absence of overfitting (Figure 4). Following common metabolomic practices, we employed stringent criteria (VIP > 1 and P < 0.05) to screen for significantly differential metabolites (Figure 4; Data S9). Subsequently, we performed bioinformatics analyses, including cluster, correlation, and pathway analyses, on these metabolites.

To comprehensively illustrate sample relationships and metabolite expression patterns, we clustered and analyzed the expression levels of all samples with the differential metabolites (Figure 4). Correlation analysis assessed the metabolic proximity of significantly different metabolites (VIP > 1, P < 0.05), visualized in chord and network diagrams derived from the correlation matrix (Figures 4–M). Various metabolites coordinate their functions within biological systems, and KEGG pathway analysis further elucidated their biological roles. Before KEGG pathway annotation, we merged differential metabolites identified from both ion modes (Data S10). To facilitate observation of each differential metabolite expression within the KEGG metabolic pathways, we selected pathways featuring more than five differential metabolites and represented them in a heatmap (Supplementary Figure S6, Data S11). KEGG enrichment analysis indicated significant enrichment in the amino acid biosynthesis pathway (Figure 4). Finally, a Differential Abundance Score pathway-based analysis captured overall metabolic changes, as presented in Figure 4.

4 Discussion

With the rising prevalence of obesity and an aging population, the incidence of diabetes has surged, posing a significant public health challenge (21). Individuals with diabetes frequently encounter a myriad of complex complications, including chronic pain stemming from diabetic ulcers and foot conditions, substantially elevating the risks of limb amputation and mortality (22). Unlike typical wounds, diabetic ulcers often progress to chronic, non-healing lesions owing to a multifaceted inflammatory microenvironment (23), characterized by factors like excessive ROS accumulation (24) and hypoxia (25). These elements collectively disrupt skin regeneration. Moreover, the impaired vasculature associated with diabetic wounds restricts blood flow, reducing oxygen supply and exacerbating inflammation around the ulcer (26). Despite available treatments for diabetic ulcers, including debridement, growth factor therapy, and topical antimicrobials, effective wound healing remains challenging (27). Further investigation into the molecular mechanisms underlying these processes may uncover potential therapeutic targets for enhancing wound healing in patients with diabetes.

Integrating NHANES database, transcriptomics, proteomics, and metabolomics has provided profound insights into the complex mechanisms underlying diabetic ulcer pathogenesis. In this study, we used normal and diabetic ulcer mice models (n = 4 for transcriptomics, n = 5 for proteomics, and n = 7 for metabolomics) to conduct a comprehensive analysis across these multi-omics platforms, aiming to uncover new dimensions in the healing landscape of DFUs. It is noteworthy that several estimates presented in Table 3 exhibit wide confidence intervals (e.g., TyG-WC Q4: 95% CI, 0.10–1.77), suggesting potential instability in the effect sizes. This variability may, in part, reflect heterogeneity in baseline metabolic profiles (e.g., BMI, glycemic status) and lifestyle factors (e.g., dietary patterns, physical activity) among the study population. Despite the implementation of stratified analyses and multivariable adjustments to account for known covariates, the influence of unmeasured confounders—such as genetic predispositions or varying environmental exposures—cannot be fully excluded. Furthermore, while standard metabolic variables (e.g., age, sex, smoking) were incorporated into the model, certain relevant factors—including indices of insulin resistance and systemic inflammation—were not comprehensively adjusted due to missing data. These limitations may have contributed to residual confounding, thereby broadening the confidence intervals. Future investigations should prioritize subgroup analyses to mitigate population-level heterogeneity and incorporate a more robust panel of metabolic and inflammatory biomarkers to enhance model stability and interpretability.

Pathway analysis of the transcriptomic data revealed that significantly altered signaling pathways were predominantly associated with endothelial cell-mediated inflammatory responses, particularly the NF-κB signaling pathway, cytokine–cytokine receptor interactions, and TNF signaling pathway. Various herbal extracts can enhance wound healing by modulating the NF-κB pathway. For instance, balsam pears facilitate diabetic ulcer healing in mouse models by influencing the RAGE/NF-κB and VEGF/VCAM-1/eNOS signaling pathways (28). Additionally, scutellarein inhibits NF-κB-mediated luciferase expression and nuclear translocation, while reducing the phosphorylation of upstream signaling enzymes. The significant suppression of Src kinase activity further underscores its potential as an anti-inflammatory agent (29). The TNF pathway is crucial in initiating and progressing inflammatory and autoimmune disorders, cancer, and cardiovascular conditions (27, 30).

Regarding cytokine–cytokine receptor interactions, miRNA-497 treatment, both in vivo and in vitro, can reduce levels of pro-inflammatory cytokines such as IL-1β, IL-6, and TNF-α, thereby promoting DFU healing (31). In our investigation, we identified the pro-inflammatory gene IGFN1 as significantly differentially expressed in both the diabetic and normal ulcer mouse models. IGFN1 was initially characterized as a protein fragment using a yeast two-hybrid assay, with the KY protein as bait (32). IGFN1-deficient clones exhibit a reduced fusion index and pronounced morphological differences compared to C2C12 myotubes, indicating IGFN1’s role in myoblast fusion and differentiation (33). Immunohistochemistry confirmed elevated IGFN1 expression in the diabetic ulcer model compared to the normal, suggesting that IGFN1 may be a critical target in diabetic ulcers. Consequently, further exploring the pathways through which IGFN1 influences these changes could have significant implications for diabetic ulcer prognosis.

Proteomics offers a powerful approach for elucidating the molecular mechanisms involved in disease progression, utilizing high-throughput, precise, sensitive, and reproducible techniques. TMT coupled with LC-MS/MS enables effective protein identification and quantification in tissue samples (34). In our proteomic analysis, the results showed significant enrichment of proteins related to DCM, oxidative phosphorylation, the PPAR signaling pathway, and the HIF-1 signaling pathway, suggesting the potential roles of these proteins in diabetic ulcers. Advanced DCM is often accompanied by severe PAD. However, effective therapeutic targets or treatments remain lacking for patients with these conditions. Therefore, discovering effective targets for DFUs may also improve DCM prognosis and treatment, representing important clinical significance. Diabetes-induced oxidative stress continually damages endothelial cells and impairs wound healing (35). The significance of HIF-1 inhibition, first noted in diabetic wound healing (36), is associated with compromised wound healing in diabetes. Enhancing HIF-1 activity has been shown to facilitate wound healing by promoting angiogenesis, proliferation, and migration of fibroblasts in diabetic mouse models (37, 38). Deferoxamine, an iron-chelating agent commonly used for treating iron overload, can mitigate oxidative stress and activate HIF-1, thus expediting diabetic wound healing (37). A recently optimized topical drug delivery system is set to undergo clinical trials (ClinicalTrials.gov registration no. NCT03137966) to evaluate its effectiveness in patients with DFUs. Notably, while HIF-1 signaling plays an important role in diabetic ulcers, a direct link between the PPAR signaling pathway and diabetic ulcers remains unexplored, presenting an avenue for future research.

Recent research has increasingly recognized that diabetes mellitus (DM) and its complications involve profound metabolic disturbances, which may be driven directly by hyperglycemia or occur independently through complex regulatory networks (39, 40). As an integrative reflection of both endogenous and exogenous factors—including genomic, transcriptomic, and proteomic layers—metabolomics offers unique insights into disease-specific molecular phenotypes. In particular, metabolic reprogramming has emerged as a critical interface between energy homeostasis and immune-inflammatory responses (41, 42). In our metabolomic analysis of DFUs, we observed significant enrichment in pathways related to amino acid biosynthesis, notably involving arginine, glutamine, and tryptophan. These amino acids are increasingly recognized as immunomodulatory metabolites. Arginine, for instance, serves as a substrate for arginase-1 (Arg1), promoting M2 macrophage polarization and suppressing chronic inflammation via inhibition of NF-κB–mediated secretion of proinflammatory cytokines such as TNF-α and IL-6 (43). Glutamine, a key energy source for lymphocytes and neutrophils, modulates immune function through the mTOR signaling axis; its depletion has been associated with impaired immune responses and enhanced inflammatory cascades (44). Tryptophan catabolism yields kynurenine, a ligand for the aryl hydrocarbon receptor (AhR), which has been shown to suppress NLRP3 inflammasome activation and mitigate oxidative stress (45).

Consistent with these findings, our data revealed decreased glutamine levels and increased kynurenine accumulation in DFU samples, suggesting that altered amino acid metabolism contributes to immune dysregulation and chronic inflammation in diabetic wound microenvironments. Moreover, metabolic stress induced by hyperglycemia and insulin resistance leads to mitochondrial dysfunction and shifts in cellular energy demands. Under such stress, hypoxia-inducible factor-1α (HIF-1α) is activated, promoting the transcription of glycolytic enzymes (e.g., LDHA) while simultaneously repressing oxidative phosphorylation (46). Notably, our pathway analysis revealed concurrent enrichment of HIF-1 signaling and amino acid biosynthesis pathways, implying a coordinated metabolic adaptation that enables immune cells to function in hypoxic, nutrient-deprived settings characteristic of chronic wounds. Taken together, these results highlight the central role of amino acid metabolic reprogramming in orchestrating inflammatory and energy responses during DFU pathogenesis.

Recent studies have further revealed the dynamic characteristics of metabolic reprogramming in DFU. Related studies have shown that metabolomics has found significant accumulation of phenylpyruvate in DFUs. Increased phenylpyruvate impairs wound healing and enhances inflammatory responses, while reducing phenylpyruvate by restricting phenylalanine in the diet can alleviate uncontrolled inflammation and benefit diabetic wounds (47). There are also studies showing that TMAO accumulates in diabetic wounds, causing macrophages to persist in the pro-inflammatory M1 phenotype, hindering tissue repair. TMAO-targeted anti-inflammatory agents can reprogram macrophage metabolism by inhibiting the IRE1α/XBP1/HIF-1α signaling pathway, prompting it to switch to anti-inflammatory M2, while improving angiogenesis and reducing oxidative stress (48). These advances together support the core conclusion of this study—the metabolism-inflammatory axis is a key target for the treatment of DFU.

Meanwhile, results from the NHANES cohort revealed a significant positive association between TyG-WC/BMI and the risk of non-healing lower extremity ulcers (NHLU). TyG-WC and TyG-BMI have been recognized as sensitive indicators of insulin resistance and visceral obesity, with elevated levels reflecting lipid metabolic disorders and a state of chronic low-grade inflammation (49, 50). Accumulation of visceral adipose tissue has been reported to promote the secretion of pro-inflammatory cytokines such as IL-6, thereby activating the NF-κB signaling pathway, impairing fibroblast migration and angiogenesis, and ultimately delaying ulcer healing (51–53). n our diabetic mouse model, we observed significant upregulation of NF-κB–related genes, including TNF and IL1B, consistent with the inflammatory phenotype observed in patients with elevated TyG-WC/BMI. These findings suggest that TyG-WC/BMI may serve as a potential biomarker for identifying individuals at high risk for NF-κB pathway overactivation and may aid in guiding the use of anti-inflammatory therapeutic strategies, such as JAK/STAT inhibitors.

This study also has a few limitations, the study relies on diabetic ulcer mouse models, which, while valuable, may not fully replicate the complexity of diabetic ulcers in humans, limiting the direct applicability of findings.

5 Conclusion

This study emphasizes the complexity and serious implications of diabetic ulcers, highlighting their association with significant comorbidities and elevated mortality rates. The integration of high-throughput multi-omics techniques—transcriptomics, proteomics, and metabolomics—enabled the development of a compound-reaction-enzyme-gene network that elucidates the underlying molecular mechanisms of diabetic ulcers. This comprehensive molecular map provides valuable insights into diabetic ulcer pathogenesis and facilitates the identification of novel therapeutic targets and strategies for improving patient outcomes.

Data availability statement

The contributions originally presented in this study are included in the article and supplementary material. Further inquiries can be directed to the corresponding authors.

Ethics statement

This study uses publicly available data from the NHANES database, which has been ethically approved with patient informed consent. All experimental procedures involving animals were conducted in full compliance with the Guide for the Care and Use of Laboratory Animals (NIH Publication No. 85-23, revised 1996) and received approval from the Laboratory Animal Center at the First Affiliated Hospital of Shandong First Medical University (Approval number: SYDWLS [2021] 002). Our study followed the ARRIVE guidelines for reporting animal research. Meanwhile, this study was conducted in accordance with the 1964 Helsinki Declaration and approved by the First Affiliated Hospital of Shandong First Medical University’s Institutional Review Board (IRB) (YXLL-KY-2022 (013)). In this study, the samples were utilized exclusively for their intended scientific purposes, without additional collection or any impact on patient privacy.

Author contributions

YL: Data curation, Formal analysis, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. HG: Data curation, Formal analysis, Investigation, Methodology, Software, Supervision, Writing – original draft, Writing – review & editing. SW: Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Writing – original draft, Writing – review & editing. HX: Methodology, Resources, Supervision, Validation, Visualization, Writing – review & editing. ZC: Project administration, Resources, Validation, Visualization, Writing – review & editing. YZ: Conceptualization, Investigation, Methodology, Project administration, Supervision, Visualization, Writing – original draft, Writing – review & editing. YG: Conceptualization, Funding acquisition, Project administration, Supervision, Validation, Writing – original draft, Writing – review & editing. XS: Supervision, Project administration, Funding acquisition, Resources, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This work was supported by the National Natural Science Foundation of China (82104876).

Acknowledgments

We would like to thank all study participants as well as all investigators of the studies. And we would like to thank Editage (www.editage.cn) for English language editing.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2025.1574858/full#supplementary-material

Abbreviations

DFUs, diabetic foot ulcers; DEPs, differentially expressed proteins; DM, Diabetes mellitus; PN, peripheral neuropathy; PAD, peripheral arterial disease; NHANES, The National Health and Nutrition Examination Survey; BMI, body mass index; WC, waist circumference; NHLU, non-healing lower extremity ulcers; TC, total cholesterol; CRP, C-reactive protein; SBP, systolic blood pressure; STZ, streptozotocin; RNA seq, RNA Sequencing; UHPLC, ultra-high-performance liquid chromatography; PCA, principal component analysis; FBG, fasting blood glucose; HbA1c, hemoglobin A1c; T2DM, type 2 diabetes mellitus; DEGs, differentially expressed genes; SNPs, single nucleotide polymorphisms; FC, fold change; PLS-DA, Partial Least Squares Discriminant Analysis; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

References

1. Cho NH, Shaw JE, Karuranga S, Huang Y, da Rocha Fernandes JD, Ohlrogge AW, et al. IDF Diabetes Atlas: Global estimates of diabetes prevalence for 2017 and projections for 2045. Diabetes Res Clin Pract. (2018) 138:271–81. doi: 10.1016/j.diabres.2018.02.023

PubMed Abstract | Crossref Full Text | Google Scholar

2. Tuttolomondo A, Maida C, and Pinto A. Diabetic foot syndrome as a possible cardiovascular marker in diabetic patients. J Diabetes Res. (2015) 2015:268390. doi: 10.1155/2015/268390

PubMed Abstract | Crossref Full Text | Google Scholar

3. Fryberg RG, Armstrong DG, Giurini J, Edwards A, Kravette M, and Kravitz S. Diabetic foot disorders: a clinical practice guideline. American College of Foot and Ankle Surgeons. J Foot Ankle Surg. (2000) 39(5):S1-60. doi: 10.1016/S1067-2516(07)60001-5

PubMed Abstract | Crossref Full Text | Google Scholar

4. Armstrong DG, Swerdlow MA, Armstrong AA, Conte MS, Padula WV, and Bus SA. Five year mortality and direct costs of care for people with diabetic foot complications are comparable to cancer. J Foot Ankle Res. (2020) 13:16. doi: 10.1186/s13047-020-00383-2

PubMed Abstract | Crossref Full Text | Google Scholar

5. Jeyaraman K, Berhane T, Hamilton M, Chandra AP, and Falhammar H. Mortality in patients with diabetic foot ulcer: a retrospective study of 513 cases from a single Centre in the Northern Territory of Australia. BMC Endocr Disord. (2019) 19:1. doi: 10.1186/s12902-018-0327-2

PubMed Abstract | Crossref Full Text | Google Scholar

6. Tuttolomondo A, Maida C, and Pinto A. Diabetic foot syndrome: Immune-inflammatory features as possible cardiovascular markers in diabetes. World J Orthop. (2015) 6:62–76. doi: 10.5312/wjo.v6.i1.62

PubMed Abstract | Crossref Full Text | Google Scholar

7. Robles-Remacho A, Sanchez-Martin RM, and Diaz-Mochon JJ. Spatial transcriptomics: emerging technologies in tissue gene expression profiling. Analytical Chem. (2023) 95:15450–60. doi: 10.1021/acs.analchem.3c02029

PubMed Abstract | Crossref Full Text | Google Scholar

8. Hanash S. Disease proteomics. Nature. (2003) 422:226–32. doi: 10.1038/nature01514

PubMed Abstract | Crossref Full Text | Google Scholar

9. Fraga-Corral M, Carpena M, Garcia-Oliveira P, Pereira AG, Prieto MA, and Simal-Gandara J. Analytical metabolomics and applications in health, environmental and food science. Crit Rev Analytical Chem. (2020) 52:712–34. doi: 10.1080/10408347.2020.1823811

PubMed Abstract | Crossref Full Text | Google Scholar

10. Montaner J, Ramiro L, Simats A, Tiedt S, Makris K, Jickling GC, et al. Multilevel omics for the discovery of biomarkers and therapeutic targets for stroke. Nat Rev Neurol. (2020) 16:247–64. doi: 10.1038/s41582-020-0350-6

PubMed Abstract | Crossref Full Text | Google Scholar

11. Hong S, Han K, and Park CY. The triglyceride glucose index is a simple and low-cost marker associated with atherosclerotic cardiovascular disease: a population-based study. BMC Med. (2020) 18:361. doi: 10.1186/s12916-020-01824-2

PubMed Abstract | Crossref Full Text | Google Scholar

12. Guerrero-Romero F, Villalobos-Molina R, Jimenez-Flores JR, Simental-Mendia LE, Mendez-Cruz R, Murguia-Romero M, et al. Fasting triglycerides and glucose index as a diagnostic test for insulin resistance in young adults. Arch Med Res. (2016) 47:382–7. doi: 10.1016/j.arcmed.2016.08.012

PubMed Abstract | Crossref Full Text | Google Scholar

13. Shi W, Xing L, Jing L, Tian Y, Yan H, Sun Q, et al. Value of triglyceride-glucose index for the estimation of ischemic stroke risk: Insights from a general population. Nutr Metab Cardiovasc Dis. (2020) 30:245–53. doi: 10.1016/j.numecd.2019.09.015

PubMed Abstract | Crossref Full Text | Google Scholar

14. Yan H, Zhou Q, Wang Y, Tu Y, Zhao Y, Yu J, et al. Associations between cardiometabolic indices and the risk of diabetic kidney disease in patients with type 2 diabetes. Cardiovasc Diabetol. (2024) 23:142. doi: 10.1186/s12933-024-02228-9

PubMed Abstract | Crossref Full Text | Google Scholar

15. Zhang Q, Xiao S, Jiao X, and Shen Y. The triglyceride-glucose index is a predictor for cardiovascular and all-cause mortality in CVD patients with diabetes or pre-diabetes: evidence from NHANES 2001-2018. Cardiovasc Diabetol. (2023) 22:279. doi: 10.1186/s12933-023-02030-z

PubMed Abstract | Crossref Full Text | Google Scholar

16. Wisniewski JR, Zougman A, Nagaraj N, and Mann M. Universal sample preparation method for proteome analysis. Nat Methods. (2009) 6:359–62. doi: 10.1038/nmeth.1322

PubMed Abstract | Crossref Full Text | Google Scholar

17. del Toro N, Shrivastava A, Ragueneau E, Meldal B, Combe C, Barrera E, et al. The IntAct database: efficient access to fine-grained molecular interaction data. Nucleic Acids Res. (2022) 50:D648–D53. doi: 10.1093/nar/gkab1006

PubMed Abstract | Crossref Full Text | Google Scholar

18. Yu CS, Lin CJ, and Hwang JK. Predicting subcellular localization of proteins for Gram-negative bacteria by support vector machines based on n-peptide compositions. Protein Sci. (2004) 13:1402–6. doi: 10.1110/ps.03479604

PubMed Abstract | Crossref Full Text | Google Scholar

19. Gu Z, Li L, Tang S, Liu C, Fu X, Shi Z, et al. Metabolomics Reveals that Crossbred Dairy Buffaloes Are More Thermotolerant than Holstein Cows under Chronic Heat Stress. J Agric Food Chem. (2018) 66:12889–97. doi: 10.1021/acs.jafc.8b02862

PubMed Abstract | Crossref Full Text | Google Scholar

20. Luo D, Deng T, Yuan W, Deng H, and Jin M. Plasma metabolomic study in Chinese patients with wet age-related macular degeneration. BMC Ophthalmol. (2017) 17:165. doi: 10.1186/s12886-017-0555-7

PubMed Abstract | Crossref Full Text | Google Scholar

21. Shou Y, Le Z, Cheng HS, Liu Q, Ng YZ, Becker DL, et al. Mechano-activated cell therapy for accelerated diabetic wound healing. Adv Mater. (2023) 35:e2304638. doi: 10.1002/adma.202304638

PubMed Abstract | Crossref Full Text | Google Scholar

22. Tu C, Lu H, Zhou T, Zhang W, Deng L, Cao W, et al. Promoting the healing of infected diabetic wound by an anti-bacterial and nano-enzyme-containing hydrogel with inflammation-suppressing, ROS-scavenging, oxygen and nitric oxide-generating properties. Biomaterials. (2022) 286:121597. doi: 10.1016/j.biomaterials.2022.121597

PubMed Abstract | Crossref Full Text | Google Scholar

23. Deng C, Zheng M, Han S, Wang Y, Xin J, Aras O, et al. GSH-activated porphyrin sonosensitizer prodrug for fluorescence imaging-guided cancer sonodynamic therapy. Adv Funct Mater. (2023) 33:2300348. doi: 10.1002/adfm.202300348

PubMed Abstract | Crossref Full Text | Google Scholar

24. Li Z, Zhao Y, Huang H, Zhang C, Liu H, Wang Z, et al. A nanozyme-immobilized hydrogel with endogenous ROS-scavenging and oxygen generation abilities for significantly promoting oxidative diabetic wound healing. Adv Healthc Mater. (2022) 11:e2201524. doi: 10.1002/adhm.202201524

PubMed Abstract | Crossref Full Text | Google Scholar

25. Chen L, Zheng B, Xu Y, Sun C, Wu W, Xie X, et al. Nano hydrogel-based oxygen-releasing stem cell transplantation system for treating diabetic foot. J Nanobiotechnology. (2023) 21:202. doi: 10.1186/s12951-023-01925-z

PubMed Abstract | Crossref Full Text | Google Scholar

26. Sharma S, Schaper N, and Rayman G. Microangiopathy: Is it relevant to wound healing in diabetic foot disease? Diabetes Metab Res Rev. (2020) 36 Suppl 1:e3244. doi: 10.1002/dmrr.3244

PubMed Abstract | Crossref Full Text | Google Scholar

27. Teh HX, Phang SJ, Looi ML, Kuppusamy UR, and Arumugam B. Molecular pathways of NF-kB and NLRP3 inflammasome as potential targets in the treatment of inflammation in diabetic wounds: A review. Life Sci. (2023) 334:122228. doi: 10.1016/j.lfs.2023.122228

PubMed Abstract | Crossref Full Text | Google Scholar

28. Fei J, Ling YM, Zeng MJ, and Zhang KW. Shixiang Plaster, a Traditional Chinese Medicine, Promotes Healing in a Rat Model of Diabetic Ulcer Through the receptor for Advanced Glycation End Products (RAGE)/Nuclear Factor kappa B (NF-kappaB) and Vascular Endothelial Growth Factor (VEGF)/Vascular Cell Adhesion Molecule-1 (VCAM-1)/Endothelial Nitric Oxide Synthase (eNOS) Signaling Pathways. Med Sci Monit. (2019) 25:9446–57. doi: 10.12659/MSM.918268

PubMed Abstract | Crossref Full Text | Google Scholar

29. Pandith H, Zhang X, Thongpraditchote S, Wongkrajang Y, Gritsanapan W, and Baek SJ. Effect of Siam weed extract and its bioactive component scutellarein tetramethyl ether on anti-inflammatory activity through NF-kappaB pathway. J Ethnopharmacol. (2013) 147:434–41. doi: 10.1016/j.jep.2013.03.033

PubMed Abstract | Crossref Full Text | Google Scholar

30. Jang DI, Lee AH, Shin HY, Song HR, Park JH, Kang TB, et al. The role of tumor necrosis factor alpha (TNF-alpha) in autoimmune disease and current TNF-alpha inhibitors in therapeutics. Int J Mol Sci. (2021) 22:2719. doi: 10.3390/ijms22052719

PubMed Abstract | Crossref Full Text | Google Scholar

31. Ban E, Jeong S, Park M, Kwon H, Park J, Song EJ, et al. Accelerated wound healing in diabetic mice by miRNA-497 and its anti-inflammatory activity. BioMed Pharmacother. (2020) 121:109613. doi: 10.1016/j.biopha.2019.109613

PubMed Abstract | Crossref Full Text | Google Scholar

32. Beatham J, Romero R, Townsend SK, Hacker T, van der Ven PF, and Blanco G. Filamin C interacts with the muscular dystrophy KY protein and is abnormally distributed in mouse KY deficient muscle fibres. Hum Mol Genet. (2004) 13:2863–74. doi: 10.1093/hmg/ddh308

PubMed Abstract | Crossref Full Text | Google Scholar

33. Li X, Baker J, Cracknell T, Haynes AR, and Blanco G. IGFN1_v1 is required for myoblast fusion and differentiation. PLoS One. (2017) 12:e0180217. doi: 10.1371/journal.pone.0180217

PubMed Abstract | Crossref Full Text | Google Scholar

34. Jeffcoate WJ, Macfarlane RM, and Fletcher EM. The description and classification of diabetic foot lesions. Diabetes Med. (1993) 10:676–9. doi: 10.1111/j.1464-5491.1993.tb00144.x

PubMed Abstract | Crossref Full Text | Google Scholar

35. Wilkinson HN and Hardman MJ. Wound healing: cellular mechanisms and pathological outcomes. Open Biol. (2020) 10:200223. doi: 10.1098/rsob.200223

PubMed Abstract | Crossref Full Text | Google Scholar

36. Botusan IR, Sunkari VG, Savu O, Catrina AI, Grünler J, Lindberg S, et al. Stabilization of HIF-1alpha is critical to improve wound healing in diabetic mice. Proc Natl Acad Sci U S A. (2008) 105:19426–31. doi: 10.1073/pnas.0805230105

PubMed Abstract | Crossref Full Text | Google Scholar

37. Duscher D, Neofytou E, Wong VW, Maan ZN, Rennert RC, Inayathullah M, et al. Transdermal deferoxamine prevents pressure-induced diabetic ulcers. Proc Natl Acad Sci U S A. (2015) 112:94–9. doi: 10.1073/pnas.1413445112

PubMed Abstract | Crossref Full Text | Google Scholar

38. Zhu Y, Wang Y, Jia Y, Xu J, and Chai Y. Roxadustat promotes angiogenesis through HIF-1alpha/VEGF/VEGFR2 signaling and accelerates cutaneous wound healing in diabetic rats. Wound Repair Regen. (2019) 27:324–34. doi: 10.1111/wrr.12708

PubMed Abstract | Crossref Full Text | Google Scholar

39. Pereira PR, Carrageta DF, Oliveira PF, Rodrigues A, Alves MG, and Monteiro MP. Metabolomics as a tool for the early diagnosis and prognosis of diabetic kidney disease. Med Res Rev. (2022) 42:1518–44. doi: 10.1002/med.21883

PubMed Abstract | Crossref Full Text | Google Scholar

40. Pillon NJ, Loos RJF, Marshall SM, and Zierath JR. Metabolic consequences of obesity and type 2 diabetes: Balancing genes and environment for personalized care. Cell. (2021) 184:1530–44. doi: 10.1016/j.cell.2021.02.012

PubMed Abstract | Crossref Full Text | Google Scholar

41. O’Neill LA and Pearce EJ. Immunometabolism governs dendritic cell and macrophage function. J Exp Med. (2016) 213:15–23. doi: 10.1084/jem.20151570

PubMed Abstract | Crossref Full Text | Google Scholar

42. SantaCruz-Calvo S, Bharath L, Pugh G, SantaCruz-Calvo L, Lenin RR, Lutshumba J, et al. Adaptive immune cells shape obesity-associated type 2 diabetes mellitus and less prominent comorbidities. Nat Rev Endocrinol. (2022) 18:23–42. doi: 10.1038/s41574-021-00575-1

PubMed Abstract | Crossref Full Text | Google Scholar

43. Rath M, Muller I, Kropf P, Closs EI, and Munder M. Metabolism via arginase or nitric oxide synthase: two competing arginine pathways in macrophages. Front Immunol. (2014) 5:532. doi: 10.3389/fimmu.2014.00532

PubMed Abstract | Crossref Full Text | Google Scholar

44. Cruzat V, Macedo Rogero M, Noel Keane K, Curi R, and Newsholme P. Glutamine: metabolism and immune function, supplementation and clinical translation. Nutrients. (2018) 10:1564. doi: 10.3390/nu10111564

PubMed Abstract | Crossref Full Text | Google Scholar

45. Platten M, Nollen EAA, Rohrig UF, Fallarino F, and Opitz CA. Tryptophan metabolism as a common therapeutic target in cancer, neurodegeneration and beyond. Nat Rev Drug Discov. (2019) 18:379–401. doi: 10.1038/s41573-019-0016-5

PubMed Abstract | Crossref Full Text | Google Scholar

46. Semenza GL. HIF-1 mediates metabolic responses to intratumoral hypoxia and oncogenic mutations. J Clin Invest. (2013) 123:3664–71. doi: 10.1172/JCI67230

PubMed Abstract | Crossref Full Text | Google Scholar

47. Lv D, Cao X, Zhong L, Dong Y, Xu Z, Rong Y, et al. Targeting phenylpyruvate restrains excessive NLRP3 inflammasome activation and pathological inflammation in diabetic wound healing. Cell Rep Med. (2023) 4:101129. doi: 10.1016/j.xcrm.2023.101129

PubMed Abstract | Crossref Full Text | Google Scholar

48. Lu L, Liao J, Xu C, Xiong Y, Zhou J, Wang G, et al. Kinsenoside-loaded microneedle accelerates diabetic wound healing by reprogramming macrophage metabolism via inhibiting IRE1alpha/XBP1 signaling axis. Adv Sci (Weinh). (2025):e2502293:1-17. doi: 10.1002/advs.202502293

PubMed Abstract | Crossref Full Text | Google Scholar

49. Wan X, Ji Y, Wang R, Yang H, Cao X, and Lu S. Association between systemic immune-inflammation index and sarcopenic obesity in middle-aged and elderly Chinese adults: a cross-sectional study and mediation analysis. Lipids Health Dis. (2024) 23:230. doi: 10.1186/s12944-024-02215-9

PubMed Abstract | Crossref Full Text | Google Scholar

50. Huang Y, Zhou Y, Xu Y, Wang X, Zhou Z, Wu K, et al. Inflammatory markers link triglyceride-glucose index and obesity indicators with adverse cardiovascular events in patients with hypertension: insights from three cohorts. Cardiovasc Diabetol. (2025) 24:11. doi: 10.1186/s12933-024-02571-x

PubMed Abstract | Crossref Full Text | Google Scholar

51. Barroso WA, Victorino VJ, Jeremias IC, Petroni RC, Ariga SKK, Salles TA, et al. High-fat diet inhibits PGC-1alpha suppressive effect on NFkappaB signaling in hepatocytes. Eur J Nutr. (2018) 57:1891–900. doi: 10.1007/s00394-017-1472-5

PubMed Abstract | Crossref Full Text | Google Scholar

52. Fain JN, Bahouth SW, and Madan AK. Involvement of multiple signaling pathways in the post-bariatric induction of IL-6 and IL-8 mRNA and release in human visceral adipose tissue. Biochem Pharmacol. (2005) 69:1315–24. doi: 10.1016/j.bcp.2005.02.009

PubMed Abstract | Crossref Full Text | Google Scholar

53. Shan B, Shao M, Zhang Q, Hepler C, Paschoal VA, Barnes SD, et al. Perivascular mesenchymal cells control adipose-tissue macrophage accrual in obesity. Nat Metab. (2020) 2:1332–49. doi: 10.1038/s42255-020-00301-7

PubMed Abstract | Crossref Full Text | Google Scholar