Identifying Therapeutic Targets and Potential Drugs for Diabetic Retinopathy: Focus on Oxidative Stress and Immune Infiltration

Introduction

Diabetic retinopathy (DR) represents a prevalent microvascular complication associated with diabetes and stands as a principal contributor to preventable vision impairment within the working-age population.¹ DR affects 22.27% of diabetics worldwide, with 103 million individuals impacted in 2020 and 165 million by 2045.² Diabetic retinopathy has emerged as a pressing global public health concern, imposing significant financial burdens on both society and the economy. Current treatments for DR, such as including intravitreal pharmacologic agents, laser photocoagulation and vitreous surgery. Anti-VEGF therapies, including intravitreal injections of bevacizumab or ranibizumab, reduce retinal edema and neovascularization by inhibiting vascular endothelial growth factor (VEGF). However, these therapies offer limited improvement in vision for DR patients, and their efficacy tends to diminish over time, requiring frequent injections.³ Moreover, these treatment modalities fail to address the underlying causes of DR. Consequently, it is essential to clarify the mechanisms of diabetic retinopathy to discover novel diagnostic markers for early intervention and new targets for clinical pharmacotherapy.

Oxidative stress arises from a disparity between the production and elimination of free radicals.⁴ High blood glucose levels lead to an excess of oxygen free radicals in the body, causing damage to retinal tissues.⁵ Moreover, oxidative stress can induce inflammatory responses, vascular abnormalities, and cell apoptosis, further promoting the development of DR.^6,7 The pathogenesis of DR is heavily influenced by oxidative stress, according to many studies.^8–10 Thus, investigating the connections between oxidative stress and DR is crucial to developing new treatments and enhancing patients’ visual acuity and quality of life. The progression of DR is a complex process closely intertwined with the immune system, where immune cells, mediators, and the complement system, especially macrophages and microglial cells, play critical roles.^11–14 Previous research has highlighted the significance of the immune response and oxidative stress in the progression and aggravation of DR.

In recent years, bioinformatics approaches have become invaluable for identifying potential therapeutic targets and biomarkers by analyzing large-scale genomic and transcriptomic data. The OS-related genes were extracted from the GeneCards database, and RNA-sequencing datasets were obtained from the GEO database. We initially identified DEGs in the DR and control groups. This was followed by the extraction of DEOSGs from the intersection of DEGs and OS-related genes. For greater clarity of the DEOSGs, we ran functional and pathway enrichment analysis on them. We developed a PPI network to figure out candidate hub genes and subsequently utilized Least Absolute Shrinkage and Selection Operator (LASSO) regression, in conjunction with Support Vector Machine-Recursive Feature Elimination (SVM-RFE) algorithm analysis, to pinpoint hub genes associated with oxidative stress. These approaches enabled us to pinpoint hub genes that play critical roles in oxidative stress pathways, which are potential therapeutic targets. Nomogram creation and ROC curve evaluation assessed the diagnostic value of hub genes. This strategy is reliable for assessing these genes as biomarkers, improving DR diagnostic tools. A single-gene GSEA was conducted. Additionally, the ImmuCellAI website was used to examine the immune microenvironment in DR and determine hub gene-immune cell infiltration association. To better understand the regulatory features of hub genes, we constructed their regulatory networks in this study. Moreover, computational techniques like molecular docking are used to predict the binding affinity of drugs to identified targets, facilitating drug discovery. Consequently, drug prediction and molecular docking analyses were conducted to investigate the therapeutic potential of the identified targets, aiming to uncover novel drug candidates and understand their interactions at the molecular level. Finally, we validated the oxidative stress-related hub genes using both dataset analysis and in vivo animal models. Our research enhances comprehension of the molecular mechanisms behind oxidative stress in diabetic retinopathy pathogenesis, facilitates the discovery of novel diagnostic biomarkers, and reveals prospective therapy targets. Given the limited effectiveness of current therapeutic strategies for DR, the findings from this study may offer promising avenues for the development of more targeted and effective treatments. These biomarkers and drug candidates have the potential to be translated into clinical applications, improving early detection and therapeutic interventions, and ultimately preventing or delaying vision loss in DR patients.

Materials and Methods

Dataset Sources and Processing

This work utilized the GEO database to obtain the transcriptome dataset pertaining to diabetic retinopathy. Analysis was conducted utilizing the DR related dataset GSE160306. Considering that the most prevalent causes of vision loss in people with diabetes are diabetic macular oedema (DME) and proliferative diabetic retinopathy (PDR), the data from different groups’ gene expression profiles: healthy control (20 samples) and NPDR (Non-Proliferative Diabetic Retinopathy)/PDR + DME (20 samples) were analyzed.^15,16 For the external validation dataset, we used GSE102485. From the GeneCards database, we extracted 1,398 genes associated with oxidative stress for further analysis.¹⁷ These genes had a relevance score of 7 or higher.

