Introduction

Diabetic retinopathy (DR) is a common chronic complication of diabetes mellitus and a leading cause of blindness in adults.1 In 2020, an estimated 103 million individuals worldwide were affected by retinopathy, with projections suggesting an increase to 130 million by 2030.2 The pathogenesis of DR is thought to involve impaired blood-retinal barrier (BRB) integrity, leading to neurovascular unit damage. If left untreated, BRB dysfunction can result in severe visual impairment.3 DR is classified into non-proliferative and proliferative stages.4 Current treatment options include panretinal laser photocoagulation, surgery, blood glucose control, blood pressure management, serum sterol regulation, and intraocular interventions. The management of DR faces significant challenges due to the extensive patient base and the substantial costs of existing treatments, which limit accessibility and compromise outcome satisfaction. Consequently, a deeper exploration of DR pathogenesis is essential. This effort is critical for pinpointing new therapeutic targets and providing actionable clues for drug discovery, which could ultimately enhance treatment efficacy and availability.5–7

Forkhead box K1 (FOXK1) is an evolutionarily conserved transcription factor (TF) recently identified as a key regulator in various cancers. Members of the FOXK family play pivotal roles in several biological processes, including cell proliferation, differentiation, apoptosis, autophagy, cell cycle progression, DNA damage response, and tumorigenesis. Dysregulation of FOXK proteins can disrupt cellular fate, promoting tumor development and progression. Regulatory mechanisms of FOXKs include post-translational modifications (PTMs), microRNA (miRNA) interactions, and protein-protein interactions.8 FOXK1 has been implicated as a significant TF in autophagy.9 Similarly, autophagy is abnormally activated in DR, as evidenced by autophagosome accumulation. Mechanistic studies have shown that DR induces PTEN expression, inhibiting AKT/mTOR phosphorylation and triggering abnormal autophagy and apoptosis.10 It is noteworthy that inflammatory signaling pathways, such as the activation of NF-κB, can interact with both autophagy and apoptosis pathways. This suggests the existence of potential crosstalk between inflammation-related mechanisms and FOXK1-mediated regulation.11–13 Our previous research demonstrated that silencing FOXK1 overexpression in a high-glucose cultured human umbilical vein endothelial cell (HUVEC) model led to increased cell apoptosis, migration, and neovascularization via activation of the p-AKT/AKT signaling pathway. These findings collectively suggest a potential role for FOXK1 in the pathogenesis and development of DR. However, further investigation is required to fully elucidate the role of FOXK1 in DR.

Based on publicly available transcriptomic data of DR, this study innovatively focused on the transcription factor FOXK1 and employed Weighted Gene Co-expression Network Analysis (WGCNA) to identify module genes significantly associated with FOXK1. Building on this, we systematically conducted screening for DR biomarkers. Furthermore, we extensively investigated the diagnostic potential of these biomarkers and their underlying molecular regulatory mechanisms. This study employed an integrated bioinformatics and experimental approach to investigate the role of the FOXK1 transcription factor in DR. Building upon the initial finding of FOXK1 upregulation in DR patients, we not only validated its functional impact on endothelial cells but also identified and experimentally confirmed SPDEF and SLC25A41 as novel FOXK1-related diagnostic biomarkers. Furthermore, the constructed regulatory network and drug prediction, particularly the highlighted potential of resveratrol, provide a multi-dimensional perspective on the underlying mechanisms. Therefore, this work extends beyond preliminary evidence by offering specific biomarkers and mechanistic clues, which are expected to enhance the understanding of DR pathogenesis and suggest promising avenues for early diagnosis and targeted therapeutic strategies.

Materials and MethodsData Provenance

GSE221521 (GPL24676 platform) served as the training cohort, and GSE189005 (GPL23126 platform) as the validation cohort. Both datasets were retrieved from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). GSE221521 included blood samples from 69 patients with DR and 50 healthy controls, while GSE189005 consisted of whole blood cell samples from 10 patients with DR and 10 control individuals. The research workflow of this study was shown in Supplementary Figure 1.

Differentially Expressed Gene Identification

Differentially expressed genes (DEGs) between patients with DR and healthy controls in GSE221521 were analyzed using DESeq2 (v 1.38.0).14 Genes with |log2Fold Change (FC)| > 0.5 and p-value < 0.05 were considered differentially expressed. Visualization of DEGs was performed with ggplot2 (v 3.3.5)15 and ComplexHeatmap (v 2.16.0).16

Weighted Gene Co-Expression Network Analysis (WGCNA)

The WGCNA package (v 1.49.4)17 was used to identify gene modules associated with FOXK1. First, FOXK1 expression levels were compared between patients with DR and normal individuals in GSE221521. The GoodSamplesGenes function was used to cluster the samples from the training cohort into distinct groups. The PckSoftThreshold function determined the optimal soft threshold for constructing the scale-free network (R2 = 0.8881339). Hierarchical clustering was then performed using the Dynamic Tree Cutting method. Pearson’s correlation analysis was applied to calculate correlation coefficients and p-values between the modules and traits to identify key module genes significantly associated with FOXK1.

Pinpointed Candidate Genes and Functional Enrichment Evaluation

Candidate genes were derived by overlapping DEGs and key module genes using the VennDiagram package (v 1.7.3).18 To explore the potential pathways and biological functions of these candidate genes, the clusterProfiler package (v 4.8.2)19 was utilized to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses, focusing on those with p-values < 0.05. GO enrichment included cellular components (CC), biological processes (BP), and molecular functions (MF). In the analysis, the false discovery rate (FDR) method was used for adjustment, and a p-adjust < 0.05 was set as the screening threshold.

Identification and Validation of Biomarkers

