Active constituents and targets of Astragalus membranaceus (Fisch.), Angelica sinensis (Oliv.) Diels, and Panax notoginseng (Burk.) F.H. Chen were identified via the Traditional Chinese Medicine Systems Pharmacology (TCMSP) platform (https://tcmsp-e.com/tcmsp.php). Criteria for filtering included oral bioavailability (OB) of at least 30% and drug-likeness (DL) of at least 0.18. Target validation was subsequently performed utilizing the UniProt database.

Target genes associated with diabetic retinopathy were retrieved by searching Genecards (https://www.genecards.org/), OMIM (https://omim.org/), PharmGKB (https://www.pharmgkb.org/), and Drugbank (https://go.drugbank.com/). The collected target genes were consolidated, and duplicates were eliminated.

The Venn package was employed to identify intersections between active ingredients in MDBD and DR-related targets. Cytoscape V3.10.0 was then used to construct a network, incorporating diseases, drugs, components, and DR-related targets as nodes, with their interrelationships represented as edges. Topological analysis was conducted to identify core components within the network.

The intersecting targets’ PPI network was retrieved from the String database (https://string-db.org), specifying Homo sapiens as the species and setting the minimum interaction score at 0.40. Using Cytoscape V3.10.0, a network was constructed with nodes representing diseases, drugs, components, and DR-related targets, while edges depicted their relationships. To pinpoint key components and core PPI targets, the MCODE plugin was applied for clustering and topological analysis.

The keyword ‘Diabetic retinopathy’ guided sample searches in the GEO database (https://www.ncbi.nlm.nih.gov/geo/), with constraints on data type (expression profiles) and species (Homo sapiens). Gene expression and clinical data were extracted to evaluate core gene expression levels related to MDBD treatment for diabetic retinopathy across normal and DR groups.

Core gene expression levels were extracted from control and DR groups, followed by differential expression analysis using R packages such as ‘limma’, ‘pheatmap’, and ‘ggpubr’ in R V4.3.3. Results were visualized via box plots and heatmaps, identifying genes with P < 0.05 as SDECGs. Perl scripts determined the chromosomal positions of these genes, which were then displayed in circular plots using the R package ‘Rcircos’. The ‘cor’ function calculated and visualized correlation coefficients for each SDECG.

The CIBERSORT command in R was employed to conduct 1000 simulations, determining the relative abundance of immune cells. A bar graph was used to illustrate immune cell content per sample. Single-sample gene set enrichment analysis (ssGSEA) was performed using R packages ‘GSVA’ and ‘GSABase’ to compare immune cell content differences between normal and DR groups, with results shown in box plots. Correlation testing was conducted between SDECGs and ssGSEA scores, and coefficients were visualized. CIBERSORT results were finalized.

Four predictive models, Random Forest (RF), Support Vector Machine (SVM), Generalized Linear Model (GLM), and Extreme Gradient Boosting (XGB), were developed based on SDECG expression data. Feature genes were selected using residual cumulative distribution plots, residual boxplots, and Receiver Operating Characteristic (ROC) curves. The optimal model was constructed incorporating feature genes and their expression levels across normal and DR groups. Decision and calibration curves were formulated to evaluate model accuracy.

A two-sample Mendelian Randomization (MR) analysis was conducted to explore the causal link between feature genes and DR risk, with SNPs designated as instrumental variables (IVs). Feature gene SNPs were sourced from the Integrative Epidemiology Unit (IEU) database (https://gwas.mrcieu.ac.uk/) as exposure variables, while DR-related SNPs were outcome variables. The “TwoSampleMR” package facilitated MR analysis, employing the inverse variance weighted (IVW) method to assess the correlation between feature gene expression levels and DR risk. Heterogeneity was examined using Cochran’s Q statistic, with P < 0.05 indicating heterogeneity in IVW results. MR-Egger regression and MR-PRESSO analysis were used to evaluate potential horizontal pleiotropy, with P < 0.05 suggesting its presence.

The 3D structures of the core components and feature genes of MDBD were retrieved from the PubChem (https://pubchem.ncbi.nlm.nih.gov) and Protein Data Bank (PDB) (http://www.rcsb.org/) databases. Molecular docking was performed using Autodock Vina to preliminarily verify the interaction between the core network pharmacological components and feature genes, and the top four docking combinations were selected and visualized with Pymol. Differential expression analysis of GEO datasets was performed to identify key genes between the normal and DR groups, providing additional validation for the molecular mechanisms.

The R package ‘ConsensusClusterPlus’ was utilized to categorize DR samples based on SDECG expression, employing the k-means clustering technique with Euclidean distance, allowing for up to nine clusters. Comparative analysis of clustering outcomes was performed using heatmaps and boxplots to evaluate expression levels. Principal Component Analysis (PCA) was conducted to discern inter-cluster differences. Subsequent ssGSEA analysis on SDECG clusters generated bar plots to depict immune cell content variations across clusters. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were executed using GMT files from the GSEA platform (http://www.gsea-msigdb.org/). Gene Set Variation Analysis (GSVA) was conducted in R V4.3.3 to assess enriched gene expression across clusters. Differential gene expression analysis, adhering to|logFC|>1 and adj.P-Value<0.05, was followed by Venn diagram visualization to identify differentially expressed genes (DEGs) among clusters.

The DEGs among SDECG clusters underwent GO enrichment analysis covering biological processes (BP), molecular functions (MF), and cellular components (CC), alongside KEGG pathway enrichment. These analyses were executed using R packages such as “clusterProfiler” and “enrichplot” in R V4.3.3, with a screening threshold of P-value <0.05. Results were depicted through circular plots and bar graphs.

An additional cluster analysis, based on DEG expression, was performed using the methodology described in Section 2.11, selecting the DEG cluster with the highest precision. The DEG clustering results were used to compare DEG expression levels across clusters, SDECG expression variations, and immune cell content differences among clusters. These findings were visualized using heatmaps and box plots.

Scores for SDECGs were computed for each sample based on expression levels, utilizing PCA by incorporating PC1 and PC2. Differential analysis of SDECGs and DEGs clustering scores was executed using R packages such as “limma” and “ggpubr” in R V4.3.3. Box plots were generated to depict SDECGs scores across clustered samples. Additionally, the “ggalluvial” package facilitated the creation of alluvial diagrams to illustrate relationships and processes among SDECGs clusters, DEGs clusters, and samples with varying SDECGs scores.

Beta-sitosterol (HY-N0171A, MedChemExpress, USA) was dissolved in DMSO (D8371, Solarbio, USA). The solution was subsequently diluted in complete culture medium for mouse retinal Müller cells (CM-M117, Procell, China) to prepare various concentrations as required for the experiments. Fresh solutions were prepared immediately prior to use.

The mouse retinal Müller cells (MMCs) (CP-M117, Procell, China) were purchased from Pricella Biotechnology Co.,Ltd. The cells were cultured in in complete culture medium for mouse retinal Müller cells (CM-M117, Procell, China), and maintained at 37.0 °C with a CO₂ concentration of 5% in a cell culture incubator. When the confluence of MMCs reached over 80%, they were passaged using trypsin (2.5% EDTA) (SH30042.01, HyClone, USA).

We induced DR cell model according to our previous study (35), utilizing a hypoxia-induced injury approach to simulate the conditions encountered in DR. MMCs were cultured under hypoxic conditions by placing them in a hypoxia chamber with a controlled oxygen concentration of 1% for a duration of 48 hours. Beta-sitosterol were used to incubated MMCs for 24h. The cells were harvested for subsequent analysis.

The RiboFECT™ CP transfection kit (C10511-05, RiboPharm, China) facilitated cell transfection in a 6-well plate. A solution of siERBB2 was prepared by mixing 3μl of 30μM stock with 120μl of 1X riboFECT™ CP Buffer, followed by the addition of 12μl of riboFECT™ CP reagent to form the transfection complex. This complex was combined with an antibody-free complete culture medium and gently agitated. Subsequent treatments led to incubation of the plates at 37 °C in a CO₂ incubator for 24 hours prior to assays.

Cell viability assessment was conducted using the Cell Counting Kit-8 (CCK-8) (K009-1000T, Zeta-life, USA). MMCs were exposed to various concentrations of Beta-sitosterol in medium at 37 °C for 24 hours following hypoxia induction. Post medium removal, 10 μL of CCK-8 solution was introduced to each well and incubated in a CO₂ incubator for 1 hour. Absorbance was recorded at 450 nm utilizing a microplate reader.

Total RNA extraction from MMCs was performed using TRIzol (DP424, Tiangen, China). cDNA synthesis followed with the RevertAid First Strand cDNA Synthesis Kit (K1622, Thermo, USA). RT-qPCR was employed to evaluate mRNA expression using 2x SYBR Green qPCR Master Mix (B21203, Selleck, USA) on a CFX96™ real-time system (Bio-Rad, USA). Relative expression levels of target genes were determined via the 2^-ΔΔCt method, based on obtained Ct values. Primers for analyzed genes included:

Erb-B2 Receptor Tyrosine Kinase 2 (ERBB2): forward, 5’-CTGTGTGACCACCTGCCCCTAC-3’, and reverse, 5’-TGCCCAGACCATAGCATACTCC-3’; Beta-actin (β-actin): forward, 5’-AGGTCATCACTATTGGCAACGAG-3’, and reverse, 5’-TTGGCATAGAGGTCTTTACGGAT-3’.

MMC lysis was achieved using the RIPA buffer (R0020, Solarbio, USA) to facilitate protein extraction for subsequent analyses. Quantification of proteins in the lysis supernatant was conducted employing the BCA Protein Assay Kit (P0011, Beyotime, China). Proteins were separated via SDS polyacrylamide gel electrophoresis (SDS-PAGE) utilizing the SDS-PAGE Gel Kit (P1200, Solarbio, USA), followed by transfer onto polyvinylidene fluoride (PVDF) membranes (IPVH00010, Millipore, Germany). After blocking and washing, the membranes underwent incubation with a primary antibody for 20 hours, succeeded by a 2-hour incubation with a secondary antibody. The proteins on PVDF membranes were processed with BeyoECL Plus (P0018S, Beyotime, China) to enhance chemiluminescence. Antibodies employed for Western blotting included ERBB2 (18299-1-AP, Proteintech), and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) (UM4002, UtiBody, China).

The quantification of interleukin-1β (IL-1β), vascular endothelial growth factor (VEGF), and angiopoietin-like 6 (ANGPTL6) in Müller cells (MMCs) was performed using specific enzyme-linked immunosorbent assay (ELISA) kits (MM-0040M1, MM-0128M1, MM-48029M1, Meimian, China). ELISA assays were conducted strictly following the manufacturer’s protocols to ensure accuracy and reproducibility. Absorbance measurements were taken at 450 nm using a calibrated microplate reader, and the concentrations were calculated by comparing the absorbance values to standard curves generated during the assay.

The statistical analysis was conducted using GraphPad Prism 10 software (GraphPad Software, USA). Quantitative data were expressed as mean ± standard error of the mean (SEM). For comparisons among multiple groups, one-way analysis of variance (ANOVA) was applied, followed by Fisher’s LSD test or Tukey’s post hoc test for pairwise comparisons. Differences between two groups were analyzed using the Student’s t-test. Statistical significance was determined by a Type I error probability (α) of less than 5%, with significance levels presented as P < 0.05, P < 0.01, or P < 0.001 to denote varying degrees of statistical difference.

Active constituents and targets of Angelica sinensis, Astragalus membranaceus, and Panax notoginseng were acquired via the TCMSP database. Following the elimination of duplicates and irrelevant entries, 30 distinct active components and 1207 associated targets were identified ( Supplementary Table S1 ).

Targets numbering 4848, 185, 64, and 4 were sourced from Genecards, Pharmgkb, OMIM, and Drugbank, respectively. After duplicate removal, a total of 4973 DR-related targets were obtained ( Figure 2 ).

The intersection between each active ingredient in MDBD and DR-related targets was determined using the venn package, resulting in 130 targets directly linked to drug and disease ( Figure 3A ). Topological analysis was conducted to pinpoint core components within the network ( Figure 3B ).

The PPI network of intersecting targets was extracted from the String database, with Homo sapiens specified as the species. The network comprised 80 nodes and 1318 edges, with an average node degree of 35.5. The MCODE plugin facilitated clustering of the PPI network to identify core PPI targets, yielding two clustered networks with a minimum interaction score of 0.40. Seventy-four core targets were identified, including MAPK8, CASP9, IL1A, IFN, EGFR, IL1B, CDKN1A, SOD1, ERBB2, NOS3, and CD40LG, which are deemed significant for MDBD treatment of DR ( Supplementary Table S2 ; Figure 4 ).

The keyword “diabetic retinopathy” was utilized to extract samples from the GEO database, with restrictions on data type and biological species. Gene expression and clinical data were obtained to assess the expression levels of core genes in MDBD-treated DR within normal and DR cohorts. Two datasets were selected: GSE160306, GPL20301. Macular tissue analysis was performed across various clinical stages: DR for diabetic retinopathy; NPDR for non-proliferative diabetic retinopathy; PDR for proliferative diabetic retinopathy; DME for diabetic macular edema, where elevated scores denote more severe lesions. The study included 20 normal samples and 39 DR samples.

Seventy-four core genes of MDBD were identified through network pharmacology analysis ( Supplementary Table S2 ). Expression levels were extracted from both the control and DR groups, followed by differential expression analysis. Fourteen genes, including NR3C2, PSMG1, INSR, HIF1A, RXRB, IKBKB, CAV1, RXRA, TP53, CCND1, ERBB2, ADRB1, PLAT, and SERPINE1, were validated as SDECGs core genes in human samples. Except for HIF1A, PSMG1, ADRB1, and PLAT, which exhibited high expression in the control group, the remaining genes were highly expressed in the DR group ( Figures 5A, B ). Chromosomal positions of MDBD core genes were mapped in RCircos ( Figures 5C, D ). Correlation analysis among SDECGs in DR samples revealed strong correlations, encompassing both positive and negative relationships ( Figures 5E, F ).

Single-sample gene set enrichment analysis (ssGSEA) was conducted to compare immune cell levels between control and DR groups ( Figure 6A ), identifying statistically significant immune cells such as macrophages M0 with elevated expression in the DR group and resting dendritic cells with high expression in the control group. Immune cell infiltration analysis determined the types and levels of immune cells in each sample, visualizing these levels ( Figure 6B ). The intersection of SDECGs with ssGSEA scores for immune cell correlation testing indicated that correlations between SDECGs and immune cells included both positive and negative relationships ( Figures 6C, D ). Among the immune cells significantly associated with SDECGs (P < 0.05), dendritic cells resting, eosinophils, macrophages M0, macrophages M1, neutrophils, activated NK cells, resting CD4 memory T cells, CD8 T cells, and regulatory T cells (Tregs) were predominantly positively correlated with the relevant SDECGs. Conversely, memory B cells, naive B cells, macrophages M2, resting mast cells, monocytes, plasma cells, and follicular helper T cells were mainly negatively correlated with the relevant SDECGs.

Expression data of SDECGs served to construct four predictive models: RF, SVM, GLM, and XGB. Analysis of residual box plots, ROC curves, and reverse cumulative distribution plots demonstrated that the RF model achieved the highest accuracy, featuring the largest area under the ROC curve and minimal residuals and reverse cumulative values ( Figures 7A–C ). The RF model, selected for further development, facilitated the construction of a predictive model to derive feature importance scores for the genes. Results indicated that SERPINE1 had the highest score, with gene importance ranked as SERPINE1, TP53, ERBB2, CCND1, INSR, PLAT, NR3C2, CAV1, PSMG1, and RXRA ( Figure 7D ). The top five genes, CCND1, ERBB2, INSR, TP53, and SERPINE1, were utilized for nomogram construction ( Figure 7E ). Finally, decision and calibration curves assessed the predictive model’s accuracy, demonstrating high accuracy based on positive outcomes ( Figures 7F, G ).

MR analysis was executed to investigate the causal relationship between specific feature genes and DR, with SNPs defined as IVs. SNPs from signature genes were utilized as exposure factors, while SNPs related to DR served as outcome factors. Information on the SNPs for the four highlighted genes CCND1, ERBB2, INSR, and SERPINE1 can be found in Supplementary Table S3 . No SNPs were identified as weak IVs. Due to insufficient data, MR analysis could not be performed on TP53. IVW analysis indicated that ERBB2 is linked to an increased risk of DR (Odds Ratio [OR]=1.70, 95% CI: 1.016 to 2.856, P = 0.04), whereas the other three genes did not exhibit a significant causal relationship with DR ( Figure 8A ; Supplementary Table S4 ). The IVW results confirmed ERBB2 as a risk factor for DR, with consistent effect directions observed in the weighted median analysis (OR = 1.860, 95% CI: 1.247-2.774, P = 0.002), as demonstrated by both scatter and forest plots ( Figures 8B, C ).

Cochran’s Q test was employed to evaluate study heterogeneity, revealing significant heterogeneity among the selected instrumental variables (P = 0.003) ( Figure 8D ; Supplementary Table S5 ). Consequently, IVW with random effects was applied in all MR analyses to mitigate heterogeneity impact. Despite Cochran’s Q test detecting some heterogeneity, MR-Egger analysis showed minimal influence from horizontal pleiotropy (P>0.05) ( Figure 8D ; Supplementary Table S4 ), indicating no directional pleiotropy and affirming stable causal relationships. The leave-one-out analysis demonstrated consistent MR analysis results upon sequential SNP removal, indicating robustness and negligible impact on the overall findings ( Figure 8E ).

Autodock Vina was utilized for molecular docking to preliminarily validate interactions between core pharmacological components and feature genes. The optimal docking combinations were selected and visualized using Pymol. Through molecular docking analysis of feature genes and core components of MDBD, it was observed that the binding energy for all docking combinations was below -6.0 kcal/mol, suggesting stable structures could be formed between feature genes and core components ( Figure 9A ; Supplementary Table S6 ).

Figures 9B–E showed the top docking combinations of four key gene proteins with the core components of MDBD, while Figure 9F presented differential expression analysis of GEO datasets between the control and DR groups. Beta-sitosterol exhibited high binding affinity with CCND1, ERBB2, INSR, and SERPINE1. Alkyl interactions with CCND1 residues CYS132, MET31, LEU56, and ALA153, along with van der Waals forces with CCND1 residues GLU162, ASP159, MET155, LYS114, MET113, HIS153, VAL27, ASN131, and ARG59 were noted, with a docking energy of -8.3 kcal/mol ( Figure 9B ). Conventional hydrogen bonds were formed with ERBB2 residues TYR100, GLY165, THR164, and SER162, along with van der Waals interactions with ERBB2 residues SER188, ASN99, TYR236, ASN235, and TYR166, and a pi-alkyl interaction with TYR100, resulting in a docking energy of -8.5 kcal/mol ( Figure 9C ). Beta-sitosterol also established a conventional hydrogen bond with INSR residues ARG1092 and a carbon hydrogen bond with GLU1207. Additional van der Waals interactions with residues GLU1094, PRO1093, SER1204, THR1203, PRO1104, MET1109, PRO1235, CYS1234, ASN1233, ASP1232, and PRO1231, as well as alkyl interactions with LUE1205 and PRO1103, were observed, with a docking energy of -8.4 kcal/mol ( Figure 9D ). Furthermore, beta-sitosterol demonstrated high binding affinity for SERPINE1 residues SER233 through conventional hydrogen bonding and ASP216 through carbon hydrogen bonding. Van der Waals interactions with SERPINE1 residues ARG265, LEU267, LEU231, ARG181, and GLN273, along with pi-alkyl interactions with residues ILE218, TYR204, MET195, VAL268, PHE352, and PRO237, were also noted, with a docking energy of -9.3 kcal/mol ( Figure 9E ).

K-means clustering, utilizing Euclidean distance and a maximum of nine clusters, was employed to categorize DR samples based on SDECG expression. The analysis revealed two distinct clusters, achieving optimal accuracy ( Figures 10A, B ) and demonstrating high cluster stability ( Figures 10C, D ). Subsequent examination of SDECG expression in these clusters identified significantly elevated levels of NR3C2, PSMG1, INSR, HIF1A, and RXRB in Cluster C1 (P < 0.05, P < 0.01, P < 0.001), whereas IKBKB, CAV1, RXRA, TP53, CCND1, ERBB2, ADRB1, PLAT, and SERPINE1 showed increased expression in Cluster C2 (P < 0.05, P < 0.01, P < 0.001) ( Figures 10E, F ). Principal component analysis (PCA) results indicated a discernible separation between Clusters C1 and C2, with Cluster C1 exhibiting a higher density ( Figure 10G ). To assess immune cell level variations between C1 and C2, ssGSEA was conducted, identifying statistically significant immune cell populations. Notably, Cluster C1 exhibited higher fractions of CD8+ T cells, resting dendritic cells, and resting mast cells (P < 0.05, P < 0.01, P < 0.001), while Cluster C2 had elevated levels of resting CD4+ memory T cells, regulatory T cells (Tregs), M0 macrophages, and eosinophils (P < 0.05, P < 0.01, P < 0.001) ( Figure 10H ). Immune cell infiltration analysis further characterized the immune cell types and levels in each cluster, visually depicted in Figure 10I .

GSVA revealed that terms related to the dicarboxylic acid biosynthetic process, pyrimidine ribonucleoside monophosphate biosynthetic process, neuron cellular homeostasis, nucleoside kinase activity, fatty acyl-CoA metabolic process, and oxoglutarate dehydrogenase complex were upregulated in C2 compared to C1. Conversely, terms associated with the secretory granule membrane, cyclin-dependent protein serine/threonine kinase inhibitor activity, lysine acetylated histone binding, polysaccharide binding, and positive regulation of histone H3 K4 methylation were downregulated ( Figure 9J ). KEGG analysis indicated that pathways such as proximal tubule bicarbonate reclamation, terpenoid backbone biosynthesis, riboflavin metabolism, pyruvate metabolism, oocyte meiosis, and xenobiotic metabolism by cytochrome P450 were upregulated in C2 relative to C1, whereas sulfur metabolism, cell cycle, circadian rhythm, ECM receptor interaction, and prion diseases were downregulated ( Figure 10K ).

Enrichment analysis was conducted to elucidate the biological functions and pathways linked to the 726 DEGs, thereby revealing the mechanisms through which MDBD influences DR from various perspectives. The DEGs underwent GO enrichment analysis for BP, CC and MF. Findings indicated significant enrichment of DEGs in BP, primarily involving neuron and cell development, encompassing neuron projection regulation, axon development, positive regulation of cell development, axonogenesis, and wound healing. Enrichment in CC included collagen-containing extracellular matrix, neuronal cell body, and endoplasmic reticulum lumen, while MF enrichment comprised growth factor binding, cytokine binding, and voltage-gated monoatomic ion channel activity ( Figures 11A, B ). KEGG enrichment analysis suggested DEG involvement in pathways such as complement and coagulation cascades, pertussis, TGF-beta signaling, focal adhesion, and the PI3K-Akt signaling pathway ( Figure 11C ).

K-means clustering, utilizing Euclidean distance with a maximum of nine clusters, was applied to classify DR samples based on the expression of 726 DEGs. This approach resulted in the division of samples into two distinct clusters, achieving precise outcomes ( Figures 12A, B ) and demonstrating high stability ( Figures 12C, D ). A subsequent analysis of DEG expression within these clusters identified 305 DEGs with significantly increased expression levels in Cluster CI and decreased levels in Cluster CII, including CD163L1, PCYT2, SLC45A1, and ELOVL6. Conversely, 421 DEGs showed elevated expression in Cluster CII and reduced levels in Cluster CI, such as LPAR6, CGNL1, SHISA9, and TNFRSF1B ( Figure 12E ). Principal Component Analysis (PCA) revealed a clear distinction between the clusters, with Cluster CI exhibiting a higher density than Cluster CII ( Figure 12F ). ssGSEA results indicated a high expression of naive B cells, CD8+ T cells, resting dendritic cells, and resting mast cells in Cluster CI (P < 0.01), whereas regulatory T cells (Tregs), activated NK cells, and M0 macrophages were more prevalent in Cluster CII (P < 0.05, P < 0.01, P < 0.001) ( Figure 12G ). Immune cell infiltration analysis characterized the types and levels of immune cells expressed in each sample from both clusters, visually represented in Figure 12H .

PCA was utilized to evaluate the SDECGs of MDBD in DR treatment, leading to the establishment of a scoring model. Differences in scores among various clusters were assessed to determine significant disparities in SDECGs and to establish the corresponding cluster relationships, ensuring the findings’ reliability. PCA results indicated a statistically significant difference between two SDECG clusters (P < 0.001) ( Figure 13A ) and two DEG clusters (P < 0.001) ( Figure 13B ), demonstrating that C1 has a higher score than C2, and CI surpasses CII. The alluvial plot suggested that C1 in SDECG clustering aligns with CI in DEG clustering, while C2 corresponds predominantly to CII ( Figure 13C ).

A cell viability assay demonstrated a significant reduction in viability in the hypoxia-induced model group compared to the control group (P < 0.001). However, treatment with beta-sitosterol at concentrations of 8μM, 12μM, 16μM, and 20μM significantly restored cell viability relative to the model group (P < 0.001 for all comparisons) ( Figure 14A ). Notably, the 12μM and 16μM beta-sitosterol treatment groups exhibited the most pronounced effects, warranting further investigation.

RT-qPCR analysis revealed a marked upregulation of ERBB2 mRNA expression in the model group compared to the control (P < 0.001). Treatment with beta-sitosterol at 12μM and 16μM significantly reversed this upregulation, resulting in a substantial downregulation of ERBB2 mRNA expression (P < 0.001 for all comparisons) ( Figure 14B ). Similarly, Western blot analysis showed a significant increase in ERBB2 protein expression in the model group relative to the control, while beta-sitosterol treatment at 12μM and 16μM significantly reduced ERBB2 protein levels (P < 0.001 for all comparisons) ( Figures 14C, D ).

ELISA analysis revealed that the hypoxia-induced model group showed a statistically significant elevation in IL-1β, VEGF, and ANGPTL6 levels compared to the control group (P < 0.001) ( Figures 14E–G ). Beta-sitosterol treatment at concentrations of 12μM and 16μM effectively mitigated these increases, with the suppressive effects being more pronounced at the higher concentration of 16μM (P < 0.001 for all comparisons) ( Figures 14E–G ).

Cell viability assays indicated a notable reduction in viability within the model group relative to the control group (P < 0.001). Conversely, siERBB2 treatment significantly enhanced cell viability compared to the model group (P < 0.001) ( Figure 15A ).

RT-qPCR analysis demonstrated that ERBB2 mRNA expression was substantially upregulated in the model group compared to the control (P < 0.001). This upregulation was effectively reversed in the siERBB2 treatment group, where ERBB2 mRNA expression was significantly downregulated compared to the model group (P < 0.001) ( Figure 15B ). Western blot analysis similarly confirmed that ERBB2 protein levels were markedly elevated in the model group versus the control group (P < 0.001). However, siERBB2 treatment significantly reduced ERBB2 protein expression compared to the model group (P < 0.001) ( Figures 15C, D ).

The ELISA results demonstrated a significant increase in IL-1β, VEGF, and ANGPTL6 levels in the hypoxia-induced model group compared to the control group (P < 0.001). Furthermore, ERBB2 inhibition notably reduced the elevated levels of IL-1β, VEGF, and ANGPTL6 in the model group, indicating its suppressive effect (P < 0.001 for both) ( Figures 15E–G ).

DM, classified within the “Xiao Ke” category in TCM, is characterized by symptoms such as excessive thirst, hunger, urination, and weight loss. Its earliest documentation appears in the Yellow Emperor’s Inner Canon (36). In TCM, DR is referred to as “Xiao Ke Mu Bing,” an ocular complication of “Xiao Ke,” driven by internal heat resulting from yin deficiency. This imbalance disrupts qi and blood flow, leading to stasis and contributing to DR progression (37, 38). MDBD addresses these pathological mechanisms by enhancing qi, nourishing blood, and regulating meridians. Clinical studies have shown that combining MDBD with calcium hydroxybenzenesulfonate reduces inflammatory and angiogenic markers such as high-sensitivity C-reactive protein, tumor necrosis factor-alpha, VEGF, and endothelin-1, effectively ameliorating DR (39, 40).

Through analysis of the “drug-component-target” network, the primary active compounds responsible for MDBD’s therapeutic effects on DR have been identified. Beta-sitosterol, a naturally occurring bioactive phytosterol can be found in Angelica sinensis and Panax notoginseng, is a key component of MDBD. Its chemical structure closely resembles cholesterol from mammalian cells (41), exhibiting antioxidant properties and anti-diabetic activity (42, 43). Beta-sitosterol exhibited high binding affinity with genes implicated in DR pathogenesis. In vitro experiments using a hypoxia-induced DR cell model demonstrated that beta-sitosterol significantly reduced IL-1β, VEGF, and ANGPTL6 secretion while improving cell viability, supporting its role in mitigating key pathological processes in DR. Moreover, a research on ADMET (absorption, distribution, metabolism, excretion, and toxicity) profiling on beta-sitosterol (44) verified that beta-sitosterol complies with Lipinski’s Rule of Five, showing high intestinal absorption and effective permeability across the blood-brain barrier. These pharmacokinetic attributes indicate that beta-sitosterol is well-suited for systemic distribution and may serve as a promising candidate for oral administration.

Additionally, compounds such as stigmasterol from Angelica sinensis and Panax notoginseng, quercetin from Astragalus membranaceus, Angelica sinensis and Panax notoginseng, along with jaranol, hederagenin, isorhamnetin, and bifendate from Astragalus membranaceus, collectively contribute to the therapeutic benefits of MDBD. Stigmasterol from displays various pharmacological properties, including anti-diabetic, anti-tumor, anti-inflammatory, and antioxidant effects, and has been shown to inhibit high glucose-induced proliferation and angiogenesis in retinal cells (45, 46). Quercetin reduces inflammation and enhances retinal layer thickness, offering therapeutic benefits in DR (47). Jaranol exhibits antitumor activity against breast and liver cancers and possesses anti-influenza properties (48, 49). Hederagenin demonstrates antitumor and anti-inflammatory effects, and has been shown to attenuate high glucose-induced fibrosis in renal cells (50–52). Isorhamnetin provides cardiovascular protection and boosts insulin secretion from pancreatic β-cells (53, 54). Bifendate is recognized for treating hepatitis, reducing hepatotoxicity, and showing antifibrotic effects (55, 56). Formononetin is reported to ameliorate type 2 diabetes progression and related complications by lowering hyperglycemia and insulin resistance in diabetic rats (57). Calycosin, a phytoestrogen, exhibits potent antioxidant and anti-tumor properties, mitigating kidney injury and cognitive impairments induced by diabetes in animal models (58–60). Kaempferol regulates angiogenesis, apoptosis, metastasis, and inflammation, enhancing diabetes management through modulation of endoplasmic reticulum stress (61, 62). Finally, Ginsenoside Rh2 shows anti-cancer properties against various cancers and improves immune function, effectively reducing elevated fasting blood glucose levels in type 1 diabetes mellitus rats (63, 64). These findings show MDBD’s multi-target therapeutic potential in DR, combining active compounds to address inflammation, oxidative stress, and angiogenesis, thereby providing a comprehensive approach to managing this complex condition.

The core targets were identified using the PPI network, followed by differential expression analysis to pinpoint SDECGs core genes, which were subsequently validated in human samples. Among the top five SDECGs, ERBB2 emerged as a significant target based on feature importance scores, demonstrating a strong association with increased risk for DR. Beta-sitosterol was found to bind to ERBB2, a receptor tyrosine kinase linked to elevated levels that significantly correlate with a higher incidence of diabetes mellitus. Individuals with high ERBB2 levels exhibit a markedly increased risk of diabetes compared to those with lower levels (65). To validate the relationship between ERBB2 and DR, in vitro experiments were conducted using a hypoxia-induced DR cell model. Results revealed significant upregulation of ERBB2 at both mRNA and protein levels. However, treatment with beta-sitosterol (12 μM and 16 μM) and ERBB2 inhibition effectively downregulated these elevated levels, affirming the interaction between ERBB2 and DR. These findings suggest that beta-sitosterol may target the MDBD-related gene ERBB2 to mitigate key pathological processes in DR. Additionally, significant alterations in inflammatory and angiogenic factors associated with DR were observed. The proinflammatory cytokine IL-1β, markedly elevated in the hypoxia-induced model group, was significantly reduced following beta-sitosterol treatment and ERBB2 inhibition, highlighting ERBB2’s role in modulating inflammation. Similarly, angiogenic factors VEGF and ANGPTL6, crucial for pathological neovascularization, were modulated. While their levels were initially increased in the hypoxia-induced model, further downregulation occurred under beta-sitosterol treatment and ERBB2 inhibition, indicating ERBB2’s involvement in angiogenesis regulation in DR.

Apart from ERBB2, the other top four SDECGs identified in our analysis include CCND1, SERPINE1, TP53, and INSR. Each of these genes plays a distinct role in the pathophysiology of DR. While ERBB2 is associated to DR risk and its strong binding affinity with beta-sitosterol, which supports the in vitro findings, the inclusion of these additional targets demonstrates their complementary roles in DR pathogenesis and reflects a more comprehensive approach to understanding the DR mechanisms. CCND1, a protein-coding gene activated by insulin, facilitates glycemic normalization via the CCND1-CDK4 pathway (66). SERPINE1 (also known as PAI-1), part of the serine protease inhibitor family, has shown elevated expression in vitreous biopsies and neovascular tissue from patients with PDR (67). TP53, a tumor suppressor protein, exhibited higher expression in blood samples from T2DM patients with DR compared to healthy individuals and T2DM patients without DR (68).

IINSR, an insulin receptor, was identified in a genome-wide association study as strongly correlated with DR in T2DM patients (69). Differential expression analysis in our study further revealed high INSR expression in the DR group, showing its significance in insulin signaling and metabolic regulation. The apparent contradiction in INSR expression across different studies warrants further investigation. Previous studies have indicated that impairments in INSR expression or functionality can lead to insulin resistance and DM (70, 71). To investigate this apparent contradiction, further analysis was conducted. Muscle tissues from DR patients and healthy individuals were analyzed, differing from detection materials in prior studies, such as skeletal muscle (71), peripheral blood lymphocytes (72), and various cell lines, including CEM T-lymphocytes and SW1990 pancreatic cells (73). Additionally, gene-by-environment interactions were found to play a critical role in complex gene regulation. Controlled environmental exposure in vitro across different cell types revealed diverse transcriptional responses (74). This variation in INSR expression across cell types may explain its complex regulation within different cellular environments.

By incorporating secondary targets alongside ERBB2, this study adopts a broader perspective on DR pathogenesis. These additional targets provide valuable context about the effect of MDBD, complementing the primary findings without detracting from the significance of ERBB2. Moving forward, a systematic and stepwise evaluation of these targets is essential. Such prioritization should integrate factors including network-level impact, confidence in target–ligand interactions, and functional assessments in retinal endothelial cells under high-glucose conditions.

The ssGSEA analysis revealed significant differences in immune cell expression between normal and DR groups. Positive correlations were observed between SDECGs and immune cells such as eosinophils, macrophages M0, and resting CD4 memory T cells, while negative correlations were found with naive B cells and plasma cells. Among these, macrophages M0 and resting dendritic cells exhibited statistically significant differences. Macrophages M0, precursors to polarized macrophages, can be induced by high glucose levels to polarize into the M1 subtype (75), playing a key role in DR progression (76). Similarly, dendritic cells, which enhance immune responses via antigen presentation, have been identified as a DR risk factor through HLA DR expression in genome-wide association studies (GWAS) (77). Further analysis revealed that macrophages M0 positively correlated with two SDECGs, TP53 and SERPINE1, while resting dendritic cells showed a positive correlation with INSR and a negative correlation with ERBB2. All these correlations demonstrated statistical significance between DR and normal samples. Resting mast cells, on the other hand, were negatively correlated with four SDECGs, TP53, SERPINE1, ERBB2, and CCND1, while positively correlated with INSR, all with statistical significance. However, no significant differences in resting mast cell levels were observed between DR and normal samples, which may be attributed to the distinct and complex regulatory mechanisms between resting mast cells and genetic samples, leading to varied outcomes.

The RF model was utilized to identify feature genes, including the top five for construction, to assess DR onset and the sensitivity and accuracy of MDBD treatment, which showed high sensitivity and accuracy. Enrichment analysis indicated that MDBD mechanisms in DR treatment are linked to immunity and inflammation, aligning with DR pathogenesis and development factors (78). In DR, elevated blood glucose levels lead to mitochondrial dysfunction, inflammation, and increased vascular endothelial growth factor secretion, causing vascular and neuronal apoptosis, and neovascularization (79). The activation of immune cells takes place earlier than neuronal dysfunction and intraretinal microvascular abnormalities (80). Immune cell activation contributes to neuronal dysfunction and intraretinal microvascular abnormalities (81, 82). Elevated inflammatory cytokines such as IFN-γ, TNF, and IL-2, secreted by lymphocytes, have been detected in DR patients (83) and in diabetic rat retinal tissue (84). Additionally, macrophages have been identified in the retinal tissue of Akimba mice, a recognized model of diabetic retinopathy, through single-cell RNA sequencing (85). These results reveal the complex relationship between immune cell activation, genetic factors, and inflammatory processes in the development of DR. The correlations identified between specific immune cells and SDECGs point to the significant role of immune regulation in disease progression. The ability of MDBD to influence these immune-related pathways suggests its therapeutic relevance for DR. Further investigation into the interaction between immune cells and genetic factors will be essential for advancing our understanding of DR and improving treatment strategies.