In the GSE160306 dataset, we selected protein-coding genes for further analyses. Gene length was defined according to the longest transcript, and TPM (transcripts per million) data were obtained. The flowchart illustrating the methodology of this study is presented in Figure 1.

Figure 1 Research flowchart.

Identification of DEGs and Oxidative Stress-Related DEGs

DEGs between 20 DR patients and 20 healthy control samples were identified using the DESeq2 R package. The screening criteria, adjusted p-value (adj. P) < 0.05 and |log2 Fold Change| > 0.58, were set to balance sensitivity and specificity, minimizing false positives and reliably identifying DEGs, following commonly accepted bioinformatics standards. For the |log2 Fold Change| > 0.58 threshold, this corresponds to a minimum 1.5-fold change in gene expression. DEOSGs were derived by intersecting the discovered DEGs with 1,398 OS-related genes. The results were visualized by generating a volcano plot with the ggplot2 package utility. The VennDiagram and ComplexHeatmap packages were used to generate a Venn diagram and heatmap, displaying the expression of DEOSGs.

GO and KEGG Enrichment Analysis of DEOSGs

The R packages clusterProfiler and GOplot are utilized to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses to elucidate gene functions and probable pathways.^18,19 Bubble and bar plots were employed to display the 15 most significant GO terms and 10 most significant KEGG terms, arranged according to GeneRatio. To show the top 10 KEGG pathways according to z-scores, a circular plot was made. Moreover, the criteria for screening significant functions and pathways were adjusted to adj. P <0.05. The outcomes are shown with the ggplot2 R package.

Hub Genes Identification Combining PPI Networks and Machine Learning

