Abstract

Background:

Diabetic retinopathy (DR) represents a major microvascular complication arising from diabetes mellitus, characterized by multifactorial pathogenesis encompassing genetic, metabolic, and microbial components. While Dihuang Yinzi (DHYZ), a traditional Chinese herbal formulation, exhibits therapeutic promise for DR management, the precise mechanistic underpinnings warrant further investigation. This research sought to elucidate critical target genes, metabolic compounds, and microbial species implicated in DHYZ's therapeutic action against DR.

Methods:

We constructed a murine DR model and procured biological specimens from 18 animals distributed across three experimental cohorts (control, disease model, and DHYZ-treated groups, n = 6 per group) for metabolomic profiling and microbiome characterization. Candidate genes emerged from overlapping DHYZ-associated targets with DR-linked genes. Through six computational algorithms within the cytoHubba plugin, we pinpointed pivotal target genes. Molecular docking studies examined binding affinity between essential target proteins and bioactive constituents. Metabolomic and microbiome datasets underwent differential expression analysis to enumerate candidate metabolites and microbial taxa, respectively. Finally, Spearman correlation-based integrative omics analysis distinguished critical metabolites and key microbial species.

Results:

We identified 110 candidate genes and five key target genes (STAT3, IL6, TNF, ESR1, and IL1B). Molecular docking analysis revealed strong binding interactions between ESR1 and six corresponding active compounds, with the highest binding affinity observed for naringenin. Additionally, metabolomic analysis identified 50 candidate metabolites, and microbiome analysis revealed 24 candidate microbes. Spearman correlation analysis further pinpointed 30 key metabolites and 18 key microbes.

Conclusion:

This study elucidates five key target genes, 30 key metabolites, and 18 key microbes through which DHYZ may exert its therapeutic effects in DR. These findings provide valuable insights and a foundational reference for understanding the multi-omics mechanism of DHYZ in the treatment of diabetic retinopathy.

1 Introduction

Diabetic retinopathy (DR), a major microvascular complication of diabetes, is the predominant cause of vision loss in the global working-age population. Epidemiological studies estimate that DR affects nearly one-third of individuals with diabetes, corresponding to more than 140 million people worldwide. The risk of developing DR exceeds 60% in those with diabetes duration longer than 10 years (1, 2). According to the latest demographic statistics, approximately 103.12 million people with diabetes were affected by DR worldwide in 2020. This number is projected to increase to 129.84 million by 2030 and 160.50 million by 2045, indicating that more than 57 million new DR cases will be added globally over the next two decades (3, 4). As one of the countries with the heaviest diabetes burden worldwide, China exhibits a particularly prominent epidemiological profile of DR. A study based on the Global Burden of Disease (GBD) 2021 data predicts that the disease burden of DR in China will continue to rise over the next 15 years, which is closely associated with the large population of patients with diabetes and the rapidly aging social structure in China (5).

A substantial diagnosis gap for DR exists globally, with more than 50% of DR cases remaining undiagnosed; this proportion is even higher in developing countries with limited medical resources, such as India (6). Such undiagnosed conditions expose numerous patients to the risk of irreversible visual impairment without any obvious subjective symptoms. Despite advances in ophthalmology, current management of DR faces significant challenges, including inadequate early screening and the limited efficacy of single-target therapeutic interventions, which often result in high recurrence rates. The pathophysiology of DR is multifactorial, involving chronic hyperglycemia-induced endothelial damage in retinal vessels (7), pericyte loss, and disruption of the blood-retinal barrier (8, 9). Consequently, a comprehensive understanding of DR is crucial for developing effective and targeted therapeutic strategies, which could facilitate early detection and improve clinical outcomes for DR patients (10, 11). In summary, the global prevalence of DR continues to expand, with a persistently high proportion of undiagnosed cases and numerous challenges in clinical management, leading to a steadily increasing disease burden. Meanwhile, against the backdrop of uneven distribution of medical resources, the gaps in the prevalence, diagnosis, and treatment of DR across regions and populations have further widened, making it a major global public health issue that urgently requires concerted international attention and response.