For biomarker screening, the least absolute shrinkage and selection operator (LASSO) method was applied using the glmnet package (v 4.1–8).20 In this process, to find the optimal balance between model complexity (lambda value) and prediction accuracy, 5-fold cross-validation was employed, and the lambda value that minimized the cross-validation error (lambda.min) was ultimately selected. The expression of the candidate biomarkers was then validated in both GSE221521 and GSE189005. The Wilcoxon test was used to compare the expression levels of candidate biomarkers, and those exhibiting consistent and significant expression trends were selected as biomarkers for further analysis (p-value < 0.05). Results were visualized using the ggplot2 package (v 3.3.5).15

Construction of Nomogram

To assess the clinical diagnostic potential of biomarkers, a predictive nomogram was constructed using the rms package (v 6.5–0) (https://CRAN.R-project.org/package=rms). The nomogram’s predictive accuracy was evaluated through calibration plots and decision curve analysis. The calibration curve was generated using the replot package (v 1.1) (https://cran.r-project.org/web/packages/regplot/index.html), while decision curve analysis was performed with the rmda package (v 1.6) (https://cran.r-project.org/web/packages/rmda/index.html).

Analysis of Immune Infiltration

To investigate immune microenvironment changes in DR, the CIBERSORT method was employed to determine the expression abundance of 22 immune cell types between patients with DR and controls in the training set. In this process, immune cells with an average expression of 0 across all samples were filtered out. The Benjamini-Hochberg multiple testing adjustment with a p-adjust < 0.05 was used as the criterion for screening significant differences. The Wilcoxon test was applied to identify immune cells with significantly different expression levels between the disease and control groups (p-value < 0.05). Spearman correlation analysis, conducted with the psych package (v 2.2.9),21 was used to explore the relationship between biomarkers and differential immune cells, further elucidating the immune system alterations in patients with DR.

Gene Set Enrichment Analysis (GSEA)

To further explore the pathways and intrinsic biological mechanisms involving the biomarkers, C2: KEGG gene sets were downloaded as background sets using the msigdbr package (v 7.5.1) (https://cran.r-project.org/web//packages//msigdbr/index.html). The corrplot package (v 0.92)22 was used to compute Spearman correlation coefficients between biomarkers and a set of additional genes, followed by sorting them in descending order. GSEA was performed using the clusterProfiler package (v 4.2.2).23 The FDR method was adopted for adjustment, with p-adjust < 0.05 set as the screening threshold.

Molecular Regulatory Networks of Biomarkers

The molecular regulatory network plays a key role in gene expression regulation, offering valuable insights into its complexity and diversity. Based on the identified biomarkers, miRNAs regulating these biomarkers were predicted using the miRTarBase database (https://mirtarbase.cuhk.edu.cn/~miRTarBase/miRTarBase_2022/php/index.php) and the TarBase database (https://dianalab.e-ce.uth.gr/tarbasev9). miRNAs co- predicted in both databases were considered as key regulators. TFs associated with the biomarkers were identified using the JASPAR database (https://jaspar.elixir.no/).

To identify additional genes associated with biomarker functions, the GeneMANIA database (http://genemania.org) was used to predict genes related to the biomarkers’ functions and associated pathways. The Cytoscape software (v 3.10.0)24 was employed to visualize the molecular regulatory networks.

Biomarker-Based Drug Prediction

To identify potential therapeutic agents for biomarkers, the DSigDB database (https://dsigdb.tanlab.org/DSigDBv1.0/) was utilized to predict drugs targeting DR-associated biomarkers. The database integrates experimentally validated data from multiple public resources, including GEO, Connectivity Map, and Library of Integrated Network-based Cellular Signatures (LINCS), and contains drug-gene regulatory relationships confirmed through cell-based, animal, or preclinical studies, ensuring high reliability. The following strict criteria were applied during screening: 1) The regulatory relationship between a compound and the target gene must be supported by at least one independent experimental study; 2) The compound must be an approved drug, be in clinical trial stages, or be a natural product/small molecule compound with clear drug development potential; 3) Compounds with known toxicity to retinal tissue based on existing evidence were excluded. Through the above process, a series of candidate compounds was identified. To visually represent the regulatory relationships between biomarkers and corresponding drugs, an interaction network between biomarkers and drugs was constructed and visualized using Cytoscape software (v 3.10.0).

RNA Extraction and Reverse Transcription-Quantitative PCR (RT-qPCR)

To validate the biomarkers, clinical samples were collected for RT-qPCR validation between January and April 2025. The samples, obtained from Wuhan Aier Hanyang Eye Hospital, included 9 DR samples and 9 healthy controls. Ethical approval for the study was granted by the Ethics Committee of Wuhan Aier Eye Hanyang Hospital, with all patients providing written informed consent. No minors were included in the study. Total RNA was extracted using TRIzol reagent (Servicebio, China), and RNA purity was assessed with a NanoDrop One spectrometer (Thermo Fisher Scientific, MA, USA). RT-qPCR was performed on a LightCycler 96 Real-Time PCR System (Hoffmann-La Roche Ltd). All the primers used were synthesized by Wuhan Tsinghua Biotechnology Co., Ltd., with detailed information shown in Table 1. β-actin was used as an internal reference. Following each qPCR reaction, a melting curve analysis was performed, and samples whose Ct values were identified as abnormal were excluded from the final data analysis. The relative expression of the target genes was calculated using the formula: relative expression = 2−ΔΔCq. P-values < 0.05 were considered statistically significant.

Table 1 Primers Sequence Table

Statistical Analysis

Statistical analyses were performed using R software (v 4.2.1). Between-group comparisons were conducted using the Wilcoxon rank-sum test for non-normally distributed data and Student’s t-test for normally distributed data. For multi-group comparisons, one-way Analysis of Variance (ANOVA) or Kruskal–Wallis test was employed followed by post hoc tests with appropriate correction (eg, Tukey’s or Dunn’s test). Correlation analysis was performed using Spearman’s method for non-parametric data. The Benjamini–Hochberg method was applied for multiple testing correction in enrichment analyses, with a FDR < 0.05 considered significant. All statistical tests were two-tailed, and a p-value < 0.05 was considered statistically significant unless otherwise specified.

Results202 Genes Identified as Differentially Expressed Genes in DR

In GSE221521, a total of 202 DEGs were identified between DR and normal samples, with 155 up-regulated and 47 down-regulated genes (Supplementary Table 1). The top five up-regulated genes were KNG1, IFNK, TAS2R30, OR5B21, and AK7, while the top five down-regulated genes were AC004832.3, CCN2, PDE1C, LY6G6E, and AC233723.1, as shown in Figure 1A. The top ten up- and down-regulated DEGs were selected for heatmap visualization (Figure 1B).

Figure 1 Identification of differentially expressed genes (DEGs) in patients with diabetic retinopathy (DR) and control samples. (A) Volcano plot of DEGs: blue dots on the left represented downregulated genes, and red dots on the right represented upregulated genes. (B) Heatmap of DEGs: red represented high expression, and green represented low expression.

Identification of FRGs by WGCNA

A gene co-expression network based on FOXK1 was constructed using WGCNA on GSE221521. Initially, FOXK1 expression was analyzed in DR and control samples, revealing significant upregulation in DR samples (Figure 2A). Hierarchical clustering of all samples showed no outliers (Figure 2B). A β value of 6 was selected for constructing a scale-free network (Figure 2C). A hierarchical clustering tree was then constructed, resulting in eight co-expression modules (Figure 2D), with 4275 genes not assigned to any group and allocated to a grey module, which was excluded from further analysis. The MEblue module (cor = 0.62, p = 6e-14) was selected as the key module for subsequent analysis, containing 3408 key module genes (Figure 2E).

Figure 2 Construction of weighted gene co-expression network analysis (WGCNA). (A) Expression of FOXK1 in DR patients and control samples: yellow boxes represented DR patients, and red boxes represented control samples. **** indicates p < 0.0001. (B) Hierarchical clustering of samples: each vertical line at the bottom represented a sample, and the vertical lines at the top represented sample clusters. (C) Selection of the optimal soft-thresholding power: the left graph represented the scale-free fitting index, and the right graph represented the average connectivity. (D) Hierarchical clustering tree: the upper part was the hierarchical clustering dendrogram of genes, and the lower part was gene modules. (E) Heatmap of the correlation between modules and the FOXK1 gene: the values outside the parentheses represented correlation coefficients (larger values indicated stronger correlations), and the values inside the parentheses represented significance.

Identification and GO and KEGG Analysis of Candidate Genes

By analyzing the overlap between DEGs and key module genes, seven candidate genes were identified (Figure 3A). GO analysis revealed that the candidate genes were enriched in 87 terms, including 10 CCs, 17 MFs, and 60 BPs. Notable CCs included presynaptic active zone membrane and presynaptic active zone, MFs included ATP transmembrane transporter activity and purinergic nucleotide receptor activity, and BPs included protein maturation and bronchus development. The top five GO terms were bronchus development, lung secretory cell differentiation, negative regulation of cell fate commitment, positive regulation of cell fate commitment, and regulation of presynaptic cytosolic calcium ion concentration (Figure 3B). Additionally, KEGG pathway analysis showed that the candidate genes were enriched in ubiquinone and other terpenoid-quinone biosynthesis, phenylalanine metabolism, tyrosine metabolism, taste transduction, and biosynthesis of cofactors (Figure 3C).

Figure 3 Functional enrichment analysis. (A) Co-expressed genes of DEGs and module genes. (B) Gene Ontology (GO) enrichment analysis: red dots in the graph represented upregulation, blue dots represented downregulation, and the outermost circle showed the number of each term. (C) Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis: five different colors represented five different enrichment results.

SPDEF and SLC25A41 Identified as Biomarkers for DR

To identify the most diagnostically promising biomarker from seven candidate genes, LASSO regression analysis was applied. The optimal regularisation parameter lambda was determined via 5-fold cross-validation (lambda_(min) = 0.01772546). At this lambda value, the model retained six non-zero coefficient genes: AC099850.2, CISD3, HPD, SPDEF, SLC25A41, and P2RY4 (Figure 4A and B). Subsequently, the expression of these six potential biomarkers was validated in both the training dataset (GSE221521) and the independent validation dataset (GSE189005). Among these, SPDEF and SLC25A41 exhibited consistent expression patterns across both datasets, with significantly higher levels in the control group compared to the DR group (Figure 4C and D). The remaining four genes failed to demonstrate stable and significant differences across both cohorts. Consequently, SPDEF and SLC25A41 were identified as the final biomarkers.

Figure 4 Identification and validation of biomarkers. (A and B) LASSO regression analysis of candidate genes. (A) LASSO regularization path plot: it showed the model fitting effects corresponding to different values of the regularization parameter (λ) in the LASSO algorithm; the abscissa represented the value of the regularization parameter, and the ordinate represented the model performance index. (B) LASSO coefficient path plot: it showed the changes in feature coefficients under different values of the regularization parameter (λ); the abscissa represented the value of the regularization parameter, and the ordinate represented the absolute value or scaled value of feature coefficients. (C) Expression of biomarkers in the training cohort. (D) Expression of biomarkers in the validation cohort. * indicates p < 0.05; *** indicates p < 0.001; **** indicates p < 0.0001.

Constructing Biomarker-Related Nomogram

To evaluate the diagnostic value of biomarkers in clinical practice, a nomogram was constructed based on SPDEF and SLC25A41 (Figure 5A). Calibration curves were used to assess the nomogram’s effectiveness, and the predicted probability curves closely overlapped with the reference line, indicating that the nomogram has strong forecasting capabilities (Figure 5B). To further assess its practical value, the nomogram’s net benefit was analyzed. As shown in Figure 5C, when the risk threshold exceeded 0.2, the nomogram demonstrated a higher net benefit compared to both the all-line and none-line models, suggesting its practical utility.

Figure 5 Construction of the nomogram. (A) Nomogram constructed based on biomarkers: key genes were on the left; the lines corresponding to each key gene on the right were marked with scales (representing the available value range of the key gene), and the length of the line reflected the contribution of the gene to DR. (B) Calibration curve of the nomogram: the abscissa represented the predicted event rate, and the ordinate represented the observed actual event rate. (C) Decision curve analysis: the “None” curve represented no treatment for all individuals (net benefit = 0); the “All” curve represented treatment for all individuals (net benefit showed a reverse slash with a negative slope).

Immune Infiltration Analysis and GSEA

The composition of immune cells in DR and normal samples is shown in Figure 6A. Analysis of the expression of 22 immune cell types between DR and normal samples revealed higher levels of macrophages M0, monocytes, and activated CD4 memory T cells in the DR group (p-value < 0.05) (Figure 6B). Additionally, a correlation analysis between biomarkers and differential immune cells revealed a significant positive correlation between activated CD4 memory T cells and SPDEF, as well as between activated CD4 memory T cells and SLC25A41 (p-value < 0.05) (Figure 6C).

Figure 6 Immune cell infiltration analysis and gene set enrichment analysis (GSEA). (A) Distribution of immune cell abundance in the cohort: the stacked plot showed the proportion of 22 types of immune cells in each sample. (B) Differences in immune cells between the control group and the DR group: the abscissa represented cell types, and the ordinate represented immune infiltration scores; dark blue boxes represented DR patients, and yellow boxes represented control samples. (C) Correlation analysis between biomarkers and differential immune cell infiltration: red represented positive correlation, blue represented negative correlation, and the size of the fan shape represented the absolute value of the correlation coefficient. (D) GSEA results of the SPDEF gene. (E) GSEA results of the SLC25A41 gene.

GSEA analysis was performed to explore the pathways associated with the biomarkers. SPDEF was enriched in 21 functional pathways, with the top five enriched KEGG pathways being ribosome, Parkinson’s disease, oxidative phosphorylation, spliceosome, and olfactory transduction (Figure 6D). SLC25A41 was enriched in 53 functional pathways, with the top five being ribosome, Parkinson’s disease, oxidative phosphorylation, Huntington’s disease, and Alzheimer’s disease (Figure 6E).

Molecular Regulatory Networks

Through the miRTarBase and TarBase databases, hsa-let-7c-5p and hsa-mir-335-5p were identified as regulators of SPDEF, with hsa-mir-335-5p also regulating SLC25A41. Notably, hsa-mir-335-5p was a shared miRNA that regulates both biomarkers (Figure 7A and B). TFs obtained from the JASPAR database, including ARID3, FOXC1, and SREBF1 for SPDEF, and CREB1, FOXC1, PPARG, NFKB1, USF2, FOXF2, RUNX2, and NR2E3 for SLC25A41, revealed that FOXC1 co-regulated both biomarkers. A TF-miRNA-mRNA network consisting of 10 TFs, 2 miRNAs, and 2 mRNAs was constructed, comprising 14 nodes and 14 edges (Figure 7C).

Figure 7 Molecular regulatory network and drug prediction. (A) Key miRNAs of SLC25A41: blue represented miRNAs predicted by the TarBase database, and pink represented miRNAs predicted by the miRTarBase database. (B) Key miRNAs of SPDEF: blue represented miRNAs predicted by the TarBase database, and pink represented miRNAs predicted by the miRTarBase database. (C) TF-miRNA-mRNA network: red nodes in the graph represented biomarkers, yellow nodes represented miRNAs, and blue nodes represented transcription factors (TFs). (D) GeneMANIA network: the two biomarkers were in the middle, and the outer circle included genes functionally related to the biomarkers; lines of different colors represented different modes of action. (E) Biomarker-drug network: red nodes were biomarkers, and blue nodes were predicted drugs.

To identify further genes related to the biomarkers, 20 genes were identified through the GeneMANIA database, with these genes interacting through seven distinct modes of action (Figure 7D).

Drug Prediction

To explore potential therapeutic drugs for the biomarkers, a search in the DSigDB database yielded 24 drugs targeting SPDEF and one drug targeting SLC25A41. An mRNA-drug network comprising 27 nodes and 25 edges was subsequently constructed (Figure 7E). Resveratrol, daunorubicin, and digitoxigenin were identified as potential therapeutic agents for DR.

Quantitative PCR Analysis

For quantitative validation, RT-qPCR was performed on SPDEF and SLC25A41, using β-Actin as the reference gene. The stability of β-Actin was confirmed by the absence of a significant difference in its mean Ct values between the NC and DR groups (p > 0.05). The intra-group coefficients of variation (CV) for β-Actin Ct values were within an acceptable range for both genes: 5.97% (NC) and 8.46% (DR) for the SPDEF cohort, and 3.35% (NC) and 4.59% (DR) for the SLC25A41 cohort. After normalization, the results confirmed the significantly lower expression of both SPDEF and SLC25A41 in peripheral blood from patients with DR (p < 0.05) (Figure 8).

Figure 8 Validation of the expression level of SLC25A4 in clinical samples. * indicates p < 0.05.

Discussion

DR is a common complication in individuals with diabetes, often leading to severe visual impairment. FOXK1, a transcriptional regulator, mediates a wide array of biological processes. Through WGCNA, 3408 key module genes most strongly associated with DR and FOXK1 expression were identified. The intersection of DEGs and key module genes from WGCNA was subjected to GO and KEGG enrichment analysis. Subsequently, candidate genes were further refined using machine learning (LASSO), resulting in the identification of two biomarkers: SPDEF and SLC25A41.

SPDEF, a member of the Ets family of transcription factors, has been implicated in regulating cell differentiation, inflammatory responses, and tissue repair in various disease models.25,26 This study is the first to report its significant downregulation in the peripheral blood of patients with DR, suggesting that SPDEF may contribute to DR pathogenesis by modulating downstream target genes. Based on existing literature, we hypothesize that SPDEF potentially influences DR progression through two primary pathways. First, concerning vascular pathology, SPDEF has been reported to modulate the TGF-β signaling pathway.27 TGF-β is a key regulator of angiogenesis and vascular permeability,28 and retinal neovascularization coupled with breakdown of the blood-retinal barrier represent core pathological features of DR.29,30 Therefore, SPDEF might affect endothelial cell function by regulating TGF-β signaling, thereby participating in aberrant angiogenesis or vascular leakage in DR. Second, regarding metabolic stress, hyperglycemia-induced oxidative stress is a key initiating factor in DR, triggering apoptosis of retinal pericytes and endothelial cells through mechanisms such as the polyol pathway and accumulation of advanced glycation end products (AGEs).31 Notably, SPDEF has demonstrated regulatory capacity against oxidative stress in other systems, for instance, by modulating the expression of antioxidant enzymes like superoxide dismutase (SOD) to reduce ROS accumulation.32 This suggests that within the hyperglycemic microenvironment of DR, SPDEF might enhance the antioxidant defense system to alleviate oxidative damage, thereby protecting retinal cells. Nevertheless, the aforementioned mechanisms remain speculative. Whether SPDEF acts directly through the TGF-β or antioxidant pathways in DR requires further experimental validation, such as modulating SPDEF expression in DR cell or animal models and observing its effects on angiogenesis, oxidative stress markers, and apoptosis.

This study revealed a significant downregulation of SLC25A41 expression in the peripheral blood of patients with diabetic retinopathy (DR). As a member of the solute carrier family 25 (SLC25), SLC25A41 functions as a calcium-independent mitochondrial transporter primarily responsible for the translocation of adenylates (eg, ATP/ADP) across the mitochondrial membrane.33 This process is crucial for maintaining the mitochondrial membrane potential, energy homeostasis, and cell survival.34 Based on its molecular function and the known pathological mechanisms of DR, we hypothesize that SLC25A41 may contribute to DR pathogenesis by modulating mitochondrial function. Chronic hyperglycemia may lead to structural and functional abnormalities in mitochondria, including enhanced oxidative stress, calcium overload, and activation of the mitochondrial permeability transition (mPT), which in turn can trigger the release of mitochondrial contents (eg, cytochrome C), initiating apoptosis and inflammatory responses.35–37 It is noteworthy that SLC25A41-mediated adenylate transport is closely associated with mitochondrial membrane stability and the regulation of mPT.34 Therefore, the downregulation of SLC25A41 might disrupt mitochondrial energy metabolism and membrane permeability, exacerbating mitochondrial damage induced by high glucose. Furthermore, mitophagy, the selective autophagy of damaged mitochondria, which is essential for maintaining oxidative phosphorylation and ATP synthesis, may represent an important cytoprotective mechanism in DR.38 SLC25A41 could potentially influence the mitochondrial membrane potential and the efficiency of ADP/ATP exchange, thereby indirectly regulating the mitophagy process. Its decreased expression might impair mitochondrial quality control, leading to the accumulation of dysfunctional mitochondria and potentially promoting the progression of DR. In summary, we speculate that SLC25A41 may influence DR through the following mechanisms: (1) regulating mitochondrial energy metabolism and membrane stability, thereby affecting mPT-related apoptotic pathways; and (2) modulating mitophagy to maintain mitochondrial quality and function. These hypotheses provide a mechanistic framework for understanding the role of SLC25A41 in DR, yet they require further validation through genetic manipulations, cellular models, and animal studies to elucidate the specific regulatory pathways involved.

Immune infiltration analysis revealed that macrophages M0, monocytes, and activated CD4 memory T cells were highly expressed in the DR group. Studies have shown that the early destruction of the BRB in DR facilitates the entry of blood monocytes into the retina. The inflammatory immune response mediated by macrophages plays a pivotal role in the pathogenesis and progression of DR. Several studies have highlighted macrophages as key drivers in the development of proliferative DR (PDR).39–41 Macrophages induce chemotaxis and fibroplasia by secreting leukotrienes and fibronectin. They also influence cell proliferation through the synthesis of VEGF, PDGF, fibroblast growth factor (FGF), and TGF-β.39 Activated macrophages secrete excessive inflammatory cytokines, such as IL-1β, TNF-α, IL-6, and IL-12,42 through the NF-κB signaling pathway. Additionally, the phagocytic function of macrophages is impaired, suggesting that macrophage dysfunction may exacerbate the inflammatory response in DR.

The increase in M0 macrophages and activated CD4 memory T cells observed in this study aligns with recent clinical evidence underscoring the central role of systemic inflammation in DR. For instance, the study by Aykut Bulu et al indicated a significant correlation between circulating levels of specific pro-inflammatory monocyte subsets and DR severity, further supporting the contribution of innate immune activation to DR pathogenesis.43 Our transcriptomic findings and such clinical evidence are mutually reinforcing, which not only strengthens the biological plausibility of immune dysregulation as a key driver of DR but also suggests that these immune features may serve as potential clinical indicators for assessing disease activity.

CD4⁺T cells are a central cell population in regulating adaptive immune responses. They can differentiate into subsets such as Th1 and Th17, secrete cytokines including IFN-γ and IL-17, and participate in inflammatory processes.44 In DR, CD4⁺T cells may contribute to disease progression through multiple mechanisms. First, they may directly engage in retinal vascular inflammation by releasing pro-inflammatory cytokines that disrupt the integrity of the blood-retinal barrier, thereby exacerbating vascular leakage.45 Furthermore, angiogenesis-related receptors such as VEGFR have been shown to modulate the expression of CD4⁺T cells,46 suggesting that these lymphocytes may be influenced by VEGFR signaling and participate in the recruitment and stabilization of new blood vessels.47 Based on the significant positive correlations observed in this study between activated CD4⁺ memory T cells and both SPDEF and SLC25A41, it is hypothesized that these two biomarkers may cooperate with CD4⁺T cells in the regulation of DR. This finding offers a new perspective on the role of CD4⁺T cells in DR and implies that SPDEF and SLC25A41 may play important roles in modulating both immune responses and angiogenesis. From a potential clinical application standpoint, the associations among immune infiltration features and these biomarkers may provide insights for immune-targeted interventions in DR. For instance, strategies aimed at modulating the polarization of M0 macrophages or restoring the functional balance of activated CD4⁺ memory T cells could help alleviate local retinal inflammation and pathological angiogenesis. Moreover, as molecules associated with the functions of these immune cells, the expression levels of SPDEF and SLC25A41 may serve as auxiliary indicators for assessing the immune-inflammatory status in DR, thereby supporting clinical evaluation of disease progression or treatment efficacy. However, these potential applications require further clinical investigation to validate their feasibility and effectiveness in actual diagnosis and therapy.

In this study, a total of 74 functional pathways were identified from the two biomarkers. Common KEGG pathways between SPDEF and SLC25A41 included the ribosome, Parkinson’s disease, and oxidative phosphorylation. Oxidative phosphorylation plays a pivotal role in energy production within organisms, and genes associated with Parkinson’s disease are also linked to oxidative phosphorylation.48 This pathway is involved in the maintenance of glucose homeostasis and intricately interacts with the glycolytic Krebs cycle. Disruption of the balance between oxidative phosphorylation and glycolysis can lead to biochemical and molecular changes associated with DR. These changes include miscommunication between the endoplasmic reticulum (ER) and mitochondria, as well as dysregulation of mitochondrial autophagy, both of which contribute to DR development.49 Genipin has been shown to protect cells and tissues from high glucose levels by regulating HIF-1α phosphorylation, which enhances oxidative phosphorylation while reducing substrate phosphorylation. Activation of the HIF-1α pathway may also offer a potential therapeutic approach for DR.50 This insight into key genes and pathways contributes to a better understanding of DR prognosis and pathogenesis.

In drug prediction, daunorubicin was identified as a potential therapeutic agent. Currently used for treating adult acute myeloid leukemia,51,52 daunorubicin has been shown to inhibit the transcriptional activity of hypoxia-inducible factor-1 (HIF-1) by blocking its binding to DNA. Intravitreal injection of daunorubicin (DXR or DNR) has been found to inhibit choroidal and retinal neovascularization (NV) and can be fabricated as nanomaterials to reduce intraocular toxic effects.53 Resveratrol, a natural antioxidant, reduces blood viscosity, inhibits platelet aggregation, induces vasodilation, and promotes blood flow. As a phenolic phytochemical, resveratrol has broad therapeutic potential, particularly in diseases related to oxidative stress. It has demonstrated efficacy in attenuating retinopathy and neuropathy.54

In the molecular regulatory network, hsa-mir-335-5p serves as a co-regulator of SPDEF and SLC25A41, and it has been investigated as a potential novel biomarker in rheumatoid arthritis.55 miRNAs are short, single-stranded non-coding RNAs that regulate gene expression by hybridizing with messenger RNA or promoter sequences.56,57 Dysregulation of miRNA expression is linked to various diseases, and specific miRNA expression profiles have been proposed as diagnostic biomarkers.58–60 Furthermore, miRNAs are increasingly recognized as therapeutic targets for several diseases, including retinal vascular diseases.61,62 Therefore, miRNAs hold significant potential as biomarkers and predictive indicators for ocular diseases.

This study, utilizing GSE221521 and FOXK1 genes, identified biomarkers associated with FOXK1 and established a nomogram with promising clinical predictive capability. The nomogram developed in this study demonstrated favorable performance in internal calibration and decision curve analysis, indicating its potential to assist in clinical decision-making. From the perspective of clinical applicability, the model is intended to serve as an auxiliary risk stratification tool, which can be used to identify patients at high risk of developing DR, thereby aiding clinicians in early screening and focused follow-up. It must be emphasized that the current performance evaluation of the model is based on the existing training and validation cohorts, and its broader clinical application still requires prospective validation in larger-scale, multi-center, external independent cohorts.

Additionally, immune infiltration analysis, GSEA enrichment analysis, molecular regulatory network analysis, and drug prediction were conducted to explore novel therapeutic mechanisms of DR based on the FOXK1 gene at the molecular level. However, this study has some limitations. First, although we identified FOXK1-associated genes through WGCNA and recognized SLC25A41 and SPDEF as potential biomarkers using bioinformatic approaches, the lack of experimental validation prevents us from establishing a causal relationship between FOXK1 and SLC25A41/SPDEF. Second, although RT-qPCR confirmed the expression levels of these biomarkers, the small sample size may limit the statistical power of the study. Furthermore, the absence of functional assays and in vivo experiments hinders the exploration of their actual roles in the pathogenesis of DR. In addition, the diagnostic and clinical value of these biomarkers remains unclear due to the lack of independent cohort validation, clinical performance evaluation, and functional characterization. The predicted potential drugs derived from drug prediction analysis also require validation to confirm the direct link between the biomarkers and drug response. In future studies, we aim to investigate the direct regulatory relationship between FOXK1 and SLC25A41/SPDEF using methods such as dual-luciferase reporter assays and co-immunoprecipitation (Co-IP). We also plan to conduct independent multi-center clinical collaborations to increase the sample size and enhance statistical power, while systematically evaluating the clinical diagnostic performance of the biomarkers and the diagnostic model. Moreover, we intend to establish DR cell and animal models for gene knockout/overexpression experiments to elucidate the functional mechanisms of these key biomarkers. In vitro and in vivo drug experiments will be carried out to validate the therapeutic effects of the predicted compounds. Combined with molecular biology techniques, these studies will further decipher the mechanism by which FOXK1 regulates DR, thereby providing a more solid experimental basis for the research, diagnosis, and treatment of DR.

Conclusion

In conclusion, this integrated bioinformatics and experimental study identifies SPDEF and SLC25A41 as promising diagnostic biomarkers for diabetic retinopathy (DR), with their expression validated. The construction of a predictive nomogram and the exploration of associated immune cell infiltration, regulatory networks (including the FOXC1/hsa-mir-335-5p axis), and key pathways like oxidative phosphorylation provide a multi-faceted perspective on their potential role in DR pathogenesis. Furthermore, the prediction of targeting drugs, such as resveratrol, offers translational potential. These biomarkers and their accompanying prediction models are expected to enhance the early detection and risk stratification capabilities for Diabetic Retinopathy (DR) in clinical practice. To translate these research findings into tangible clinical applications, future work must prioritize rigorous validation through large-scale, prospective, multi-ethnic cohort studies, coupled with functional investigations to clearly establish the roles of these biomarkers in the pathophysiological mechanisms of DR.

Abbreviations

DR, Diabetic retinopathy; FOXK1, Forkhead box K1; WGCNA, Weighted Gene co-expression network analysis; LASSO, Least absolute shrinkage and selection operator; GSEA, Gene set enrichment analysis; RT-qPCR, Reverse transcription quantitative polymerase chain reaction; TFs, Transcription factor microRNA; miRNAs, microRNA; PTM, Post-translational modification; HUVECs, Human umbilical vein endothelial cells; GEO, Gene expression omnibus; DEGs, Differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto encyclopedia of genes and genomes; CC, Cellular components; BC, Biological processes; MF, Molecular functions; FGF, Fibroblast growth factor; HIF-1, Hypoxia-inducible factor-1.

Data Sharing Statement

The data that support the findings of this study are available from the first author, Xin Li, upon reasonable request.

Ethics Approval and Consent to Participate

This study was performed in accordance with the Declaration of Helsinki and was approved by the Ethics Committee of Wuhan Aier Eye Hanyang Hospital (IEC-AD-HYEYE2024008). All participants provided written informed consent.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Funding

The Funding for Scientific Research Projects from Wuhan Municipal Health Commission (No. WX23A27).

Disclosure

No conflicts of interest, financial or otherwise, are declared by the authors.

References

1. Zhu J, Huang J, Sun Y, Xu W, Qian H. Emerging role of extracellular vesicles in diabetic retinopathy. Theranostics. 2024;14(4):1631–17. doi:10.7150/thno.92463

2. Teo ZL, Tham YC, Yu M, et al. Global prevalence of diabetic retinopathy and projection of burden through 2045: systematic review and meta-analysis. Ophthalmology. 2021;128(11):1580–1591. doi:10.1016/j.ophtha.2021.04.027

3. Li SY, Zhao N, Wei D, et al. Ferroptosis in the ageing retina: a malevolent fire of diabetic retinopathy. Ageing Res Rev. 2024;93:102142. doi:10.1016/j.arr.2023.102142

4. Kollias AN, Ulbig MW. Diabetic retinopathy: early diagnosis and effective treatment. Dtsch Arztebl Int. 2010;107(5):75–83. quiz 4. doi:10.3238/arztebl.2010.0075

5. Amoaku WM, Ghanchi F, Bailey C, et al. Diabetic retinopathy and diabetic macular oedema pathways and management: UK consensus working group. Eye. 2020;34(Suppl 1):1–51. doi:10.1038/s41433-020-0961-6

6. Tomkins-Netzer O, Niederer R, Lightman S. The role of statins in diabetic retinopathy. Trends Cardiovasc Med. 2024;34(2):128–135. doi:10.1016/j.tcm.2022.11.003

7. Tan TE, Wong TY. Diabetic retinopathy: looking forward to 2030. Front Endocrinol. 2022;13:1077669. doi:10.3389/fendo.2022.1077669

8. Liu Y, Ding W, Ge H, et al. FOXK transcription factors: regulation and critical role in cancer. Cancer Lett. 2019;458:1–12. doi:10.1016/j.canlet.2019.05.030

9. Bowman CJ, Ayer DE, Dynlacht BD. Foxk proteins repress the initiation of starvation-induced atrophy and autophagy programs. Nat Cell Biol. 2014;16(12):1202–1214. doi:10.1038/ncb3062

10. Gong Q, Luo D, Wang H, et al. Inhibiting autophagy by miR-19a-3p/PTEN regulation protected retinal pigment epithelial cells from hyperglycemic damage. Biochim Biophys Acta Mol Cell Res. 2023;1870(7):119530. doi:10.1016/j.bbamcr.2023.119530

11. Liu M, Gong S, Sheng X, Zhang Z, Yang J. Mechanism of TNF- α inducing apoptosis and autophagy of chondrocytes by activating NF- κ B signal pathway. Cell Mol Biol. 2023;69(15):95–98. doi:10.14715/cmb/2023.69.15.16

12. Tan B, Chen T, Song P, et al. Interaction of inflammatory factors related to pulmonary infection and TLR4/NF-κB signaling pathway in patients with spontaneous intracerebral hemorrhage. Cytokine. 2025;186:156851. doi:10.1016/j.cyto.2024.156851

13. Yang S, Li F, Lu S, et al. Ginseng root extract attenuates inflammation by inhibiting the MAPK/NF-κB signaling pathway and activating autophagy and p62-Nrf2-Keap1 signaling in vitro and in vivo. J Ethnopharmacol. 2022;283:114739. doi:10.1016/j.jep.2021.114739

14. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550. doi:10.1186/s13059-014-0550-8

15. Gustavsson EK, Zhang D, Reynolds RH, Garcia-Ruiz S, Ryten M. ggtranscript: an R package for the visualization and interpretation of transcript isoforms using ggplot2. Bioinformatics. 2022;38(15):3844–3846. doi:10.1093/bioinformatics/btac409

16. Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32(18):2847–2849. doi:10.1093/bioinformatics/btw313

17. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinf. 2008;9(1):559. doi:10.1186/1471-2105-9-559

18. Chen H, Boutros PC. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinf. 2011;12(1):35. doi:10.1186/1471-2105-12-35

19. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation. 2021;2(3):100141. doi:10.1016/j.xinn.2021.100141

20. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010;33(1):1–22. doi:10.18637/jss.v033.i01

21. Silva-Martins M, Canário AC, Abreu-Lima I, Krasniqi L, Cruz O. Psychometric properties of the Portuguese version of the physical activity parenting practices questionnaire. BMC Psychol. 2023;11(1):417. doi:10.1186/s40359-023-01444-4

22. Wang L, Wang D, Yang L, et al. Cuproptosis related genes associated with Jab1 shapes tumor microenvironment and pharmacological profile in nasopharyngeal carcinoma. Front Immunol. 2022;13:989286. doi:10.3389/fimmu.2022.989286

23. Qin Y, Zu X, Li Y, et al. A cancer-associated fibroblast subtypes-based signature enables the evaluation of immunotherapy response and prognosis in bladder cancer. Iscience. 2023;26(9):107722. doi:10.1016/j.isci.2023.107722

24. Szklarczyk D, Morris JH, Cook H, et al. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 2017;45(D1):D362–d8. doi:10.1093/nar/gkw937

25. Rajavelu P, Chen G, Xu Y, et al. Airway epithelial SPDEF integrates goblet cell differentiation and pulmonary Th2 inflammation. J Clin Invest. 2015;125(5):2021–2031. doi:10.1172/jci79422

26. Fang R, Uchiyama R, Sakai S, et al. ASC and NLRP3 maintain innate immune homeostasis in the airway through an inflammasome-independent mechanism. Mucosal Immunol. 2019;12(5):1092–1103. doi:10.1038/s41385-019-0181-1

27. Kramer EL, Hardie WD, Madala SK, Davidson C, Clancy JP. Subacute TGFβ expression drives inflammation, goblet cell hyperplasia, and pulmonary function abnormalities in mice with effects dependent on CFTR function. Am J Physiol Lung Cell Mol Physiol. 2018;315(3):L456–l65. doi:10.1152/ajplung.00530.2017

28. Zhao M, Hu Y, Jin J, et al. Interleukin 37 promotes angiogenesis through TGF-β signaling. Sci Rep. 2017;7(1):6113. doi:10.1038/s41598-017-06124-z

29. Arrigo A, Aragona E, Bandello F. VEGF-targeting drugs for the treatment of retinal neovascularization in diabetic retinopathy. Ann Med. 2022;54(1):1089–1111. doi:10.1080/07853890.2022.2064541

30. Wang X, He X, Li Z, et al. Insight into dysregulated VEGF-related genes in diabetic retinopathy through bioinformatic analyses. Naunyn Schmiedebergs Arch Pharmacol. 2025;398(6):7199–7217. doi:10.1007/s00210-024-03638-y

31. Hussain A, Ashique S, Afzal O, et al. A correlation between oxidative stress and diabetic retinopathy: an updated review. Exp Eye Res. 2023;236:109650. doi:10.1016/j.exer.2023.109650

32. Kojima T, Dogru M, Ibrahim OM, et al. Effects of oxidative stress on the conjunctiva in Cu, Zn-Superoxide dismutase-1 (Sod1)-knockout mice. Invest Ophthalmol Vis Sci. 2015;56(13):8382–8391. doi:10.1167/iovs.15-18295

33. Traba J, Satrústegui J, Del Arco A. Characterization of SCaMC-3-like/slc25a41, a novel calcium-independent mitochondrial ATP-Mg/Pi carrier. Biochem J. 2009;418(1):125–133. doi:10.1042/bj20081262

34. Brenner C, Subramaniam K, Pertuiset C, Pervaiz S. Adenine nucleotide translocase family: four isoforms for apoptosis modulation in cancer. Oncogene. 2011;30(8):883–895. doi:10.1038/onc.2010.501

35. Adebayo M, Singh S, Singh AP, Dasgupta S. Mitochondrial fusion and fission: the fine-tune balance for cellular homeostasis. FASEB j. 2021;35(6):e21620. doi:10.1096/fj.202100067R

36. Vercellino I, Sazanov LA. The assembly, regulation and function of the mitochondrial respiratory chain. Nat Rev Mol Cell Biol. 2022;23(2):141–161. doi:10.1038/s41580-021-00415-0

37. Bonora M, Giorgi C, Pinton P. Molecular mechanisms and consequences of mitochondrial permeability transition. Nat Rev Mol Cell Biol. 2022;23(4):266–285. doi:10.1038/s41580-021-00433-y

38. Qiao S, Dennis M, Song X, et al. A REDD1/TXNIP pro-oxidant complex regulates ATG4B activity to control stress-induced autophagy and sustain exercise capacity. Nat Commun. 2015;6(1):7014. doi:10.1038/ncomms8014

39. Pavlou S, Lindsay J, Ingram R, Xu H, Chen M. Sustained high glucose exposure sensitizes macrophage responses to cytokine stimuli but reduces their phagocytic activity. BMC Immunol. 2018;19(1):24. doi:10.1186/s12865-018-0261-0

40. Al-Rashed F, Sindhu S, Arefanian H, et al. Repetitive intermittent hyperglycemia drives the M1 polarization and inflammatory responses in THP-1 macrophages through the mechanism involving the TLR4-IRF5 pathway. Cells. 2020;9(8):1892. doi:10.3390/cells9081892

41. Jia Y, Zhou Y. Involvement of lncRNAs and macrophages: potential regulatory link to angiogenesis. J Immunol Res. 2020;2020(1):1704631. doi:10.1155/2020/1704631

42. Cheng CI, Chen PH, Lin YC, Kao YH. High glucose activates Raw264.7 macrophages through RhoA kinase-mediated signaling pathway. Cell Signal. 2015;27(2):283–292. doi:10.1016/j.cellsig.2014.11.012