Proteins were integrated into PPI networks according to their interactions with other proteins. The DEOSGs were evaluated for PPI utilizing the STRING database (http://string-db.org; version 11.5). The subsequent interactions within the network were illustrated utilizing Cytoscape. To identify candidate hub DEOSGs, we employed two Cytoscape plugins: CytoHubba and MCODE. For the MCODE plugin, we utilized default parameters. The highest-scoring network was selected. Furthermore, we evaluated the top 10 genes with the highest scores using the MCC algorithm from the CytoHubba plugin. Following this, these genes had been combined with the highest-scoring network to identify candidate hub genes.

We utilized two different machine learning analytics approaches, namely LASSO logistic regression analysis and SVM-RFE analysis, in order to increase the accuracy of selecting target hub genes. These methods were implemented to screen hub genes. The outcomes of the two machine learning analyses were subsequently intersected, leading to the identification of overlapping genes as the primary OS-related hub genes associated with DR. The ggplot2 and Venn Diagram packages were used to display the results graphically.

Nomogram Construction and ROC Curve Evaluation

The rms package in R was used to produce the nomogram. The model’s accuracy was assessed by the analysis of calibration curves. The evaluation of the ROC curve was conducted on the GSE160306 dataset to validate the hub genes. These genes were considered to have diagnostic significance if their area under the ROC curve exceeded 0.7. ROC analysis with pROC and ggplot2 visualization.

Single-Gene GSEA for Hub Genes in DR

We used clusterProfiler and enrichplot to perform a single-gene GSEA analysis to elucidate significant pathway in DR. Based on expression levels, hub genes in the DR group were categorised as high- or low-expression. Significant gene enrichment was defined as p-value <0.05.

Correlation Analysis Between Hub Genes and Immune Cell Infiltration

The normalized gene expression matrix derived from the GSE160306 dataset was employed for subsequent analysis of immune infiltration. To facilitate a study on immune infiltration, the normalized gene expression matrix was sent to ImmuCellAI (https://guolab.wchscu.cn/ImmuCellAI)²⁰ In order to compare the groups, the Wilcoxon rank sum test was employed. Using R’s ggplot2 package, Spearman correlation analysis was performed on infiltrating immune cells and hub genes. Grouped violin graphs showed DR and control immune cell abundance disparities.

Regulatory Networks Construction and Potential Drugs Prediction

Transcription factors and miRNAs were predicted using the JASPAR and TarBase databases via NetworkAnalyst (https://www.networkanalyst.ca/)^21,22 Uploading the identified hub genes to the DSigDB using the Enrichr platform (https://amp.pharm.mssm.edu/Enrichr/) for potential drug prediction analysis.²³

Molecular Docking

We used the OPENBABEL software to convert the 3D structures of the anticipated medicinal candidates from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) to PDB format. The RCSB PDB (https://www.rcsb.org/) was used for acquiring the crystal structures of the hub genes proteins. These protein structures were processed through the open-source version of PyMOL and AutoDockTools 1.5.7.²⁴ Subsequently, AutoDockTools 1.5.7 was utilized to generate a sufficiently large docking box for each target, covering the entire macromolecular protein to guarantee that all potential functional pockets were incorporated, facilitating the molecular docking of the candidate drugs.

To conduct the molecular docking, QuickVina-W was utilized.²⁵ The conformations with the lowest binding energies were chosen after the binding energies of the molecules had been evaluated. The complexes of the target proteins and drugs were then outputted in PDB format using PyMOL. Finally, all complexes were uploaded to the PLIP web tool for the analysis of additional non-covalent interactions between the target proteins and drugs.²⁶ Molecular docking results were visualized with PyMOL.

Dataset Validation of Hub Genes

We normalized gene expression data from the GSE102485 dataset using the DESeq2 R package. An unpaired t-test compared gene expression levels between diabetes and normal groups. The ggpubr package was used for visualization.

Establishment of Animal Models

From Jiangsu Huachuang Sino Pharma Technology (Jiangsu, China), 8-week-old male C57BL/6J mice were obtained. Standard rodent feed for mice. To create the Streptozotocin (STZ) solution, the drug was dissolved in a sodium citrate buffer with a concentration of 0.1 mol/L and a pH of 4.2. Five days of 50 mg/kg STZ intraperitoneal injections were given to mice. Blood glucose levels were assessed via tail vein sampling seven days after the final injection. A diabetic mouse model was deemed successfully established if blood glucose exceeded 16.7 mmol/L. Untreated age-matched normal control mice were utilized.

Real-Time Quantitative Polymerase Chain Reaction

TRIzol was utilized in order to successfully isolate whole cellular RNA at the desired concentration. Total RNA was reverse-transcribed into complementary DNA using Roche Premix for quantitative PCR (qPCR) from Roche and Accurate Biology. Subsequently, RT-qPCR was performed using the TB Green Fast qPCR Mix (Takara). The forward and reverse primers can be located in Supplementary Table 1. The 2-ΔΔCT method, normalized to β-actin levels, was employed to assess the expression levels of target mRNA.

Western Blot

Retinal tissue protein was extracted according to manufacturer directions (Beyotime Biotechnology, China). The protein samples underwent denaturation and were then separated via SDS-PAGE. Standard procedures were then used to put them onto PVDF membranes (G2154-1L, Servicebio, China). Primary antibodies targeting CCL4 (1:1000, Bioss, bs-2475R, China), FCGR2B (1:1000, Abclonal, A12553, China) and β-actin (1:1000, ZSGB-BIO, TA-09, China) overnight at 4°C. ImageJ software was used to evaluate band intensity, with β-actin acting as an internal reference.

Statistical Analysis

Each bioinformatics analysis was executed with R software. The data was analyzed using GraphPad Prism 9.5 and is shown as the mean ± SEM from three separate experiments. The degrees of significance are stated as *** P<0.001, ** P<0.01, * P<0.05.

Results

Identification of DEOSGs

Among the 772 DEGs that were found in the GSE160306 dataset, 588 were revealed to be upregulated and 184 were discovered to have decreased expression in the DR group. A representation of the distribution of DEGs may be seen on the volcanic map (Figure 2A). To further identify DEOSGs, we intersected 1398 oxidative stress-related genes with the DEGs. This study found 91 DEOSGs in total, with 75 of them being upregulated and 16 of them being downregulated in DR. We defined oxidative stress-related DEGs as DEOSGs (Figure 2B). The heatmap illustrates that the relative gene expression of DEOSGs varies between the DR and control groups (Figure 2C).

Figure 2 Identification of DEOSGs. (A) Volcano plot of DEGs in GSE160306. (B) Venn diagram of oxidative stress-related genes and DEGs in GSE160306. (C) Heat map of DEOSGs. Abbreviations: DEGs, differentially expressed genes; OS, oxidative stress.

Functional Enrichment Analyses of DEOSGs

Enrichment analysis of the differentially expressed DEOSGs was conducted to better understand their molecular functions and signaling pathways. GO enrichment study indicated that in the biological processes (BPs), DEOSGs were predominantly enriched in response to oxidative stress, oxygen levels and decreased oxygen levels (Figure 3A). In the cellular components (CCs), DEOSGs are primarily enriched in the endoplasmic reticulum lumen, neuronal cell body, and early endosome (Figure 3A). In molecular functions (MFs), activity as a signaling receptor activator, receptor ligand, enzyme inhibitor and peptidase regulator are most abundant in DEOSGs (Figure 3A). According to the findings of KEGG enrichment analysis, DEOSGs were shown to be significantly enriched in the PI3K-Akt signaling pathway as well as Cytokine-cytokine receptor interaction (Figure 3B). Analysis using z-scores indicated that DEOSGs exhibited significant enrichment in the PI3K-Akt signaling pathway, Staphylococcus aureus infection, Hepatocellular carcinoma, and Complement and coagulation cascades (Table 1 and Figure 3C).

Table 1 The Top 10 KEGG Pathways According to z-Scores

Figure 3 The functional enrichment analyses of the DEOSGs. (A) represents the GO enrichment analysis plot of DEOSGs. (B) represents the KEGG enrichment analysis plot of DEOSGs. (C) KEGG analysis of the DEOSGs. The inner ring is a bar plot where height displays the significance of the term, and the outer ring displays scatter plots which indicate the expression levels (log FC) for the genes in each term. Abbreviations: BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Hub Genes Identified Combining PPI and Machine Learning

Applying the STRING database, we conducted a study of the PPI for the 91 DEOSGs. Cytoscape was used to show the results as a network (Figure 4A). Identifying significant gene clusters was accomplished with the usage of the MCODE plug-in. One particular module, comprising 12 nodes and 43 edges, was categorized as significant and includes genes such as ALOX5, CD4, CR2, CSF3, FCGR2B, FOXP3, ITGAL, ITGAM, SELL, TFRC, TIMP1, TLR3 (Figure 4B). By utilizing the MCC algorithm within the CytoHubba plug-in, the 10 highest-scoring genes were determined, including CCL4, CD4, CSF1, CSF3, CXCR4, FOXP3, ITGAM, SPP1, TP53, and VCAM1 (Figure 4C). Combining the results, 18 candidates of oxidative-stress-related hub genes were identified, including CCL4, CSF1, CXCR4, SPP1, TP53, VCAM1, ALOX5, CD4, CR2, CSF3, FCGR2B, FOXP3, ITGAL, ITGAM, SELL, TFRC, TIMP1, TLR3 were eventually obtained.

Figure 4 Screening candidate hub genes by Protein-Protein Interaction (PPI) Network. (A) PPI of DEOSG. (B) A key cluster with 12 genes was further chosen as candidate hub genes by MCODE. (C) Top 10 candidate hub genes explored by CytoHubba.

Initially, 6 candidate hub genes were discovered via LASSO regression (Figure 5A), while the SVM-RFE algorithm selected 13 candidate hub genes (Figure 5B). Intersecting these results with a Venn diagram revealed four hub genes: CCL4, CR2, FCGR2B, and FOXP3 (Figure 5C).

Figure 5 Screening hub genes by machine learning. (A) LASSO regression algorithm. (B) SVM-RFE algorithm. (C) Venn diagrams for two algorithms. Abbreviations: Lasso, Least Absolute Shrinkage and Selection Operator; SVM-RFE, Support Vector Machine-Recursive Feature Elimination.

Diagnostic Value Assessment

To improve the accuracy of predicting DR progression, a nomogram with four hub genes was developed (Figure 6A). The calibration curve results confirmed the nomogram’s high predictive reliability for DR patient outcomes (Figure 6B). The diagnostic efficacy of the 4 hub genes was assessed using ROC curve analysis; promising diagnostic markers were hub genes with an AUC value higher than 0.7. Results revealed that within the GSE160306 dataset the AUC values for FOXP3, FCGR2B, CR2, and CCL4 in relation to DR were 0.887, 0.855, 0.840, and 0.780 correspondingly (Figure 6C). The findings suggest that the 4 hub genes possess significant diagnostic potential for DR.

Figure 6 Nomogram construction and the diagnostic value evaluation. (A) The visible nomogram for diagnosing DR. (B) Calibration curves of nomogram for predicting DR the GSE160306 dataset. (C) ROC curves of the hub genes in the GSE160306 dataset. Abbreviation: AUC, area under the curve.

Single Gene GSEA of Hub Genes

We performed a single-gene GSEA-KEGG analysis to clarify the potential pathways related to the hub genes. The top 6 pathways that are significantly enriched are displayed based on the p-value. CCL4 was found to be positively correlated with pathways related to Alpha-linolenic acid metabolism, Leishmania infection and T cell receptor signaling pathway. Conversely, it showed a negative correlation with pathways related to Oocyte meiosis and Progesterone-mediated oocyte maturation (Figure 7A). CR2 exhibited positive correlations with T cell receptor signaling, Allograft rejection, Gap junction, Glutathione metabolism, and Pathogenic Escherichia coli infection, but negatively correlated with Oxidative phosphorylation (Figure 7B). FCGR2B showed positive associations with Hematopoietic cell lineage and Natural killer cell-mediated cytotoxicity. FCGR2B showed positive associations with immune pathways, indicating its potential role in immune regulation (Figure 7C). FOXP3 demonstrated negative correlations with Tyrosine and Tryptophan metabolism, Allograft rejection, Melanoma, Oxidative phosphorylation, and Glyoxylate and dicarboxylate metabolism, suggesting its role in immunomodulation (Figure 7D).

Figure 7 Single-gene GSEA of hub genes. (A-D) represents the signaling pathway associated with the hub genes. Abbreviation: KEGG, Kyoto Encyclopedia of Genes and Genomes.

Immune Cell Infiltration and Correlation of Hub Genes With Immune Cells

Enrichment analysis revealed that DEOSGs were mostly associated with signaling pathways that pertain to the immune system. Consequently, the immune infiltration in DR was investigated using the ImmuCellAI algorithm. Comparing the DR and control groups, the results revealed that 12 distinct types of immune cells were significantly different (Figure 8A and B). It was noted that DC cells, Macrophage, NK cells, CD4 T cells, Tr1 cells, nTreg cells, iTreg cells, Th1 cells and Th17 cells were mainly enriched in DR. However, in the control group, B cell, Neutrophil and Mucosal-Associated Invariant T cells (MAIT) were enriched. To discover the association between hub genes and immune cell infiltration, we employed Spearman correlation analysis (Figure 8C). Analyses of correlations showed that CCL4 correlated positively with NK cells and Macrophage cells and negatively with CD8(+) T cells and Mucosal-Associated Invariant T cells. CR2 correlated most positively with B cells and negatively with Effector Memory Cells. FCGR2B showed the most positive correlation with Macrophage cells and the most negative correlation with B cells. Although FOXP3 was negatively correlated with B cells, it was positively correlated with Natural Regulatory T cells.

Figure 8 Immune cell infiltration and correlation analysis of Hub genes with immune cells. (A) represents the expression of 24 types of immune cells in the training set GSE160306 of the DR and control group. (B) represents a stacked plot of the expression of 24 types of immune cells in each sample. (C) represents the correlation analysis plot between hub genes and 24 types of immune cells.

Construction of Regulatory Networks and Drug Prediction

Through the JASPAR database, 24 TFs were identified, among which 9 exhibited a degree ≥ 2, including PPARG, NR3C1, STAT3, HNF4A, TFAP2C, YY1, GATA2, USF2, FOXC1 (Figure 9A). Possible miRNAs were identified through the TarBase database, revealing 13 miRNAs with a degree of ≥2 (Figure 9B). In this study, potential therapeutic compounds targeting the hub genes were predicted using the DSigDB database. On the basis of the adjusted p-values, the top five possible chemical compounds were determined to be shown in Table 2. The results revealed that (+)-chelidonine (PubChem CID 197810), oxazolone (PubChem CID 1712094), and eugenol (PubChem CID 3314) were the three most significant compounds associated with CCL4 and FCGR2B. Additionally, AGN-PC-0JHFVD (PubChem CID 71581418) was linked to FCGR2B and FOXP3, while simvastatin (PubChem CID 54454) was associated with CCL4 and FOXP3. The compound AGN-PC-0JHFVD, a compound identifier from Angene Chemical (https://www.angenechemical.com/), is identified by PubChem CID 71581418 and has a molecular formula of C₂₆H₂₄Cl₂N₄O₅S₂.

Table 2 Candidate Drug Predicted Using DSigDB

Figure 9 Regulatory networks. (A) Interaction network of TFs and genes for the hub genes. (B) Interaction network of miRNAs and genes for the hub genes.

Molecular Docking Analysis

Molecular docking scores between the candidate drugs and target proteins are presented in Table 3. A stronger and more desirable binding relationship is suggested by a lower negative binding energy. Requiring less than −7 kcal/mol for binding energy suggests a stable connection.^27,28 The binding energies of (+)-chelidonine with both CCL4 and FCGR2B were less than −7 kcal/mol. Additionally, AGN-PC-0JHFVD exhibited binding energies less than −7 kcal/mol with FCGR2B and FOXP3, indicating extremely stable binding. After analyzing noncovalent interactions, details including hydrogen bonds, hydrophobic interactions, π-stacking, and relative amino acid residues are exhibited in Figure 10.

Table 3 Docking Results of Available Proteins With Candidate Drugs

Figure 10 Molecular docking simulation diagram. (A) (+)-chelidonine binds to CCL4. (B) oxazolone binds to CCL4. (C) eugenol binds to CCL4. (D) simvastatin binds to CCL4. (E) (+)-chelidonine binds to FCGR2B. (F) oxazolone binds to FCGR2B. (G) eugenol binds to FCGR2B. (H) AGN-PC-0JHFVD binds to FCGR2B. (I) AGN-PC-0JHFVD binds to FOXP3. (J) simvastatin binds to FOXP3. CCL4 is marked in red, FCGR2B in wheat, FOXP3 in green, and candidate drugs in cyan. Blue solid lines represent hydrogen bonds, grey dotted lines represent hydrophobic interactions, and green dotted lines represent π-stacking interactions.

Dataset Validation

The validation of the dataset revealed that the levels of CCL4 and FCGR2B expression were greater in DR samples compared to control samples (Figure 11A-D). This conclusion is in line with the findings that were obtained from the GSE160306 dataset.

Figure 11 Dataset validation of the Hub genes. (A-D) Expression differences of the hub genes in GSE102485 dataset. *** represents P<0.001, ** represents P<0.01 and ns, no significance.

Verification of Relative Expression of Hub DEOSGs by RT-qPCR and WB

To validate the bioinformatic analysis findings, we analyzed the expression of CCL4 and FCGR2B in STZ-induced 8-week diabetic mice retina samples. RT-qPCR and WB were performed on retinal samples from these diabetic mice. RT-qPCR determined that CCL4 expression was substantially higher in DR samples than in control samples; nonetheless, FCGR2B expression was not different between the groups (Figure 12A). Consistently, WB analysis showed a marked elevation of CCL4 protein levels in DR samples, while FCGR2B protein expression was similar between DR and control groups (Figure 12B).

Figure 12 External validation of CCL4 and FCGR2B. (A) The mRNA levels of CCL4 and FCGR2B were evaluated in animal models by RT-qPCR. (B) The protein levels of CCL4 and FCGR2B were evaluated animal models by WB. ** represents P<0.01, * represents P<0.05 and ns, no significance. Abbreviation: STZ 8w, Streptozotocin 8 weeks.

Discussion

DR is triggered by a confluence of circumstances.^29,30 Nonetheless, the precise processes have yet to be comprehensively elucidated. Limited therapeutic options are currently available for DR, highlighting the urgent need to discover new molecular pathways to help treat and diagnose this disease. Inflammation, vascular dysfunction, and oxidative stress have all been shown to have a major effect on the course of DR.^6,31,32 The onset of diabetes-associated retinal diseases is significantly influenced by immune cell infiltration. The pathophysiological process of DR involves immune reactions such as leukocyte aggregation, neutrophil and macrophage infiltration, as well as complement and microglia activation in the retina.^11,33,34 Oxidative stress is considered to be a crucial component of the pathophysiological processes of DR, potentially interacting with pathological processes such as inflammation.^7,9,35 The molecular mechanism of oxidative stress in DR, however, continues to be ambiguous and requires additional investigation.

In this study, for the first time, a transcriptome dataset was retrieved from the GEO database using bioinformatics methods, identifying 772 DEGs under the specified threshold, with 588 exhibiting upregulation and 184 demonstrating downregulation. To focus on oxidative stress, we intersected these DEGs with oxidative stress-related genes, identifying 91 overlapping genes, which we designated as DEOSGs. Further GO and KEGG enrichment analyses of DEOSGs have shown their participation in significant biological processes and pathways, particularly those associated with oxidative stress and the PI3K-Akt signaling pathway. At various stages in the DR process, the PI3K-Akt signaling pathway displays distinct activation levels.³⁶ Research indicates that in early stages DR, the PI3K/Akt/mTOR signaling pathway can decrease reactive oxygen formation and augment antioxidant capacity, thus diminishing the apoptosis of retinal pericytes.³⁷ In the later stages, there is a significant contribution made by the VEGFA/PI3K/Akt pathway to the process of compensatory angiogenesis in DR.^38,39 To further narrow down the key players, we employed PPI network analysis in combination with machine learning techniques, which allowed us to identify hub genes with potential significance in DR. PPI networks were constructed and combined with machine learning algorithms to determine a set of four hub genes related to oxidative stress, namely CCL4, CR2, FCGR2B, and FOXP3. CCL4 (MIP-1β) is a C-C motif chemokine ligand synthesized by many cells, especially distinct immune cell types.⁴⁰ CCL4 plays a crucial role in recruiting inflammatory cells to migrate towards inflamed tissues during inflammatory responses.⁴¹ Research indicates that stimulation by CCL4 significantly increases the production of reactive oxygen species. Additionally, CCL4 activates the PI3K-Rac1 cascade, which in turn promotes cell adhesion reactions between endothelial cells and monocytes.⁴² Furthermore, Dai et al demonstrated that that PDR patients had higher vitreous CCL levels.⁴³ Chang et al found that inhibiting CCL4 boosts the function of endothelial progenitor cells and facilitates angiogenesis triggered by ischemia in diabetic models.⁴⁴ According to research, CCL4 recruits monocyte lineage cells produced from bone marrow, which may be important for the physiological revascularization of hypoxic avascular retinas.⁴⁵ These studies suggest a significant correlation between CCL4 and both inflammation and DR. CR2, a membrane glycoprotein that bind activation and processing fragments of the complement system, has been suggested to be involved in immune dysregulation and inflammatory processes.^46,47 The complement system is intricately linked to the onset and advancement of DR.^48,49 Although there are limited studies directly linking CR2 to DR, CR2 regulates immune response and inflammatory pathways that affect vascular damage and neuroinflammation.^50–52 Furthermore, studies indicate that targeting complement inhibitors to sites of complement activation and C3d deposition can be achieved using CR2 fragments is effective in treating choroidal neovascularization (CNV) in age-related macular degeneration (AMD).^53–55 FCGR2B, Fc gamma receptor IIb, is a low-affinity immunoglobulin γ receptor in the Fc region and belongs to the fragment crystallizable receptor (FcR) family. FCGR2B, as a key regulator of the immune response, might play a role in the development of DR through its involvement in the activation of immune cells and inflammatory mechanisms.^56–58 Although there is limited evidence linking FCGR2B dysregulation directly to DR, studies have shown that it may be involved in autoimmune diseases and inflammatory disorders, which have been shown to be contributing factors to retinal damage and vascular dysfunction in DR.^59,60 FCGR2B has been implicated in modulating immune cell activation and may affect the immune response in diabetic retinopathy, thereby influencing disease progression.^61,62 CD4+ regulatory T cells (Tregs) are the main cells that express FOXP3 (Forkhead Box P3). Coordinated regulation of Treg formation, function, and homeostasis is accomplished by FOXP3 via its action on target genes.^63,64 Tregs are crucial for maintaining immune balance, including the maintenance of tolerance to self-antigens, minimizing excessive immune responses and maintaining tissue homeostasis.⁶⁵ The X-linked syndrome (IPEX), which stands for X-linked Immune Dysregulation, Polyendocrinopathy, and Enteropathy, is a rare and severe autoimmune disease. It results from mutations in the FOXP3 gene, disrupting the internal balance and functionality of Tregs.⁶⁶ FOXP3 exhibits conditional expression in macrophages within stroke lesions. FOXP3+ macrophages demonstrate enhanced clearance capabilities in the context of stroke lesions.⁶⁷ An increasing consensus suggests that Foxp3 plays a significant role in controlling internal balance, inflammation, and oxidative stress.^68–70

Subsequently, the hub genes were employed in the construction of a nomogram, accompanied by ROC curve analysis, for the purpose of assessing their diagnostic utility. The findings were further validated through this approach, thereby affirming the value of the hub genes as diagnostic indicators. The robust performance of these genes in diagnostic models underscores their potential as reliable biomarkers for the detection and risk stratification of DR.

Expanding upon these discoveries, we conducted single-gene GSEA to delve deeper into the functional functions of these hub genes. Additionally, we also investigated the relationships between the hub genes and different immune cells and performed immune infiltration study to better understand how oxidative stress interacts with immune response in DR. We employed the ImmuCellAI online tool for a comprehensive assessment of immune infiltration in DR. From an overall distribution perspective, compared to the control group, the DR group exhibits an increased proportion of DC cells, Macrophages, NK cells, CD4 T cells, Tr1 cells, nTreg cells, iTreg cells, Th1 cells, Th2 cells, and Th17 cells. Conversely, the proportion of B cells, Neutrophils, and MAIT cells decreases. Meanwhile, we performed Spearman correlation analysis to evaluate the relationship between hub genes and the various immune cell infiltrations. Studies have shown that there is a significant interaction between immune cells and DR.^71,72 T cells were the main infiltrating cells in DR Samples. Some studies have found that, Th1/Th2 imbalance seems to influence DR development, since Th1 cytokine release is elevated and Th2 secretion is decreased during DR.⁷³ Th1 and Th2 cytokines also play an important role in regulating angiogenesis.⁷⁴ Additionally, research found that Th17 cells infiltrate the retina of the mouse model of DR.⁷⁵ Possible correlation between DR and interleukin (IL)-17A) level disorders.⁷⁶ The transient increase of regulatory T cells in retinopathy can reduce neovascular retinopathy in mice.⁷⁷ There is a link between MAIT cells and metabolic disorders.⁷⁸ The quantity of MAIT cells in circulation is dramatically decreased in type 2 diabetic patients.⁷⁹

Recognizing the regulatory complexity of these hub genes, we established regulatory networks that included TFs and miRNAs, which offered deeper insights into the molecular mechanisms underlying DR. Drug prediction analysis identified potential therapeutic agents targeting these hub genes, further validated by molecular docking studies. In our study, drug screening and molecular docking analyses revealed potential interactions between (+)-chelidonine and CCL4 and FCGR2B, as well as between AGN-PC-0JHFVD and FCGR2B and FOXP3. The pharmacological actions of (+)-chelidonine, an alkaloid produced from the Chelidonium majus, include anti-inflammatory, antibacterial, and anticancer characteristics.^80–82 Previous research has shown that the efficacy of (+)-chelidonine in modulating related NF-κB pathways, making it a promising candidate for further investigation.^83,84 AGN-PC-0JHFVD is a compound with specific pharmacological properties yet to be fully characterized in the literature. Initial studies suggest that it may have significant therapeutic potential due to its strong binding affinity with key target proteins. The binding of AGN-PC-0JHFVD to FCGR2B and FOXP3 suggests a potential therapeutic application. Drug predictions are crucial for providing new treatments for DR. Docking scores and binding affinity demonstrated strong interactions, indicating that these drugs might effectively modulate the activity of key proteins. Targeting these proteins might be a promising therapy option for DR, according to this discovery.

The expression levels of hub genes were confirmed in the dataset relevant to diabetic retinopathy, demonstrating that only the findings for CCL4 and FCGR2B agreed with the analytical results from GSE160306. To confirm the outcomes of bioinformatics analysis, the CCL4 and FCGR2B were validated using RT-qPCR and Western blot. The results from RNA-seq did not match up with the results from RT-qPCR and Western blot for FCGR2B in this study. Potential reasons for these inconsistencies include variations in sample preparation, where RNA-seq utilizes whole transcriptome sequencing covering the entire transcriptome, while RT-qPCR selectively targets specific genes. This divergence may arise from differences in sample handling and preparation steps. Additionally, RT-qPCR’s selective amplification of specific gene fragments could introduce bias due to primer selectivity, contrasting with RNA-seq, which sequences the entire transcriptome, offering a more comprehensive view. Moreover, RNA-seq’s capacity to detect low-expression genes across the entire genome contrasts with RT-qPCR’s potentially higher sensitivity, which may, however, be compromised by technical noise, particularly at lower expression levels. It is important to note that the consistency between the Western Blot (WB) and RT-qPCR results suggests that FCGR2B does not exhibit differential expression in the diabetic mouse model. We ultimately identified CCL4 as oxidative stress-related biological markers. Chemokine-based therapies have shown promise in treating inflammatory diseases by targeting immune cell recruitment. For example, studies have demonstrated that treatment with mNOX-E36–3’PEG, an anti-CCL2 L-enantiomeric RNA aptamer, can improve renal function and reduce kidney damage in the db/db mouse model of type 2 diabetic nephropathy.⁸⁵ The CCL4/CCR5 axis is pivotal in immune cell recruitment to inflamed tissues.⁸⁶ In DR, targeting this pathway could potentially reduce retinal inflammation by limiting both oxidative stress and the infiltration of immune cells.

Our study has successfully identified oxidative stress-related biomarkers in DR for the first time. Nevertheless, it is imperative to recognize the inherent limitations that accompany this investigation. While it offers valuable preliminary insights into underlying mechanisms, it falls short of providing definitive evidence. Primarily, our methodology relied on database mining, lacking adequate experimental validation. Furthermore, our validation process was confined to mouse models, with a constrained sample size. It’s essential to note the significant genomic and biological differences between mice and humans. To bridge the gap toward clinical application, it is crucial to validate these findings using human clinical samples, such as vitreous humor or retinal tissues from DR patients. This will help establish a direct connection between the identified biomarkers and their role in disease pathophysiology. Moreover, incorporating these biomarkers into diagnostic panels or predictive models for DR progression may enhance their clinical utility. Developing therapeutic strategies that target the oxidative stress pathways linked to these biomarkers could pave the way for novel interventions. Further studies are warranted to explore how the identified transcription factors, miRNAs, and therapeutic agents could be integrated into clinical practice. These efforts would significantly enhance the translational impact of our study, providing a clear path from biomarker discovery to clinical intervention. Moreover, we did not proceed with validation following the successful prediction of transcription factors, miRNAs, and therapeutic agents. Lastly, there remains a necessity for further verification concerning the co-expression relationship between DEOSGs and immune cells.

Conclusion

This study conducted comprehensive bioinformatic analysis of gene activity in patients with DR, identifying 4 hub genes closely associated with oxidative stress in DR. Additionally, it revealed the relationship between CCL4, CR2, FCGR2B, and FOXP3 genes and immune cell infiltration. Experimental validation confirmed the importance of CCL4 as a biomarker for oxidative stress in DR, demonstrating its considerable clinical translational potential as both a diagnostic marker and a therapeutic target. Targeting CCL4 could represent a promising novel strategy for the treatment of DR. Overall, these findings improve our understanding of DR’s molecular processes in connection to oxidative stress and provide potential biomarkers for diagnosis and treatment.