To address the complexity of DR pathogenesis, network pharmacology has emerged as a valuable interdisciplinary approach that integrates principles and analytical methods from systems biology and network science (12). It examines complex disease and drug action networks from a systems-level perspective to investigate multi-target pharmacological effects (13, 14). By modeling these intricate interactions, it facilitates the comparative analysis of pathological and drug-response systems, enabling the identification of critical disease drivers and therapeutic targets. Key advantages of this approach include a more comprehensive understanding of drug mechanisms, promotion of novel drug development, support for personalized treatment design, enhancement of treatment efficacy and safety, and efficient discovery of drug targets and modes of action. As a result, network pharmacology provides more accurate and reliable guidance for both drug discovery and clinical translation (15). Given the multi-component and multi-target characteristics of traditional Chinese medicine, network pharmacology represents an ideal approach to systematically investigate the therapeutic mechanisms of herbal formulae such as DHYZ in the treatment of DR.

Dihuang Yinzi (DHYZ) is a traditional Chinese medicinal formula consisting of 15 herbs, with a history of use in the treatment of neurodegenerative diseases such as Alzheimer's disease. The pharmacological actions of its constituent herbs are as follows: Rehmannia glutinosa and Cornus officinalis are known to tonify the kidney and replenish essence; Aconitum carmichaelii, Cinnamomum cassia, Morinda officinalis, and Cistanche deserticola function to warm and invigorate kidney yang; Ophiopogon japonicus, Dendrobium nobile, and Schisandra chinensis work to nourish yin and astringe bodily fluids; Acorus tatarinowii and Polygala tenuifolia help resolve phlegm and open orifices, thereby clearing turbid phlegm from visual pathways; Poria cocos promotes diuresis, eliminates dampness, strengthens the spleen, and harmonizes the stomach; Mentha haplocalyx exhibits antibacterial and anti-inflammatory properties; and Zingiber officinale and Ziziphus jujuba serve to harmonize the stomach, tonify the middle energizer, and regulate gastric qi. Together, these components act synergistically to enhance kidney essence, resolve phlegm, and improve visual function (16, 17). Elucidating the biological mechanisms of DHYZ in DR is of great importance, as it may uncover novel therapeutic targets and pave the way for innovative treatment strategies.

Multi-omics integrated analysis combines biological data from multiple levels—such as metabolomics and microbiomes—enabling researchers to observe gene functions and regulatory mechanisms from diverse perspectives and to elucidate complex molecular interactions and signaling pathways. By leveraging complementary omics datasets, this approach enables in-depth exploration of biological processes, integrating data from phenotypic manifestations to molecular mechanisms, thereby providing novel insights into fundamental biological mechanisms and disease pathogenesis. Applying multi-omics integrated analysis to uncover the complex pathogenesis of DR holds promise for providing critical new insights into its molecular underpinnings and facilitating the development of targeted therapeutic strategies (18).

This study employed network pharmacology to identify key targets of DHYZ in the intervention of DR. Using metabolomic and microbiome sequencing data from control, model, and DHYZ-treated groups, we conducted integrated bioinformatic analyses to identify candidate metabolites and microbes linked to the therapeutic mechanisms of DHYZ across omics layers. Subsequent multi-omics integration highlighted critical metabolites and key microbial taxa, offering novel theoretical and empirical support for future clinical research and treatment development in DR.

2 Materials and methods2.1 Animals

Eighteen male db/db mice at 9 weeks of age, SPF-grade, with mean body weight of 40.88 ± 2.15 g, underwent random allocation into model and DHYZ treatment cohorts via random number table methodology. An additional nine age-matched male db/m mice (mean body weight: 28.39 ± 1.39 g) constituted the control cohort. All experimental subjects were procured from Changzhou Cavens Laboratory Animal Co., Ltd. [License No.: SCXK (Jiangsu) 2021-0013] and housed under specific pathogen-free (SPF) standards at the Animal Experiment Center of Shanxi University of Chinese Medicine. Environmental conditions were maintained at 23 ± 1 °C with 50%−65% relative humidity throughout the study period. This experimental protocol received ethical approval from the Animal Ethics Committee of Shanxi University of Chinese Medicine (Approval No.: 2021DW238).

