The drug targets of Evodia rutaecarpa were studied via the Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP) (https://tcmspw.com/tcmsp.php) (14). In addition, the OMIM database (https://www.omim.org/) (15) and GeneCards database (https://www.genecards.org/) (16) were used to screen the pathological targets of DFU. A Venn diagram was used to integrate the drug targets of Chinese medicine and the pathological targets of DFU to determine the potential mechanism of action of traditional Chinese medicine in the treatment of DFU.

The series matrix file of GSE147890 was downloaded from the NCBI GEO public database (https://www.ncbi.nlm.nih.gov/geo/info/datasets.html) (17), and the annotation file was GPL571. A total of 24 samples of expression spectrum data were included, including 12 cases in the control group and 12 cases in the disease group. The GSE165816 dataset from the NCBI GEO public database was downloaded, and a total of 22 datasets were included for single-cell analysis.

Consensus clustering was performed on the basis of the drug–disease intersection gene expression matrix to define DFU subtypes. The analysis involved 50 iterative computations, each utilizing a random 80% subset of the samples. The determination of the optimal cluster count was informed by the characteristics of the consensus matrix heatmap and the trajectory of the consensus cumulative distribution function curve (18).

To investigate signalling pathway heterogeneity between the identified subtypes, GSEA was conducted. The background gene set was the version 7.0 annotation gene set downloaded from the MsigDB database (https://www.gsea-msigdb.org/gsea/msigdb) as the annotation gene set for the subtype pathways. Gene sets demonstrating significant enrichment (adjusted P < 0.05) in the intersubtype differential expression analysis were subsequently ranked on the basis of their normalized enrichment scores (19).

LASSO is a compression estimate. It constructs a penalty function to obtain a more refined model, which compresses some coefficients and sets some coefficients to zero. Therefore, the advantages of subset shrinkage are retained, and it is a biased estimate for processing data with complex collinearity. We used the Lasso logistic regression algorithm to select features for diagnostic markers of diseases. The Lasso algorithm uses the “glmnet” package.

The R package “ClusterProfiler” was used to perform functional annotations on important genes to explore the functional relevance of these genes fully. KEGG was used to evaluate relevant functional categories. KEGG-enriched pathways with p values and q values less than 0.05 were considered significant categories.

First, the expression profile was read through the Seurat package, and low-expression genes were screened out (nFeature_RNA > 200 & percent.mt <= median+3MAD & nFeature_RNA <= median+3MAD & nCount_RNA <= median+3MAD & percent.ribo <= median+3MAD); the data were standardized, homogenized, and subjected to principal component analysis (PCA). In turn, the optimal number of pcs was observed through ElbowPlot, and the positional relationship between each cluster was obtained through UMAP analysis. The clusters were annotated by marker genes, and some cells that were important for the occurrence of the disease were annotated.

For the key genes, the corresponding protein 3D structure was obtained from the RCSB PDB database (https://www.rcsb.org/), and the drug component structure was obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). Molecular docking was performed via PyMOL software and AutoDock software.

Statistical analysis was performed via R language (version 4.2.2). P < 0.05 was considered statistically significant.

HUVECs and RAW264.7 cells were purchased from Haixing Biotechnology Co., Ltd. (Suzhou, China). HUVECs were cultured in Endothelial Cell Basal Medium supplemented with 5% fetal bovine serum (FBS), 1% Endothelial Cell Culture Supplement, and 1% penicillin/streptomycin. RAW264.7 cells were cultured in Dulbecco’s Modified Eagle’s Medium (5.5 mM D-glucose) supplemented with 10% FBS and 1% penicillin/streptomycin. Both cell lines were incubated at 37°C under 5% CO₂ and cultured separately.

RAW264.7 cells and HUVECs were seeded into 96-well plates, respectively. After 24-h incubation with various concentrations (0, 6.25, 12.5, 25, 50, 100 μM) of Que, 50 μL of 3-(4, 5)-dimethylthiahiazo(-z-y1)-2,5-diphenyltetrazolium bromide (MTT) (Solarbio, Beijing, China) (2 mg/mL) was added per well and then incubated for 3 h. The formazan crystals were dissolved in 200 μL of dimethyl sulfoxide (DMSO) and the absorbance was read at 490 nm using a microplate reader.

HUVECs in the logarithmic growth phase were pretreated with Que or DMSO at 37°C for 30 min. The samples were then aliquoted and subjected to a temperature gradient (38, 42, 48, 54, 62, 68°C) using a thermal cycler, followed by immediate cooling on ice. After centrifugation at 12,000 g and 4°C for 15 min, the protein expression level of hypoxia-inducible factor-1 α (HIF1α) (Wanleibio, Shenyang, China) was quantified by western blotting.

Based on a previously described method (20), a co-culture system was established using RAW264.7 cells and HUVECs. Specifically, HUVECs were seeded in the lower chamber of Transwell inserts, while RAW264.7 cells were placed in the upper chamber. The co−culture system was treated with LPS (100 ng/mL) or varying concentrations of Que for 24 h. After incubation, HUVECs from the lower chamber were collected, and the protein expression levels of HIF1α and vascular endothelial growth factor (VEGF) (Wanleibio, Shenyang, China) were evaluated by western blotting.

At the endpoint of the 7-day or 14-day intervention period, the rats were first anesthetized with 4% isoflurane for induction, followed by maintenance of anesthesia with 2% isoflurane. After serum collection, the isoflurane concentration was increased to 5% to ensure humane euthanasia and minimize animal suffering as much as possible. Serum levels of insulin (INS), low-density lipoprotein cholesterol (LDL-C), triglyceride (TG), total cholesterol (T-CHO), tumor necrosis factor α (TNF-α), and interleukin-1β (IL-1β) were assessed according to the manufacturers’ instructions for commercial kits (Jiancheng, Nanjing, China).

The wounds were photographed via a digital camera, and the wound healing rate was calculated by measuring the wound area with ImageJ software. The reduction in wound size was calculated via the following equation: reduction in wound size (%) = (A₀ − A_t)/A₀ × 100%, where A₀ is the initial wound area and A_t is the wound area on days 0, 3, 7, 10 and 14 after surgery. The wound area at 7 d and 14 d was considered an important indicator of wound healing (26). Wound tissues were collected on days 7 and 14 after intervention, preserved in 4% paraformaldehyde solution, dehydrated, and processed with paraffin. The solidified paraffin-embedded samples were then subjected to serial sectioning. Histopathological evaluation was performed via hematoxylin–eosin staining and Masson’s trichrome staining.

Following collection into sterile cryovials from each experimental group, the preservation of fecal samples involved immediate flash freezing in liquid nitrogen (1 h) prior to transfer to a -80°C freezer. For characterization of the gut microbial community, 16S rRNA gene sequencing was performed by Biomarker Technologies Co., Ltd. (Beijing, China) (21).

After the intestinal tissues of the rats were lysed, the protein concentration was quantified via a BCA assay. The expression levels of Zonula occludens-1 (ZO-1), Occludin and Claudin 1 (Proteintech, Wuhan, China) were detected via western blotting.

The wound images of each group were taken and recorded at 0, 3, 7, and 14 days. ImageJ software was used to measure and calculate the wound healing rate. Wound tissues from the rats were collected for histological analysis via hematoxylin–eosin staining, Masson’s trichrome staining, and cluster of differentiation 31 (CD31) immunohistochemical staining, followed by microscopic examination to evaluate remodelling and pathological changes. RNA was also extracted from these tissues to analyze the mRNA expression levels of VEGF, HIF1α, and IL-6 (Table 4). The expression levels of VEGF and HIF1α were detected via western blotting.

To investigate the potential impact of Que on the progression of chronic diabetic wound healing, we used type 1 diabetic SD rats to further evaluate wound healing rates and quality (Figure 6A). As shown in Figures 6B, C, on the 7th and 14th days, compared with the sham group, the model group exhibited a significant decrease in body weight, and the upper plasma presented a milky white and turbid appearance. In contrast, Que intervention attenuated the decrease in body weight of diabetic rats, and the plasma was clear and pale yellow. Oral glucose tolerance test data revealed a rapid increase in serum glucose levels in all groups following glucose administration, peaking at approximately 30 min, followed by a gradual decrease. Compared with the sham group, the model group presented significantly elevated serum glucose levels at all time-matched time points (Figure 6D). Que treatment also significantly reduced T-CHO, TG, and LDL-C levels (Figures 6E-G) while increasing INS concentrations (Figure 6H). Compared with those in the untreated model groups, the levels of inflammatory cytokines, particularly TNF-α and IL-1β, were markedly lower in the Que-treated groups (Figures 6I, J).

Figures 7A-C shows that, compared with that in the sham group, the wound healing process in the model group was prolonged, and the wound healing rate was lower. After treatment with Que, the wound healing rates improved within 14 days. To investigate the potential advantages of Que at different stages of chronic diabetic wound healing, we focused on detecting neovascularization within 7 days and collagen fibre deposition within 14 days. On day 7, diabetic wounds treated with Que displayed an increased relative number of of neovessels, occasional granulation tissue with new blood vessels, and collagen formation characterized by well-arranged fibroblasts. On day 14, the collagen fibres in the Que group were not only densely packed but also regularly arranged, indicating the initiation of the remodelling phase (Figure 7D). Masson staining further confirmed that, compared with the model group, the newly synthesized collagen fibers in the Que group exhibited a significantly higher density and more ordered structure (Figures 7E, F).

To further interpret the changes in metabolites in the serum of the rats in the disease state, the serum of the rats in each group on day 14 was analyzed via untargeted metabolomics. Spearman’s rank correlation coefficient (r > 0.8) was obtained for all pairwise sample comparisons, indicating high methodological consistency and robust biological reproducibility (Figure 8A). PCA revealed distinct segregation of the model group from both the sham and Que-treated groups but partial overlap with the metformin intervention group. Significant differences were observed between the Que-treated group and the model group along both the PC1 and PC2 dimensions (Figure 8B). OPLS-DA modelling revealed that both the sham vs. model comparison and model vs. Que-treated group comparisons presented a Q2Y > 0.5, with permutation tests demonstrating positive slopes in the corresponding validation regression lines. These parameters confirm the statistical validity and biological reliability of the OPLS-DA models for differential metabolites screening (Figures 8C-F). Univariate statistical analysis of all detected metabolites was performed. As demonstrated by volcano plots, significantly differential metabolites were identified between the sham and model groups, as well as between the model and Que-treated groups, under both positive and negative ionization modes (Figures 8G, H). In Supplementary Table S1, metabolites with FC ≥ 2 or ≤ 0.5, VIP ≥ 1, and P < 0.05 in the t test were selected as significantly differential metabolites. Comparative analysis among the sham, model, and Que-treated groups revealed 673 differentially abundant metabolites with significant variations in relative abundance. These included 334 metabolites detected in negative ionization mode and 339 in positive ionization mode, including 313 upregulated and 360 downregulated species.

Owing to the important role of the gut bacteria in disease, the fecal samples of the rats were analyzed. A total of 1464722 CCS sequences were retrieved from 24 samples after sequencing and barcode identification (Table 6), and the alpha diversity indices were calculated via Mothur version 1.30 (Table 7). According to the results of the analysis, the sequencing depth of each sample was adequate and saturated, as shown by the Shannon and rarefaction analyses. Most of the microbial diversity information and the structure of the bacterial communities were reflected (Figures 9A, B). The Venn diagrams constructed on the basis of distinct groupings revealed that the four-way intersection contained 797 shared features (Figure 9C). PCA was subsequently performed on the samples, revealing distinct differences among the experimental groups. Although the Quercetin-treated group did not completely overlap with the Sham group in the PCA space, it exhibited a distinct separation from the Model group. This trajectory suggests that Quercetin intervention drives the gut microbiota towards an alternative beneficial equilibrium via structural remodeling, likely by specifically regulating key functional taxa rather than purely reversing the community structure to the naive state (Figure 9D). Compared with the sham group, the model group presented significantly elevated levels of Clostridia, Negativicutes, and Desulfovibrionia at the class taxonomic rank (P < 0.05). Notably, Clostridia and Desulfovibrionia abundances were significantly reduced following Que intervention (P < 0.05) (Figure 9E). At the genus level, the model group presented significantly greater abundances of Lactobacillus, Ruminococcus, unclassified_Lachnospiraceae, uncultured_rumen_bacterium, and unclassified_Ruminococcaceae than did the sham group (P < 0.05), whereas the abundances of Romboutsia, Allobaculum, and Ligilactobacillus were significantly lower (P < 0.05). Following Que intervention, a marked decrease in Romboutsia, uncultured_rumen_bacterium, and unclassified_Ruminococcaceae was observed (P < 0.05), accompanied by a significant increase in Ligilactobacillus abundance (P < 0.05) (Figure 9F). Further functional prediction and annotation of the 16S rRNA sequencing data were performed via PICRUSt2. Functional prediction analysis of the microbiota revealed that Que significantly modulated multiple pathways in diabetic wound rats, including those related to bacterial infections and immune regulation. This further highlights the critical role of inflammatory responses in the disease progression (Figures 9G, H).

On the basis of the aforementioned results, we further conducted an integrated analysis of the microbiota and metabolome by correlating microbial taxa (at both the phylum and genus levels) with metabolomic data from holistic and differential perspectives. Therefore, the differential metabolites-differential microbes chord diagram (Figures 10A, B) and the differential metabolites-differential microbes network diagram (Figures 10C, D) were constructed. The results showed that after Que intervention, Anaerosporobacter was involved in regulating the abnormal metabolism of diseased rats.