Abstract Diabetic retinopathy (DR) is a common microvascular complication of diabetes and a leading cause of vision loss worldwide. Although several mechanisms have been implicated in the pathogenesis of DR, emerging evidence suggests a link between ferroptosis and DR. Unfortunately, the exact mechanism underlying this connection is not clear. Therefore, investigating the role of ferroptosis in diabetic retinopathy holds promise for advancing our understanding of this complex disease and developing innovative treatments. We have identified differentially expressed genes (DEGs) and differentially expressed marker genes (DEMGs) from open-source single-cell RNA sequencing datasets by using in depth in silico approach. Subsequently, ferroptosis-associated DEGs (FA-DEGs), ferroptosis-associated DEMGs (FA-DEMGs), and ferroptosis-associated Hub Genes (FAHGs) were identified. The FDA-approved drugs for our target proteins were also identified, and their ADMET properties were assessed. Molecular docking and simulation were utilized to explore the interaction stability of the compounds with the target proteins. Overall, we identified 63 FA-DEMGs that were significantly enriched in Peroxiredoxin activity, Ferroptosis, Mitophagy, and Autophagy. Further analysis predicted that PRDX1 and UBC are candidate target proteins. Molecular docking results showed that dexamethasone has a high binding affinity for both PRDX1 and UBC. Additionally, molecular dynamics simulations revealed that dexamethasone (which showed the best hit in the docking analysis) exhibited a ‘stable effect’ on both PRDX1 and UBC. To summarize, this study showed that PRDX1 and UBC could be suitable therapeutic targets for dexamethasone, which might be helpful in the advance of DR treatments in the future. Keywords: Diabetic retinopathy, Ferroptosis, Ferroptosis associated differentially expressed markers genes, Molecular docking, Molecular dynamics simulations, Therapeutic targets Highlights * • Ferroptosis-associated hub genes identified through single-cell RNA sequencing analysis. * • Identified marker genes associated with ferroptosis, oxidative stress, and inflammatory pathways. * • Fast molecular docking identified compounds targeting ferroptosis proteins. * • Dexamethasone identified as a potential drug for treating DR-associated ferroptosis. 1. Introduction Diabetic retinopathy (DR) is one of the most common microvascular complications associated with type I and type II diabetes. It causes irreversible vision loss and may result in complete blindness if left untreated [[31]1]. The alarming rate at which diabetes incidence is increasing has led to great concern about diabetes-associated complications among the scientific community. According to a recent meta-analysis, in 2020, 103.12 million (22.27 %) diabetic individuals suffered from DR, and 18.83 million (6.17 %) had vision-threatening DR. These figures are expected to increase to 160.50 million (55.6 %) and 44.82 million (57.0 %), respectively, by 2045 [[32]2]. DR progresses from initial non-proliferative diabetic retinopathy (NPDR) to sight-threatening proliferative diabetic retinopathy (PDR). NPDR exhibits mild, moderate, and severe phenotypes depending on the clinical hallmarks it possesses [[33]3,[34]4]. Although previous reports have linked DR to various pathways affecting mitochondria and the endoplasmic reticulum (ER), leading to aggravated oxidative stress and inflammation, the exact mechanism leading to disease pathogenesis remains unclear. Among the different factors that govern the pathogenesis of DR, oxidative stress plays a key role. Additionally, in response to extensive lipid peroxidation-mediated stress in membrane bilayers, cells commit to a nonapoptotic mode of cell death known as ferroptosis [[35]5]. The cells undergo mitochondrial condensation associated with the disappearance of cristae, followed by the rupture of the outer mitochondrial membrane [[36]6]. Ferroptosis is mainly driven by iron-dependent accumulation of lipid peroxides and inhibition of GPX4 (an enzyme that utilizes GSH to scavenge lipid ROS) [[37]7]. Furthermore, GPX4 can be inhibited either directly by the organic compound (1S,3R)-RSL3, indirectly by the inhibition of system Xc- (a cystine/glutamate antiporter) or indirectly by the inhibition of the rate-limiting enzyme in GSH production [[38]8]. Oxidative stress likely acts as a triggering factor for ferroptosis [[39]9]. Thus, oxidative stress may induce ferroptosis-mediated death of retinal cells and aggravate DR pathogenesis. Recent studies have established a connection between DR and ferroptosis. DR-associated ferroptosis has been shown to cause immense damage to pigmented epithelial cells, which leads to deterioration of the blood‒retinal barrier and an increase in the permeability of capillary endothelial cells in the human retina [[40]10]. Unfortunately, the metabolic pathways involved in ferroptosis in the context of DR pathogenesis are not yet fully understood. To alleviate the complications associated with DR, it is necessary to predict the therapeutic targets of ferroptosis. This study utilized in silico approach to identify differentially expressed marker genes (DEMGs) by analysing openly available single-cell RNA sequencing (scRNA-Seq) DR mouse model datasets. To obtain ferroptosis-associated differentially expressed marker genes (FA-DEMGs), a comparison of DEMGs with a freely available database of ferroptosis-associated genes (FerrDB database) was performed. In addition, pathway enrichment analyses were conducted to predict the pathways associated with FA-DEMGs, after which network analysis was carried out to identify ferroptosis-associated hub genes. To explore the interactions of the drug compounds with the therapeutic targets of DR-associated ferroptosis, we performed molecular docking, and the SwissADME and ProTox-II web servers were used to predict the pharmacokinetic and drug-likeness properties of the compounds. Finally, molecular dynamics (MD) simulations were performed to determine the dynamic behavior and interaction stability of the compounds with their target proteins. Overall, this study provides new insight into identifying drug targets and inhibitors to ameliorate ferroptosis-associated pathological processes in DR. 2. Materials and methods 2.1. Raw data acquisition Single-cell RNA sequencing (scRNA-seq) of DR datasets were retrieved from two different sequencing database archives. The E-MTAB-9061 (a total of 6 mice samples, of which 4 were Akimba (DR) and 2 were WT control) mouse dataset was downloaded from the ArrayExpress database at EMBL-EBI ([41]https://www.ebi.ac.uk/arrayexpress), while the PRJNA653629 (a total of 6 mice samples, of which 3 were db/db mice and 3 were db/m mice) dataset was retrieved from the European Nucleotide Archive (ENA) ([42]https://www.ebi.ac.uk/ena/browser/home) [[43]11,[44]12]. The E-MTAB-9061 data consisted of 10x Genomics-based single-cell transcriptome profiles of DR from wild-type control and Akimba DR mouse samples. Sequencing was carried out on an Illumina HiSeq 4000 platform with a paired-end configuration. In contrast, the PRJNA653629 dataset comprised a high-throughput single-cell transcriptional landscape of retinal cells in the DR from three control db/m mice and three leptin receptor-deficient diabetic db/db mice. The raw reads were generated on an Illumina NovaSeq 6000 platform with a paired-end configuration ([45]Table 1). Table 1. Summary of single cell DR datasets. Accession ID Platform Sample E-MTAB-9061 Ilumina Normal Retina- 02 Diabetic Retinopathy (DR)- 04 PRJNA653629 Ilumina Normal Retina- 03 Diabetic Retinopathy (DR)- 03 [46]Open in a new tab 2.2. Data preprocessing and quality control analysis The raw reads of scRNA-seq were downloaded from the database in Fastq file format, and the quality was checked with the help of FastQC tools ([47]https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). STAR aligner, a tool of Cell Ranger software (version 7.0, 10x Genomics, USA), was used to align the raw sequence with mm10, the reference genome of mice. The DEGs were identified with Cell Ranger software (version 7.0, 10x Genomics), which uses the “Cell Ranger count” function to establish single-cell feature counts in the form of features, barcodes, and single-cell expression matrix files. The data were further processed by the “Seurat package” (version 4.3.0; [48]https://cran.rproject.org/web/packages/Seurat/index.html) of “R” studio (version 4.2.1; [49]https://cran.r-project.org/bin/windows/base/old/4.2.1/). Briefly, in the Seurat package, the features, barcodes, and matrix files of both datasets were refiltered, and cells with a specific range of parameters [>200 genes and <8000 genes; >400 unique molecular identifiers (UMIs) and <5 % mitochondrial RNA (mtRNA)] were retained. The filtered matrix files obtained from both datasets were normalized to obtain the normalized counts with the help of the “LogNormalize” function of the Seurat package. While the “FindVariableFeatures” function was used to identify the highly variable features, the “ScaleData” function applied a linear transformation to the data before proceeding with dimensionality reduction or principal component analysis. The top 20 principal components were considered by the t-distributed stochastic neighbor embedding (t-SNE) technique in an attempt to reduce the dimensionality. Finally, cells were clustered using a graph-based clustering algorithm approach followed by the Louvain modularity optimization algorithm at a resolution of 0.5 [[50]13]. Annotation of different cell types was performed based on known retinal cell-specific markers and by choosing cell markers from the Human Proteome Atlas [[51]12,[52]14]. A brief workflow of the scRNA-seq pipeline is provided in [53]Fig. 1. Fig. 1. [54]Fig. 1 [55]Open in a new tab Outline of the workflow executed for the analysis of single-cell RNA sequencing (scRNA-seq) data. ENA- European Nucleotide Archive; RefGen- Reference Genome; DEGs- Differentially Expressed Genes; DEMGs- Differentially Expressed Marker Genes; FerrDB- Ferroptosis Database; GSEA- Gene Set Enrichment Analysis; KEGG- Kyoto Encyclopedia of Genes and Genomes; PPIs-Protein-Protein Interactions. 2.3. Identification of differentially expressed genes (DEGs), differentially expressed marker genes (DEMGs), ferroptosis-associated DEGs (FA-DEGs), and ferroptosis-associated DEMGs (FA-DEMGs) The screening of differentially expressed genes (DEGs) from the two abovementioned datasets was conducted by the “Cell Ranger counts” function present in Cell Ranger tools, and the identification was based on the Model-based Analysis of Single-cell Transcriptomics (MAST) algorithm. Differentially expressed genes (DEMGs) were screened by the Seurat package “R”, and the markers were identified by the Wilcoxon rank sum test. The test utilizes the “FindAllMarkers” function to determine the DEMGs between different clusters of cells, and the identification of DEMGs for each group of cells was made based upon a certain threshold of parameters (p_val_adj≤ 0.05, min.pct = 0.25 and |log2 Fold Change|≥ 1). Ferroptosis-associated genes were retrieved from the FerrDb database ([56]www.zhounan.org/ferrdb/) [[57]23]. To identify FA-DEGs (ferroptosis-associated differentially expressed genes) and FA-DEMGs (ferroptosis-associated differentially expressed marker genes), the intersection between the DEG datasets and the DEMGs bearing ferroptosis-associated genes was determined through the “InteractiVenn” tool ([58]http://www.interactivenn.net), and the results were visualized via a Venn diagram [[59]6,[60]15]. 2.4. Enrichment analysis of ferroptosis-associated DEMGs In an attempt to understand the biological process, molecular function, and cellular distribution of the FA-DEMGs, an enrichment analysis was performed with the help of the Gene Ontology (GO) database. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was used to determine the biological utility, molecular-level information, and cellular function of the genes in a biological system. In this study, we used ShinyGo V0.77 ([61]http://bioinformatics.sdstate.edu/go/), a web server, to identify Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and performed gene set enrichment analysis (GSEA) of ferroptosis-associated DEMGs (FA-DEMGs). 2.5. Analysis of the protein‒protein interaction (PPI) network and screening of ferroptosis-associated hub genes (FAHGs) The Search Tools for Retrieval of Interacting Genes/Proteins (STRING) database (version 11.5) ([62]https://string-db.org/) (PPI enrichment adjusted P < 0.05, confidence score ≥0.400) was used to construct a protein‒protein interaction (PPI) network of ferroptosis-associated differentially expressed marker genes (FA-DEMGs). Cytoscape software (version 3.9.1) ([63]https://cytoscape.org/) was used to construct and visualize the network of PPIs, and the CytoHubba plug-in was subsequently used to determine the hub genes and bottleneck genes based on the degree and bottleneck parameters [[64]16]. The top 10 genes were screened as candidate FAHGs. 2.6. Identification of candidate gene expression in the retinal atlas The expression of candidate genes was identified in different retinal cell types from two mouse datasets (E-MTAB-9061 and PRJNA653629) using different R packages (the Seurat package and ggplot2). 2.7. Docking analysis of targeted compounds The protein crystal structures of the candidate genes Prdx1 (PDB: [65]1QQ2) and Ubc (PDB: [66]3AUL) were retrieved from the RCSB Protein Data Bank database [67]https://www.rcsb.org/. The ligand libraries of FDA-approved anti-inflammatory drugs was downloaded from ChEMBL database [68]https://www.ebi.ac.uk/chembl/which contain 110 drugs respectively, The 3D chemical structure of 110 FDA-approved anti-inflammatory drug was retrieved from the ZINC15 database [69]https://zinc15.docking.org/substances/home/. Open Babel 3.1.1 software was utilized to convert the 3D chemical structure file format into PDB format and hydrogeneted by AutoDock Tools 4.2.6 software. Using PyMol 3.1 software, all the water molecules and heteroatoms in the PDB of ligand and protein structures were deleted before importing AutoDockTools software for hydrogenation and other operations and saved in PDB file format for docking. Fast molecular docking computations were conducted using AutoDock Vina software. Exhaustiveness was set to 10 iterations for optimal binding pose prediction. For PDB: [70]1QQ2, The grid box was centered at the standard point (X = 58.258, Y = 43.912, Z = 32.173 Å) and its dimensions were set to X = 80 Å, Y = 80 Å, and Z = 80 Å, for PDB: [71]3AUL The grid box was positioned at the standard coordinates (X = 22.904, Y = 26.473, Z = 22.612 Å) and assigned dimensions of X = 80 Å, Y = 80 Å, and Z = 80 Å.and results of docking were visualized by PyMol 3.1 software. Screening of the best-hit compounds mainly relies on their binding pattern and binding energy. Among the top 5 hits, dexamethasone was revealed to be a best hit candidate for both targets. The best-fit 3D structure of the docking complex was evaluated based on the binding energy of the ligand. 2.8. Prediction of the ADMET properties of the compound The different physicochemical and pharmacokinetic properties, including oral absorption, blood‒brain barrier permeability, and drug likeness, of the candidate drug compounds were analysed using SwissADME ([72]http://www.swissadme.ch/) [[73]17]. The ProTox-II ([74]https://tox-new.charite.de/protox_II/) online web server was utilized to predict the toxicity of compounds, including LD50, organ toxicity, and toxicity endpoints [[75]18]. 2.9. Molecular dynamics simulation Molecular dynamic simulations of the two identified target proteins (PRDX1 (PDB ID: [76]1QQ2) and UBC (PDB ID: [77]3AUL)) with the FDA-approved drug dexamethasone (ZINC ID: ZINC3875332) were conducted using the CHARMM36 force field in the GROMACS 2022 package. The CGenFF ([78]https://cgenff.umaryland.edu/initguess/) online server was used to generate the topology of the compound and force field parameter [[79]19]. To simulate the protein‒ligand complexes, the target proteins and docked protein‒ligand complexes were placed at the center of the cubic box, following which the periodic boundary condition (PBC) was set to 2 nm [[80]20]. A simple point charge (SPC) water model was used to recreate aqueous conditions around the box of protein‒ligand complexes [[81]21]. To maintain the electroneutrality of the system, water molecules were replaced with an equal number of Na^+ and Cl^− counterions. Furthermore, the steepest descent algorithm was employed for energy minimization. The ‘final production run’ of the equilibrated system was carried out for 100 ns. Xmgrace was used for analysing all the generated results [[82]22]. 3. Results 3.1. Identification of retinal cell types in two mouse model DR datasets Single-cell analysis was performed on openly available DR datasets from two different mouse models. In the E-MTAB-9061 dataset, a total of 7258 cells (control: 2009 cells and diabetic: 5249 cells) were examined following quality control and preprocessing of the data. A total of 100 % valid unique molecular identifiers (UMIs) from each cell were identified, and a total of 13,489 differentially expressed genes (DEGs) were detected. The cells were classified into 21 different transcriptionally distinct clusters based on specific parameters (1:20 dimension and 0.4 resolution) ([83]Fig. 2A, B, and 2C; [84]Fig. S1A). In this process, 7 different retinal cell types were identified and annotated accordingly. Thus, glial cells were represented by clusters 0, 5, 10, 12, 14, and 20; cone cells were represented by cluster 16; rod cells were represented by clusters 1, 3, 4, and 13; bipolar cells were represented by clusters 2, 6, and 15; fibroblasts were represented by clusters 7 and 8; pericyte cells were represented by cluster 9; endothelial cells were represented by clusters 8, 11, and 17; and retinal pigment epithelial cells were represented by clusters 18 and 19. While retinal ganglion cells and astrocytes were not identified, amacrine cells were present in relatively lower numbers in the total population of retinal cells. Fig. 2. [85]Fig. 2 [86]Open in a new tab Retinal atlas of scRNA-Seq analysis from E-MTAB-9061 Dataset. (A) tSNE plot showing 7258 isolated cells and 21 clusters that were identified from E-MTAB-9061 Dataset. (B) tSNE plot showing 7258 isolated cells and cell clusters grouped together as per each cell type. (C) tSNE plot of 7258 cells isolated from E-MTAB-9061 dataset. The color coding is done on the basis of sample source (control/diabetic). (D) Bar plot represents the percentage of affected cells among the different retinal cell types under control vs diabetic conditions. tSNE - t-Distributed Stochastic Neighbor Embedding; ctrl-control. Similarly, in the case of the PRJNA653629 dataset, a total of 133,036 cells (control: 76,139 cells and diabetic: 56,897 cells) were identified with valid unique molecular identifiers following quality control and data preprocessing. Here, a total of 17,836 DEGs were obtained, and the cells were categorized into 26 different clusters based on specific parameters (1:20 dimension and 0.5 resolution) ([87]Fig. 3A, B and 3C; [88]Fig. S1B). Seven retinal cell types were identified and annotated based on a cluster's well-derived marker expression. Clusters 4, 6, 17, 20, and 24 represented glial cells; Clusters 0, 1, and 2 represented cone cells; Clusters 0, 1, and 2 represented rod cells; Clusters 7, 8, 9, 10, 11, 13, 15, 16, 22, and 25 represented bipolar cells; Cluster 19 represented endothelial cells and pericytes; Clusters 12 and 14 represented amacrine cells; and Clusters 18 and 21 represented horizontal cells. Fig. 3. [89]Fig. 3 [90]Open in a new tab Retinal atlas of scRNA-Seq analysis from PRJNA653629 Dataset. (A) tSNE plot showing 133,036 isolated cells and 26 clusters that were identified from PRJNA653629 Dataset. (B) tSNE plot showing 133,036 isolated cells and cell clusters grouped together as per cell type. (C) tSNE plot of 133,036 cells isolated from PRJNA653629 dataset. The color coding is done on the basis of sample source (control/diabetic). (D) Bar plot represents the percentage of affected cells among the different retinal cell types under control vs diabetic conditions. tSNE-t-Distributed Stochastic Neighbor Embedding; ctrl-control. Differences in the cell populations of different retinal cell types in the E-MTAB-9061 and PRJNA653629 datasets were also identified to determine which cell types are primarily affected under diabetic conditions ([91]Fig. 2, [92]Fig. 3D). 3.2. Identification of DEGs, DEMGs, ferroptosis-associated DEGs and DEMGs DR-associated DEGs and DEMGs, along with upregulated and downregulated genes, were identified in both datasets (E-MTAB-9061 and PRJNA653629). Analysis of the E-MTAB-9061 dataset revealed 11,491 DEGs, of which 6362 were upregulated, while 5129 were downregulated (threshold parameters: p_val_adj≤ 0.05, min.pct = 0.25 and |log2-fold fold change|≥ 1). A total of 1946 DEMGs were also identified, of which 1428 remained upregulated and 518 were downregulated (threshold parameters: p_val_adj≤ 0.05, min.pct = 0.25 and |log2 Fold Change|≥ 1). Similarly, 13,332 DEGs were detected in the PRJNA653629 dataset (threshold parameters: p_val_adj≤ 0.05, min.pct = 0.25 and |log2-fold fold change|≥ 1), of which 3799 were upregulated and 9533 were downregulated. A total of 2109 DEMGs were identified; among them, 1860 were upregulated, and 249 were downregulated (p_val_adj≤ 0.05, min.pct = 0.25 and |log2 Fold Change|≥ 1) ([93]Table 2). Table 2. Summary of Differentially Expressed genes (DEGs) and Differentially Expressed Markers genes (DEMGs) in each datasets. Accession ID Up regulated Down regulated Total DEGs Up regulated Markers Down regulated Markers Total DEMGs E-MTAB-9061 6362 5129 11,491 1428 518 1946 PRJNA653629 3799 9533 13,332 1860 249 2109 [94]Open in a new tab ∗∗Threshold parameter: p_val_adj≤ 0.05, min.pct = 0.25 and |log2 Fold Change|≥ 1, Significant DEGS and DEMGs. After the identification of DEGs and DEMGs, 259 ferroptosis-associated genes were retrieved from the FerrDB database ([95]www.zhounan.org/ferrdb/) [[96]23]. We identified upregulated and downregulated DEGs and DEMGs from both datasets using a Venn diagram. The intersection between the 259 ferroptosis-associated genes and the DEGs from both datasets revealed 152 FA-DEGs. In addition, the intersection of the 259 ferroptosis-associated genes with the DEMGs in the two datasets revealed 63 FA-DEMGs ([97]Fig. 4A–D) that were considered for further analysis. Fig. 4. [98]Fig. 4 [99]Open in a new tab Venn diagram of Differentially Expressed Genes (DEGs). (A) Venn diagram represents the overlap between Differentially Expressed Genes (DEGs) from E-MTAB-9061 and PRJNA653629. (B) Venn diagram represents the overlap between Differentially Expressed Genes (DEGs) from all the three datasets (E-MTAB-9061, PRJNA653629 and FerrDB); the region of overlap represents 152 ferroptosis associated DEGs out of which 95 are up-regulated and 57 are down-regulated. (C) Venn diagram represents the overlap between Differentially Expressed Marker Genes (DEMGs) from E-MTAB-9061 and PRJNA653629. (D) Venn diagram represents the overlap between Differentially Expressed Marker Genes (DEMGs) from all three datasets (E-MTAB-9061, PRJNA653629 and FerrDB); the region of overlap represents 63 ferroptosis associated DEMGs out of which 52 are up-regulated and 11 are down-regulated. 3.3. Enrichment analysis of ferroptosis-associated DEMGs Gene Ontology and pathway enrichment analysis predicted the enrichment of FA-DEMGs in different biological processes, including chemical stress and oxidative stress-mediated cellular responses. GO enrichment analysis also revealed that the spatial distribution of the genes was mostly towards the apical and basal parts of the cell and its associated plasma membrane. The enrichment analysis also indicated that FA-DEMGs were probably associated with the execution of different molecular functions, such as protein kinase binding, peroxiredoxin activity, protein serine/threonine kinase inhibitor activity, oxidoreductase activity, and antioxidant activity. KEGG pathway analysis predicted the significant involvement of genes involved in different biochemical pathways, such as ferroptosis, mitophagy, autophagy, HIF1 signalling, mTOR signalling, MAPK signalling, PPAR signalling, NOD-like receptor signalling, and AGE-RAGE signalling, in diabetic complications and pathways related to neurodegeneration ([100]Fig. 5A–D). Gene Set Enrichment Analysis of FA-DEMGs also revealed that many of these genes were involved in various molecular, cellular, and biological processes that can be associated with the pathogenesis of DR. Detailed information on the top 15 pathways is shown in ([101]Table 3). Fig. 5. [102]Fig. 5 [103]Open in a new tab Gene ontology and Gene Set Enrichment Analysis predicted enrichment of Ferroptosis-associated DEMGs. (A) Dot plot represents enrichment of biological process of FA-DEMGs. (B) Dot plot represents enrichment of cellular components of FA-DEMGs. (C) Dot plot represents enrichment of molecular function of FA-DEMGs. (D) Bar plot represents enriched KEGG pathways of FA-DEMGs. Color coding of all the above figures has been done based on the extent of false discovery rate (-1og10 FDR). BP- Biological Process; CC- Cellular Component; MF- Molecular Function. Table 3. Detailed information of the top 15 pathways. ID Fold Enrichment Pathway mmu04216 78.16071429 Ferroptosis mmu04137 45.97689076 Mitophagy mmu04140 22.33163265 Autophagy mmu03320 19.51578384 PPAR signaling pathway mmu04066 18.44500632 HIF-1 signaling pathway mmu04933 13.75766148 AGE-RAGE signaling pathway in diabetic complications mmu04370 11.97865353 VEGF signaling pathway mmu05167 11.31007752 Kaposi sarcoma-associated herpesvirus infection mmu04621 10.31824611 NOD-like receptor signaling pathway mmu04217 9.982211275 Necroptosis mmu04150 8.907203907 MTOR signaling pathway mmu04010 8.270975057 MAPK signaling pathway mmu05012 6.629407488 Parkinson disease mmu05022 5.184790334 Pathways of neurodegeneration mmu05200 5.127394131 Pathways in cancer [104]Open in a new tab 3.4. Network construction and identification of FAHGs Following the enrichment analysis, a protein‒protein interaction network (PPI) was constructed for the FA-DEMGs. The network was found to possess 63 nodes and 207 interactions with confidence (0.400). The top ten Hub genes identified were Hif1a, Jun, Ubc, Atf3, Prdx1, Epas1, Hspb1, Bnip3, Slc2a1 and Dusp1, whereas the bottleneck genes were Hif1a, Prdx1, Ubc, Jun and Fth1. The screening of the Hub genes and Bottleneck genes from the PPI network was conducted based on high degree and bottleneck values, respectively ([105]Fig. 6A–D) ([106]Table 4). To identify potential candidates for ferroptosis-associated pathogenesis in DR, genes common to both the ‘Hub gene’ and ‘Bottleneck gene’ (i.e., Prdx1, Ubc, Hif1a, and Jun) were considered for further studies. However, Prdx1 and Ubc are of primary interest because of their novelty in the context of DR-associated ferroptosis. Fig. 6. [107]Fig. 6 [108]Open in a new tab Prediction of Interaction network and identification of Hub genes and Bottleneck genes from FA-DEMGs. (A) Interaction network of FA-DEMGs obtained from STRING Database represents 63 nodes and 177 interactions. (B) Interaction network of FA-DEMGs obtained from Cytoscape represents 56 nodes and 354 interactions. The color coding is based on the extent of protein-protein interactions (purple indicates the highest degree of interactions followed by light purple, grey, cyan blue, green, light green and yellow). (C) Protein-protein interactions represent the top 5 hub genes identified from Cytoscape network. The color coding is based on the degree of interaction (Red represents the highest interaction). (D) Protein-protein interactions represent the top 5 bottleneck genes identified from Cytoscape network. The color coding is based on bottleneck value (Red represents the highest value). Table 4. Tropological Parameters of Candidate genes (Hub genes and Bottleneck genes) found in PPIs network. Types of genes Genes Degree Hub Genes Hif1a 50 Jun 40 Atf3 30 Ubc 26 Epas1 26 Bnip3 22 Hspb1 20 Prdx1 16 Slc2a1 14 Dusp1 14 Bottleneck Genes Hif1a 9 Prdx1 9 Ubc 8 Jun 6 Fth1 6 [109]Open in a new tab 3.5. Expression of candidate genes associated with retinal heterogeneity Candidate gene expression was evaluated in the context of retinal heterogeneity by creating a violin plot. While the Prdx1 and Ubc genes exhibited high expression in all retinal cell types ([110]Fig. 7A–D), Hif1 alpha and Jun showed variable expression ([111]Figs. S2A and S2B). Hence, the unvarying expression of Prdx1 and Ubc among all retinal cell types provided a suitable reason to further select them as candidate genes associated with ferroptosis in DR. Fig. 7. [112]Fig. 7 [113]Open in a new tab Violin plot representing the expression of candidate genes among the different retinal cell types. (A) Violin plot represents Prdx1 expression among the different retinal cell types obtained from E-MTAB-9061 dataset. (B) Violin plot represents Ubc expression among the different retinal cell types obtained from E-MTAB-9061 dataset. (C) Violin plot represents Prdx1 expression in retinal cell types of PRJNA653629 dataset. (D) Violin plot represents Ubc expression in retinal cell types of PRJNA653629 dataset. The retinal cell types are color coded in the graph. 3.6. Screening of potential drug compounds and molecular docking Following the selection of Prdx1 and Ubc as potent candidate genes for DR-associated ferroptosis, screening of FDA-approved drugs were conducted, and dexamethasone was identified to be a candidate drug compound. The search for drug compounds was carried out from ChEMBL database ([114]https://www.ebi.ac.uk/chembl/), after which the 3-D structures of 110 drug compounds were retrieved from the ZINC15 database [115]https://zinc15.docking.org/substances/home/. The probable binding affinities of the drugs for the protein structures of PDB: [116]1QQ2 and PDB: [117]3AUL were investigated by fast molecular docking computations. Top 5 hits were listed on ([118]Table 5, [119]Table 6) based on the binding affinity with the target proteins. Among the top 5 hits, dexamethasone exhibited increased binding affinity for both PDB: 1QQ2and PDB: [120]3AUL. The following docking results revealed the binding affinity, hydrogen bonding, intermolecular energy (Van der Waals forces and Hydrogen bonding), s and electrostatic interaction between the amino acid residue of 1QQ2 and dexamethasone are ΔG = −8.7 kcal/mol, −7.55 kcal/mol, −0.3 kcal/mol Chain A: VAL149 ARG151 SER152 ASP154 GLU155 ARG158 GLY161. Chain B: VAL149 SER152 ASP154 GLU155 ARG158 PRO174respectively. Bond Distance between the amino acid residue of 1QQ2 and dexamethasone are ARG151–1.979, ARG158–1.894, GLY161–2.031).In the case of 3AUL, the binding affinity and hydrogen bonding, intermolecular energy (Van der Waals forces and Hydrogen bonding), and electrostatic interaction with the respective amino acid residues were as follows: dexamethasone are ΔG = −8.0 kcal/mol, −7.16 kcal/mol, −0.13 kcal/mol, chain A: THR22 GLU24 ASN25 LYS27 ALA28 GLN31 ASP32 PRO38 ASP39 LEU50 GLU51 ASP52 LYS27 LYS29 LYS33 chain B: LYS6 LEU8 ILE23 GLU24 LYS27 LYS29 LYS33 PRO38 ASP39 GLN41 ARG42 ILE44 GLY47 LYS48 GLN49 LEU50 GLU51 ASP52 GLY53 HIS68 LEU69 VAL70 ARG72 respectively. Bond Distance between the amino acid residue of 3AUL and dexamethasone are LYS29–1.762, LYS33–1.572. ([121]Fig. 8, [122]Fig. 9, [123]Table 7). Table 5. List of Top 5 ligands that have low binding energy with Prdx1. ZINC ID/DRUG NAME CHEMICAL STRUCTURE MOLECULAR WEIGHT (g/mol) BINDING ENERGY (kcal/mol) ZINC3875332/DEXAMETHASONE Image 1 392.5 −8.7 ZINC4212809/DEFLAZACORT Image 2 441.5 −8.2 ZINC3977978/FLUOCINONIDE Image 3 494.5 −8.0 ZINC4212851/DESONIDE Image 4 416.5 −7.9 ZINC4213357/FLUROMETHOLONE Image 5 376.5 −7.8 [124]Open in a new tab Table 6. List of Top 5 ligands that have low binding energy with Ubc. ZINC ID/DRUG NAME CHEMICAL STRUCTURE MOLECULAR WEIGHT (g/mol) BINDING ENERGY (kcal/mol) ZINC3875332/DEXAMETHASONE Image 6 392.5 −8.0 ZINC4212809/DEFLAZACORT Image 7 441.5 −7.6 ZINC3977978/FLUOCINONIDE Image 8 494.5 −7.4 ZINC4212851/DESONIDE Image 9 416.5 −7.3 ZINC4213357/FLUROMETHOLONE Image 10 376.5 −7.2 [125]Open in a new tab Fig. 8. [126]Fig. 8 [127]Open in a new tab Molecular docking model represents the 3-dimensional physical association between the drug compound (A) Dexamethasone and Prdx1 (B) Dexamethasone and Ubc. (Blue dotted line: Hydrogen bonding, Grey dotted line: Hydrophobic interaction, Yellow dotted line: Electrostatic interaction). Fig. 9. [128]Fig. 9 [129]Open in a new tab Molecular docking model represents in 3-dimensional docking interaction showing bond distance between amino acid residue and ligand (A) Dexamethasone with 1QQ2 (B) Dexamethasone with 3AUL. Table 7. List of ligands that have low binding energy, Intermolecular energy and Electrostatic energy with proteins. ZINC ID/DRUG NAME CHEMICAL STRUCTURE AutoDock PBE (kcal/mol)/(Target protein) Residue Bond Distance (Å) Intermolecular Energy (VdW + Hbond) Residue/(kcal/mol) Hydrophobic residue Electrostatic Energy Residue/(kcal/mol) ZINC3875332/DEXAMETHASONE Image 11 −8.7 kcal/mol (Prdx1) ARG151 ARG158 GLY161 1.979 1.894 2.031 ARG151 ARG158 GLY161 −7.55 ALA58 VAL153 GLU155 0.3 ASP154 Image 12 −8.0 kcal/mol (Ubc) LYS29 LYS33 1.762 1.572 LYS29 LYS33- 7.16 ASN25 ALA28 ILE23 LEU50 LYS27 ARG72 0.13 ASP52 [130]Open in a new tab 3.7. ADMET properties of the compound At an early stage of drug development, physiologically undesirable drug candidates can be identified in a significantly less time-consuming and more cost-effective manner with the help of molecular docking-based virtual screening. The pharmacokinetic properties of the ligands were predicted with the help of in silico tools. This study evaluated the pharmacokinetic, drug likeness, and toxicity of dexamethasone. The pharmacokinetic properties of chemical compounds indicated that GI tract absorption is high for dexamethasone. The intestinal absorption of dexamethasone was >85 %, while that of dexamethasone was impermeable to the blood‒brain barrier. Dexamethasone did not inhibit cytochrome P450 enzymes such as 1A2, 2C19, 2D6, 2C9, or 3A4. Dexamethasone was found to be a good skin permeant. The selected drugs followed Lipinski's rule of five without any violations, and the bioavailability score was 0.55. The ADME results of the chemical compounds are listed in ([131]Table 8). Table 8. ADME evaluation of chemical compounds. Drug Pharmacokinetics (ADME) Drug likeness Compound GI absorption BBB permeant Pgp substrate CYP1A2 inhibitor CYP2C19 inhibitor CYP2C9 inhibitor CYP2D6 inhibitor CYP3A4 inhibitor log Kp(cm/s) Lipinski violations Bioavailability Dexamethasone High No Yes No No No No No −7.32 No 0.55 [132]Open in a new tab The toxicity evaluation of the drug compounds revealed that dexamethasone exhibited active immunotoxicity. The evaluation of the toxicity of these chemical compounds is shown in ([133]Table 9). Table 9. Toxicity evaluation of chemical compounds. Compound Toxi city LD[50] Hepato toxicity Carcino Genicity Immuno toxicity Muta genicity Cyto toxicity Dexamethasone Class 5 3000 mg/kg 0.99 0.72 0.99 0.69 0.72 (inactive) (inactive) (active) (inactive) (inactive) [134]Open in a new tab 3.8. Molecular dynamics simulation A molecular dynamics simulation study was conducted to better understand the protein‒ligand complex's binding mechanism and structural stability. Using GROMACS software, target proteins (PRDX1 (PDB ID: [135]1QQ2) & UBC (PDB ID: [136]3AUL)) and docked complexes (PRDX1 with dexamethasone & UBC with dexamethasone) were simulated for more than 100 ns. Conventional analytical techniques were used to generate MD simulation trajectories. The simulation trajectories were analysed by using the root mean square deviation (RMSD) and root mean square fluctuation (RMSF), which revealed the deviation and fluctuation of proteins (PRDX1 and UBC) and protein‒ligand complexes (PRDX1 with dexamethasone and UBC with dexamethasone). To analyse the simulation trajectory, the radius of gyration (Rg) obtained from the simulation was used to evaluate the structural stability and compactness of proteins (PRDX1 and UBC) and their ligand-bound complexes (PRDX1 with dexamethasone and UBC with dexamethasone). The Rg values, which ranged between 1.58 nm and 1.65 nm for PRDX1 and between 1.17 nm and 1.22 nm for UBC, remained steady throughout the simulation for the proteins as well as their docked complexes with dexamethasone. The average Rg of unbound PRDX1 was 1.58 nm, and that of UBC was 1.20 nm. The average Rg values of PRDX1 and UBC were found to be similar to those of their complex forms with dexamethasone. Overall, the Rg of the PRDX1-dexamethasone complex was 1.65 nm, whereas the Rg of the UBC-dexamethasone complex was lower, at 1.22 nm ([137]Fig. 10A). The RMSD fluctuations of unbound PRDX1 and UBC were greater than those of their complex forms. Thus, the RMSD plot revealed that the curve pertaining to the PRDX1-dexamethasone complex fluctuated at intervals of 20–30 ns, 50–60 ns and 80–90 ns, after which it finally flattened. Fig. 10. [138]Fig. 10 [139]Open in a new tab Interaction analysis graphs of unbound PRDX1 (PDB ID: [140]1QQ2) and UBC (PDB ID: [141]3AUL) along with their docked complex with Dexamethasone following Molecular Dynamics Simulations. (A) Radius of gyration vs Time (B) RMSD vs Time (C) RMS Fluctuation vs Atom (D) RMS Fluctuation vs Residue. RMSD- Root Mean Square Deviation; RMSF- Root Mean Square Fluctuation; Unbound Protein and Protein-ligand complex are color coded in the plot. On the other hand, the UBC-dexamethasone complex exhibited fluctuations at 10 ns but remained flat throughout the run ([142]Fig. 10B). The RMSD values of the PRDX1-and UBC- complexes with dexamethasone were found to be 0.1–0.5 nm and 0.1–0.3 nm, respectively. The RMSD values remained more stable throughout the simulation for both complexes. These data indicated that dexamethasone might be a potentially effective inhibitor of selected ferroptosis-associated protein targets due to the formation of stable complexes. The RMS fluctuation (RMSF) values revealed the protein dynamic equilibrium and mobility rate of all the amino acid residues. Higher RMSF values denoted unstable regions, whereas protein-stable regions showed lower RMSF values. During simulation, the flexible residues of ligand binding areas were transformed simultaneously based on ligand recognition to maintain equilibrium. The RMSF values of the PRDX1-dexamethasone and UBC-dexamethasone complexes remained unaltered throughout the simulation. At particular intervals during simulation, the RMSFs of the UBC-dexamethasone complex and the PRDX1-dexamethasone complex remained constant ([143]Fig. 10C and D). The selected ligands consistently displayed the best binding affinity with the target proteins. Therefore, PRDX1 and UBC could be potent targets for dexamethasone. 4. Discussion Diabetic retinopathy (DR) is a chronic neurovascular complication associated with diabetes that primarily affects the eyesight of working-age people and is one of the major causes of visual impairment worldwide [[144]24]. Reports predict that the global incidence of DR will increase at an alarming rate by 2045, when the number of DR cases will reach approximately 160.45 million [[145]2]. NPDR is the initial phase of DR progression characterized by weakened retinal blood vessels, leading to fluid leakage and macular swelling, eventually leading to vision blurriness. PDR is a more severe phase of DR and is characterized by unusual neovascularization accompanied by leakage of blood into the back of the eye, ultimately leading to permanent vision loss [[146]25]. Therefore, understanding the pathogenesis of DR is critical for alleviating the risk of blindness in diabetic patients [[147]26]. Ferroptosis is an iron-induced apoptotic event that occurs because of the accumulation of lipid peroxides in the cell membrane and is very different from conventional cell death events such as autophagy, necrosis, and apoptosis [[148]27]. Recent studies have indicated the importance of ferroptosis in the pathogenesis and development of DR [[149]28]. In vivo and in vitro studies revealed that ferroptosis is a new mode of iron-induced programmed cell death associated with metabolic pathways and the pathobiological conditions of several eye disorders [[150]29]. To date, various hallmarks of ferroptosis have been characterized; the key cellular features are a reduction in mitochondrial volume and cristae without disruption of the integrity of the cell membrane or nucleus. Although the depletion of intracellular glutathione (GSH) and the reduction of glutathione peroxidase-4 (GPX4) are considered the main causes of ferroptosis, much remains unknown regarding the mechanism of DR-associated ferroptosis. This study is essential because it helps to better understand ferroptosis in the context of DR. In this study, in silico methods were used to analyse single-cell datasets to discover candidate genes and inhibitors to target ferroptosis-associated pathogenesis in DR. We identified 63 FA-DEMGs from all the datasets. Downstream analysis of the FA-DEMGs revealed associations between genes and pathways involved in ferroptosis, oxidative stress, reactive oxygen species metabolism, HIF1 signalling, and MTOR signalling, all of which have been reported to be related to DR pathogenesis [[151]26,[152][30], [153][31], [154][32], [155][33], [156][34], [157][35]]. The results of the enrichment analysis further validated the pathways associated with the previously identified FA-DEMGs identified in our study. The PPI network revealed the hub genes and bottleneck genes related to DR-mediated ferroptosis. The genes common to the ‘Hub gene’ and ‘bottleneck gene’ were Hif1a, Jun, Ubc, and Prdx1. However, we focused on Prdx1 and Ubc, primarily because these genes are novel because no studies have reported their association with DR-related ferroptosis. In addition, these genes were also predicted to be highly expressed among all retinal cell types. In contrast, other genes, such as Hif1a and Jun, are well-known biomarkers of diabetes-associated ferroptosis [[158]5,[159]36]. Hif1a is generally known as the key member that orchestrates the platform necessary for cells to adapt to a hypoxic environment [[160]37], whereas the transcription factor Jun has been found to be involved in mitophagy-related pathways [[161]38]. Our results support previous findings showing that the genes of interest, i.e., Prdx1 and Ubc, remained upregulated in DR [[162][39], [163][40], [164][41], [165][42], [166][43]]. Prdx1 belongs to the selenium-independent glutathione peroxidase family, a group of antioxidant enzymes that are involved in the reduction of hydrogen peroxide levels and aid in ferroptosis suppression. In addition to being an antioxidant, Prdx1 can function as an ROS-scavenging enzyme, chaperone, tumor suppressor protein, and signal transducer [[167]44]. Prdx1 remains distributed in the cytosol and associates with the lipid domain of the plasma membrane, which helps regulate growth factor signalling receptors [[168]45]. Prdx1 protects telomeres from oxidative damage, regulates lipid peroxidation in corneal endothelial cells, and plays a significant role in tumorigenesis and cancer progression by regulating mTOR signalling pathways and inflammatory cytokine production [[169]29,[170][46], [171][47], [172][48]]. Ubc plays a significant role in cellular processes such as heat shock, translational impairment, and oxidative stress [[173]49]. Ubc acts as a ubiquitin precursor of the ubiquitin‒proteasome system and thereby helps maintain steady levels of the protein. The ubiquitin‒proteasome system is crucial for inflammation, autophagy, and cancer progression [[174]50]. Our work showed that high expression of Prdx1 and Ubc in different retinal cell types might be linked to ferroptosis-associated DR. However, their exact role in DR progression and pathogenesis remains unclear. We further validated the role of the drug compounds Prdx1 and Ubc as therapeutic targets for DR-mediated ferroptosis. Docking analysis of these critical genes with the drug compound dexamethasone predicted high binding affinity, confirming the role of these two genes as proper therapeutic targets of DR. Molecular dynamics simulation revealed that dexamethasone had a stable effect on both PRDX1 and UBC. Dexamethasone is an FDA-approved corticosteroid drug known for its anti-inflammatory properties. It can manifest active immunotoxicity in certain conditions, complicating its use due to its mechanism and immunosuppressive properties [[175]51,[176]52]. It can contribute to its immunotoxic effects when used at high doses or for prolonged periods. However, a lot of clinical research and trials have corroborate the therapeutic potential of dexamethasone in various diseases [[177][53], [178][54], [179][55], [180][56], [181][57], [182][58]]. Dexamethasone as a therapeutic drug is known for its well-documented efficacy in managing inflammatory and hyperimmune conditions.It is often applied topically or "intravitreally to alleviate postoperative inflammation in the eyes and provide relief in patients with dry-eye disease. Recent studies have emphasized the application of dexamethasone for treating various eye disorders, such as diabetic and nondiabetic macular edema and uveitis [[183]59]. This can be established through careful assessment, proper dosing, and rigorous observation, which are crucial for ensuring safe and effective use of this medication and its recommendation as a therapeutic drug. Although dexamethasone has proven to be a promising option for treating retinal disorders, our study indicates that dexamethasone is suitable for treating DR-associated ferroptosis and deserves in-depth research context in the coming years. Our study is important because it reports for the first time the predicted role of several novel genes responsible for ferroptosis in the of DR. However, our research has several limitations. First, our results are based on in silico analysis of openly available datasets and data mining and require validation by wet laboratory experiments. Moreover, the biological and molecular functions of the candidate genes and the therapeutic activity of the compounds need to be confirmed by in vitro and in vivo experiments. Future studies should focus on validating reported genes and the therapeutic activity of a suitable drug compound using in vitro and in vivo DR models. 5. Conclusions The present study reports that ferroptosis-associated marker genes (Prdx1 and Ubc) are associated with the development and pathogenesis of DR by employing in-depth in silico analysis. These genes are predicted to be critical biomarkers of ferroptosis and therapeutic targets for DR. Molecular docking analysis revealed a high binding affinity between drug compounds and their protein targets. Molecular dynamics simulations were performed to calculate the dynamics and stability of protein‒ligand complexes (dexamethasone with PRDX1 and dexamethasone with UBC), which might prove to be helpful in the development of DR treatments in the future. Overall, our research provides a novel understanding of the pathogenesis of ferroptosis, along with a theoretical basis for investigating novel markers and therapeutic approaches for DR. CRediT authorship contribution statement Md. Maqsood Ahamad Khan: Formal analysis, Conceptualization. Ananya Ganguly: Writing – review & editing, Methodology. Shubhrajit Barman: Writing – review & editing, Validation, Methodology. Chirasmita Das: Formal analysis. Senthil Kumar Ganesan: Writing – review & editing, Supervision, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization. Ethics approval and consent to participate Not applicable. Consent for publication All authors read and approved the paper for publication. Availability of data and materials The datasets that contributed to the outcome of this study are publicly available in ArrayExpress ([184]https://www.ebi.ac.uk/biostudies/arrayexpress) and European Nucleotide Archive (ENA) ([185]https://www.ebi.ac.uk/ena/browser/home), Reference number (E-MTAB-9061 and PRJNA653629). Funding This work was supported by the Indian Council of Medical Research (ICMR), India (Grant no: 5/4/6/8/OPH/2020-NCD-II and Grant no: 5/4/6/30/OPH/2021-NCD-II) and institutional Grant no P07-MLP-120 from CSIR‐ Indian Institute of Chemical Biology, Kolkata, India. Declaration of competing interest The authors declare that there is no conflict of interest. Acknowledgments