Following anesthesia via isoflurane inhalation, three mice from each group underwent cardiac perfusion. Perfusion was continued until the liver turned pale, indicating complete fixation. Eyeballs were then enucleated and fixed in 4% paraformaldehyde (PFA) for 24 h. The remaining six mice per group were used for serum and fecal sample collection. These samples were subsequently analyzed using metabolomic profiling and 16S rRNA gene sequencing of the gut microbiota. Corresponding sample identifiers for both metabolomics and Microbiome sequencing are provided in the Supplementary Table.

2.2 Preparation for DHYZ

The corresponding sample identifiers are provided in the Supplementary Table. All herbs were supplied by Beijing Tongrentang Co., Ltd. The mixture was soaked in water for 30 min, decocted twice for 30 min each time, filtered, and then concentrated to a final concentration of 2 g/ml (crude drug). The equivalent dose for mice was calculated based on body surface area conversion. The adult human daily dose of DHYZ is 99 g, corresponding to a dose of 30.03 g/kg for mice in the treatment group. The control and model groups received an equal volume of normal saline. All treatments were administered continuously for 28 days. After anesthesia and euthanasia, serum and colonic content samples were collected from six mice per group. These samples were subjected to metabolomic profiling and 16S rRNA gene sequencing of the gut microbiota.

2.3 Retinal morphological observation

Following fixation, retinal tissues were dehydrated, cleared, and embedded in paraffin. Sections were then sequentially stained with hematoxylin to visualize nuclei, differentiated in hydrochloric acid alcohol, blued, counterstained with eosin for cytoplasmic visualization, dehydrated through a graded ethanol series, cleared in xylene, and finally mounted with neutral balsam. Stained retinal sections from each group were scanned using a digital slide scanner, and pathological alterations were examined and documented.

2.4 Prediction of target genes for DHYZ

DHYZ is a formula comprising 15 herbs (Table 1). To identify its active compounds, the TCMSP database (https://tcmsp-e.com/tcmsp.php) was first searched by entering each of the 12 herb names individually in the “Herbname” field, excluding Dendrobium officinale Kimura & Migo (Shi Hu), Ophiopogon japonicus(L.f) Ker-Gawl. (Mai Dong), and Polygala tenuifolia Willd. (Yuan Zhi). The active compounds of these 12 herbs were then mapped to their target genes using the UniProt database (19) (https://www.uniprot.org/). After merging and deduplication, these targets were designated as target gene 1.

Botanical plant nameChinese nameRehmannia glutinosa Libosch.Di HuangMorinda officinalis HowBa Ji TianCornus officinalis Siebold.et ZuccShan Zhu YuDendrobium nobile Lindl.Shi HuCistanche deserticola Y.C.MaRong Cong RongAconitum carmichaelii Debx.Fu ZiSchisandra chinensis (Turcz.) Baill.Wu Wei ZiCinnamomum cassia PreslRou GuiPoria cocos (Schw.) WolfFu LingOphiopogon japonicus (L.f) Ker-Gawl.Mai DongAcorus calamus L.Chang PuPolygala tenuifolia Willd.Yuan ZhiMentha haplocalyx Briq.Bo HeZingiber officinale Rosc.Sheng JiangZiziphus jujuba Mill.Da Zao

Next, the active compounds and corresponding target genes of the remaining three herbs—Shi Hu, Mai Dong, and Yuan Zhi—were retrieved from the HERB database (20) (http://herb.ac.cn). These targets were merged and deduplicated to form target gene 2.

Finally, target set 1 and target gene 2 were combined and deduplicated to obtain the complete set of predicted target genes for DHYZ. To ensure data reliability, TCMSP and HERB databases were used in a complementary manner. TCMSP served as the primary source for the 12 herbs, while HERB supplemented the three herbs not fully covered by TCMSP. The ETCM database was not included as it lacks complete information for key herbs such as Rehmannia glutinosa. The combined use of TCMSP and HERB thus provided comprehensive coverage of active compounds and targets for all 15 herbs in DHYZ.

2.5 Identification and enrichment analysis of candidate genes

DR-associated target genes were retrieved from Genecards (https://www.genecards.org/), OMIM (https://www.omim.org/), and Therapeutic Target Database (21) (https://db.idrblab.net/ttd/) employing “diabetic retinopathy” as the query. Targets from OMIM and TTD were merged with the top 500 relevance-ranked genes from Genecards. Following consolidation and deduplication, we established the DR-associated target gene set. Candidate genes were determined by intersecting DHYZ-predicted targets with DR-associated genes via the VennDiagram package (v 1.7.3) (22).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted using clusterProfiler package (v 4.2.2) (23) to elucidate biological functions and pathways of candidate genes, applying P < 0.05 as the significance cutoff. Results underwent p-value-based ascending sorting. For each GO domain-biological process (BP), cellular component (CC), and molecular function (MF)-the five most significantly enriched terms were chosen, alongside the ten leading KEGG pathways, for subsequent analysis and graphical representation.

2.6 Identification of key target genes

Subsequently, candidate genes were imported into the STRING database (24) (http://string-db.org) for protein-protein interaction (PPI) network construction using a confidence threshold exceeding 0.90. Cytoscape software (v 3.9.1) (25) visualized interactions among the top 15 candidates prioritized by the Closeness algorithm. The cytoHubba plugin facilitated additional refinement by assessing genes across six topological parameters: Closeness, Degree, Edge Percolated Component (EPC), Maximal Clique Centrality (MCC), Maximum Neighborhood Component (MNC), and Radiality. From each ranking approach, the top 10 genes were extracted, with overlapping genes designated as key targets via the UpSetR package (v 1.4.0) (26). Finally, Cytoscape integrated and visualized key targets, their cognate bioactive compounds, and source herbs as a comprehensive network.

2.7 Molecular docking analysis

Molecular docking was employed to examine binding interactions between pivotal target proteins and their cognate bioactive compounds. Three-dimensional compound structures were obtained from PubChem database (https://pubchem.ncbi.nlm.nih.gov/), while protein architectures were sourced from UniProt database. CB-Dock web server (27) (https://cadd.labshare.cn/cb-dock/php/blinddock.php) executed docking simulations. Binding energy calculations quantified ligand-target interaction strength, wherein lower energies signify enhanced binding affinity, with values below −5 kcal/mol indicating favorable interactions. PyMOL software (v 2.5) rendered visualization of docking outcomes (28).

2.8 Metabolome sequencing and data pre-processing

The collected samples were thawed on ice, and metabolites were extracted using pre-chilled 50% methanol. After vortexing and incubation at room temperature, the samples were stored at −80 °C. The extracts were then centrifuged, and the supernatant was transferred to a new 96-well plate and stored at −80 °C until LC-MS analysis. To ensure data quality, quality control (QC) samples were included in each batch.

Chromatographic separation was performed using an ultra-high-performance liquid chromatography (UPLC) system equipped with a dedicated column under a specific gradient elution program. Mass spectrometry analysis was carried out on a high-resolution TripleTOF 5,600 plus mass spectrometer, which alternated between positive and negative ion modes. Data were acquired in information-dependent acquisition (IDA) mode with standard calibration. System stability was monitored by injecting QC samples throughout the acquisition process.

The raw data were converted to mzXML format using Proteowizard (v 3.0) (29). Peak detection, alignment, and integration were subsequently performed using XCMS (30). For initial metabolite identification, the accurate m/z values were queried against the Human Metabolome Database (HMDB) and KEGG database using MetaX software (v 2.0.0) (31). For further annotation, secondary spectra generated from ion fragmentation were matched against a curated spectral library.

2.9 Principal component analysis (PCA) and orthogonal partial least squares discriminant analysis (OPLS-DA)

Following this, PCA was employed using the stats (v 4.3.1) package to assess the quality of metabolome data among the DHYZ, model, and control groups. To illustrate the relationship between metabolite expression levels and sample categories, while maximizing sample separation and predicting sample categories, OPLS-DA was applied to the metabolome sequencing data for the model and control groups, as well as for the DHYZ and model groups, in both positive and negative ion modes using the ropls (v 1.38.0) package. The Variable Importance in Projection (VIP) values for each metabolite were then obtained. Additionally, multiple comparisons (999 permutations) were conducted, and the results were corrected using the false discovery rate (FDR) method to assess the robustness of each OPLS-DA model.

2.10 Differential expression analysis and identification of candidate metabolites

Differentially expressed metabolites (DEMs) were identified using the limma package (v 3.58.1) (32) with the criteria of |log2 Fold Change (FC)| > 1.0, P < 0.05, and VIP >1. The analysis was performed under two ion modes for specific group comparisons: In positive ion mode, comparisons between the model and control groups yielded DEMs1, and between the DHYZ and model groups yielded DEMs2. In negative ion mode, the corresponding comparisons produced DEMs3 (model vs. control) and DEMs4 (DHYZ vs. model). Volcano plots for DEMs1–DEMs4 were generated using the ggplot2 package (v 3.5.1), labeling the top 10 up- and down-regulated metabolites based on |log2 FC|.

Venn diagrams were constructed using the VennDiagram package (v 1.7.3) to identify metabolites with reversed expression trends. For the positive ion mode data, the intersection between up-regulated metabolites in DEMs1 and down-regulated metabolites in DEMs2 was defined as intersection metabolites 1. Conversely, the intersection between down-regulated metabolites in DEMs1 and up-regulated metabolites in DEMs2 was defined as intersection metabolites 2. The union of these two sets constituted candidate metabolites 1. A parallel analysis was performed on the negative ion mode data (DEMs3 and DEMs4) to identify intersection metabolites 3 and 4, the union of which formed candidate metabolites 2. The overall candidate metabolite set was defined as the union of candidate metabolites 1 and 2. The abundance profiles of these final candidate metabolites across all samples were visualized in a heatmap using the pheatmap package (v 1.0.12). To explore the KEGG pathways associated with the candidate metabolites, an enrichment analysis was performed using the MetaboAnalyst platform, with significance set at P < 0.05. The results were visualized as a bubble chart using ggplot2.

2.11 Microbiome sequencing and data pre-processing

Amplification targeted conserved 16S rDNA/ITS2 gene regions employing total genomic DNA as template material. Target sequences underwent 35-cycle amplification via one-step PCR utilizing Phusion High-Fidelity DNA Polymerase. Subsequently, a secondary PCR step facilitated attachment of indexed sequencing adapters for library construction. Amplicon libraries experienced purification through AxyPrep PCR Cleanup Kit, underwent quantification via Quant-iT PicoGreen dsDNA Assay Kit on the Promega QuantiFluor platform, and were combined in equimolar proportions. Paired-end sequencing (2 × 250 bp) was ultimately executed on an Illumina platform following established protocols. Raw sequence outputs underwent processing through the DADA2 pipeline, encompassing quality filtration, denoising, and paired-end read consolidation, succeeded by chimeric sequence elimination via Vsearch to yield amplicon sequence variant (ASV) tables. ASV taxonomic classification was accomplished by sequence alignment to the SILVA database (https://www.arb-silva.de/) employing the Mothur algorithm. Sequencing depth sufficiency was evaluated through rarefaction curve generation for all three-group samples using ggplot2 package (v 3.5.1), examining the correlation between observed species diversity and sequencing intensity.

2.12 Alpha and beta diversity analyses and species composition profiling

Alpha diversity, quantifying species richness, diversity, and evenness within localized environments, was determined by computing ACE, Chao1, Shannon, and invSimpson metrics via the vegan package (v 2.6.4). Inter-group statistical variations in these metrics across DHYZ, model, and control cohorts were assessed through Wilcoxon rank-sum testing, applying P < 0.05 significance criteria. Violin plots generated via ggplot2 (v 3.5.1) provided result visualization. Beta diversity, reflecting inter-community compositional variation, was examined using Bray-Curtis and Jaccard distance measures. Principal coordinates analysis (PCoA) facilitated dissimilarity visualization, while permutational multivariate analysis of variance (PERMANOVA) tested statistical significance of group clustering. Microbial community compositional profiling across cohorts involved calculating relative abundance of microbial taxa per sample. The ggplot2 package generated community composition bar charts displaying the top 10 most prevalent taxa at phylum and genus hierarchical levels. Additionally, genus-level taxonomic entities were compared among groups to discern group-specific or shared taxa, visualized through the Venn Diagram package (v 1.7.3).

2.13 Functional and linear discriminant analysis effect size (LEfSe) analyses

Microbial community functional variations among the three cohorts were forecasted based on KEGG pathway annotations using PICRUSt2 (https://github.com/picrust/picrust2). Inter-group disparities in predicted KEGG pathway relative abundances were evaluated via Kruskal-Wallis testing (P < 0.05). For identifying microbial taxa exhibiting significant differential prevalence across groups, Linear Discriminant Analysis Effect Size (LEfSe) methodology was implemented via the microeco package (v 1.8.0). This assessment employed significance criteria of P < 0.05 coupled with linear discriminant analysis (LDA) scores exceeding 2. Genera satisfying these statistical thresholds were designated as candidate microbes. LEfSe analytical outcomes were visualized through circular clustering diagrams produced with the circlize package (v 0.4.16) and bar charts created using ggplot2 (v 3.5.1).

2.14 Omics integration analysis

Spearman correlation analysis via the psych package (v 2.4.3) investigated associations between candidate metabolites and candidate microbes. Correlations exhibiting absolute coefficients (|ρ|) exceeding 0.30 with P < 0.05 were deemed statistically significant and displayed in correlation heatmaps generated through the ComplexHeatmap package (v 2.21.1).

To emphasize the most robust associations, more stringent criteria (|ρ| > 0.60, P < 0.05) were implemented to delineate subsets of key metabolites and key microbes. Sankey diagrams were subsequently constructed to depict interaction networks among these pivotal elements utilizing the ggsankey package (v 3.7.1).

Lastly, to discern biological pathways co-enriched by both key target genes (from antecedent analysis) and key metabolites, joint pathway examination was executed via the Joint-pathway module of the MetaboAnalyst platform (P < 0.05). Outcomes were rendered using the ggplot2 package (v 3.5.1).

2.15 Statistical analysis

R software (v 4.2.2) executed all statistical assessments. Inter-group comparisons employed either Wilcoxon rank-sum testing (two-group scenarios) or Kruskal-Wallis testing (multi-group scenarios). Statistical significance was established at P-values below 0.05.

3 Results3.1 Results of retinal histopathological staining

In the control group, the retinal tissue exhibited well-defined laminar architecture, compact cellular arrangement, and showed no signs of edema or vascular abnormalities. In contrast, the model group displayed marked retinal edema, a reduction in ganglion cells, increased spacing with sparse cellularity in the outer nuclear layer, disorganized cell arrangement, and the presence of neovascularization (Figure 1A).

Histopathological staining confirms the therapeutic efficacy of DHYZ on DR, and network pharmacology reveals the multi-target mechanisms of DHYZ against DR. (A) Representative hematoxylin and eosin (HE) staining images of retinal tissue (Magnification × 400; Scale bar = 50 μm). (B) Venn diagram showing 110 candidate genes from the intersection of 703 DHYZ-related targets and 699 DR-related targets. (C) Gene Ontology (GO) enrichment analysis. (D) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis.

3.2 Recognition of 703 target genes for DHYZ

The active compounds of the 12 herbs in DHYZ were identified through the TCMSP, with the following numbers of compounds: 76 for Shudihuang, 174 for Baji-tian, 226 for Shanzhuyu, 75 for Roucongrong, 65 for Fuzi, 169 for Wuweizi, 100 for Rougui, 34 for Fuling, 105 for Shichangpu, 164 for Bohe, 265 for Shengjiang, and 133 for Dazao. A total of 556 target genes 1 were then retrieved from the Uniprot database.

In addition, the HERB database provided data for 3 additional herbs: 58 active compounds and 158 target genes for Shi Hu, 55 active compounds and 3 target genes for Mai Dong, and 101 active compounds and 108 target genes for Yuan Zhi yielding a total of 235 target genes 2. The intersection of the 556 target genes 1 and the 235 target genes 2 resulted in 703 target genes for DHYZ.

3.3 Identification and enrichment analysis of 110 candidate genes

Following explorations in the OMIM and TTD databases, 234 and 30 target genes associated with DR were identified, respectively. The 234 target genes from the OMIM database, 30 target genes from the TTD database, and 500 genes from the Genecards database were merged and deduplicated, yielding 699 DR-related target genes. Subsequently, the intersection of 703 target genes for DHYZ and 699 DR-related target genes resulted in 110 candidate genes (Figure 1B).

Further analysis of these 110 candidate genes using GO analysis identified 2,584 terms (P < 0.05), including 2,331 BP terms, 77 CC terms, and 176 MF terms. At the GO level, these candidate genes were primarily associated with functions such as “response to decreased oxygen levels” (BP), “membrane microdomain” (CC), and “cytokine activity” (MF) (P < 0.05) (Figure 1C). Moreover, 160 KEGG pathways were identified, including the “IL-17 signaling pathway”, “HIF-1 signaling pathway”, and “AGE-RAGE signaling pathway in diabetic complications” (P < 0.05) (Figure 1D). These analyses provided a strong foundation for understanding the functional importance of candidate genes in the progression of DR.

3.4 Selection of 5 key target genes

The PPI network of the top 15 candidate genes revealed that STAT3, EGFR, and TP53 were hub nodes with extensive connections (Figure 2A). Five key target genes (STAT3, IL6, TNF, ESR1, and IL1B) were subsequently identified as the consensus core set by intersecting the top-ranked genes from six distinct topological algorithms (Figure 2B). To visualize the relationships between these key targets, their corresponding active compounds, and the source herbs, a comprehensive “compound-target-herb” network was constructed (Figure 2C).

Identification of key targets and construction of the compound-target-herb network. (A) Protein-protein interaction (PPI) network of the top 15 candidate genes. Each node represents a protein, and the size of the node usually reflects its degree value, while the thickness of the lines typically indicates the confidence level or strength of the interaction. The node labels directly indicate the official gene symbols of the proteins (such as IL10, ESR1, etc). (B) Identification of 5 key target genes based on the overlap from six topological analysis methods. This figure consists of two parts. The left part shows the number of genes in each set, while the right part indicates the number of genes that are common to the combinations of different algorithms. (C) Integrated network diagram incorporating key targets, their corresponding active compounds, and related herbal constituents. The green nodes represent traditional Chinese medicines, the blue nodes represent active ingredients, and the yellow nodes represent candidate target genes.

3.5 Molecular docking analysis of ESR1 with bioactive compounds

Molecular docking was performed to characterize the binding interactions between the protein ESR1 and six corresponding active compounds (salidroside, schisanhenol, naringenin, 6-gingerol, hyperin, and leonuride). All six compounds exhibited strong binding affinities to ESR1, with naringenin showing the highest affinity (Table 2, Figures 3A–F). The detailed binding parameters are as follows: Naringenin displayed the strongest binding energy (-8.9 kcal/mol), forming hydrogen bonds with residues ARG-394 and GLU-353 (Figure 3C). Salidroside bound at −8.1 kcal/mol, interacting via hydrogen bonds with ARG-394 and GLU-353 (Figure 3A). Hyperin and leonuride showed binding energies of −7.9 kcal/mol and −7.8 kcal/mol, interacting with residues ARG-394/GLY-521 and ARG-394/LEU-346/LEU-387, respectively (Figures 3E, F). Schisanhenol and 6-gingerol had binding energies of −6.7 kcal/mol and −5.6 kcal/mol, forming hydrogen bonds with ARG-412 and LEU-462/SER-468, respectively (Figures 3B, D). These results confirm robust binding between ESR1 and the bioactive compounds, supporting their potential functional relevance.

GeneActive compoundsBinding energyESR1salidroside−8.1 kcal/molESR1Schisanhenol−6.7 kcal/molESR1naringenin−8.9 kcal/molESR16-gingerol−5.6 kcal/molESR1hyperin−7.9 kcal/molESR1leonuride−7.8 kcal/mol

Binding energies of ESR1 and 6 active compounds.

Molecular docking analysis reveals interactions between ESR1 and active compounds. (A–F) Binding modes of ESR1 with salidroside (A), schisanhenol (B), naringenin (C), 6-gingerol (D), hyperin (E), and leonuride (F), respectively. (The protein is depicted as a yellow cartoon structure, while the compound is represented by red stick models).

3.6 Quality assessment of metabolomic profiles

PCA showed clear separations among the control, model, and DHYZ groups, indicating distinct metabolic profiles (Figure 4A). OPLS-DA further demonstrated clear metabolic separations between the model and control groups, and between the DHYZ and model groups, in both positive and negative ion modes (Figures 4B, C). The robustness of the OPLS-DA models was validated by permutation tests (999 iterations), which confirmed the models were statistically significant and not overfitted (Figures 4D, E). The distinct group separations and validated models demonstrate the high quality and reliability of the metabolomic data, providing a solid foundation for subsequent analyses.

Metabolomic data quality assessment demonstrates model stability and reliability. (A) PCA plot showing clear separation among control, model, and DHYZ groups. (B, C) Orthogonal projections to latent structures-discriminant analysis (OPLS-DA) score plots. (D, E) Permutation test results validating the OPLS-DA models.

3.7 Identification and functional exploration of candidate metabolites

Differential expression analysis (thresholds: |log2FC| > 1.0, P < 0.05, VIP > 1) identified 82 DEMs1 [25 up-, 57 down-regulated, Differential metabolite 1 (MC1)], 76 DEMs2 [46 up-, 30 down-regulated, Differential metabolite 2 (DM1)], 94 DEMs3 [26 up-, 68 down-regulated, Differential metabolite 3 (NMC)], and 63 DEMs4 [48 up-, 15 down-regulated, Differential metabolite 4 (NDM)] in the respective group comparisons (Figure 5A). Venn analysis was used to identify metabolites with reversed expression trends. In the positive ion mode, the intersection of DEMs1 and DEMs2 yielded 8 intersection metabolites 1 and 22 intersection metabolites 2, which were combined into 30 candidate metabolites 1 (Figure 5B). Similarly, analysis of DEMs3 and DEMs4 in the negative ion mode identified 3 intersection metabolites 3 and 17 intersection metabolites 4, together forming 20 candidate metabolites 2 (Figure 5C). The union of these two sets resulted in a final list of 50 candidate metabolites, whose expression patterns across the three groups are displayed in a heatmap (Figure 5D).

Metabolite identification and functional analysis. (A) Volcano plot illustrating up- and down-regulated differential metabolites. Each point in the figure represents a metabolite. The red points represent up-regulated metabolites, the blue points represent down-regulated metabolites, and the gray points represent metabolites with no significant difference. (B) The intersection situation of up-regulation and down-regulation for MC1 and DM1. MCU and MCD, respectively represent the up-regulated genes and down-regulated genes of MC1; DMU and DMD, respectively represent the up-regulated and down-regulated genes of DM1. (C) The intersection situation of up-regulation and down-regulation of NMC and NDM. NMCU and NMCD, respectively represent the up-regulated genes and down-regulated genes of NMC; NDMU and NDMD, respectively represent the up-regulated genes and down-regulated genes of NDM. (D) Expression heat map of key metabolites. The horizontal axis represents the samples, and the vertical axis represents the key metabolites. The darker the color, the higher the expression level of the metabolite in that sample. (E) KEGG pathway enrichment analysis bubble chart for key metabolites. The vertical axis represents the description of each pathway, and the larger the circle, the more key metabolites are enriched in that pathway.

Enrichment analysis of the 50 candidate metabolites revealed four significant KEGG pathways (P < 0.05): arachidonic acid metabolism, alpha-Linolenic acid metabolism, glycerophosp