Two independent PDR datasets were downloaded from the GEO database. (http://www.ncbi.nlm.nih.gov/geo/), which included GSE102485 and GSE60436 ( Table 1 ). The first transcriptome dataset was the test dataset (GSE102485) and the second microarray dataset was the validation dataset (GSE60436). For GSE102485, we selected a subset of 22 neovascular proliferative membrane specimens and three normal retina samples (14). We downloaded and analyzed the raw data from the GSE102485 dataset, processed and normalized the protein-coding genes using the R package “DESeq2” for further analysis. Additionally, we selected the GSE60436 microarray dataset as the validation set, comprising 6 PDR samples and 3 normal samples (15). The raw data from the GSE60436 dataset were downloaded, underwent ID conversion, normalization, and background correction. Additionally, we identified 328 ER stress-related genes with a relevance score greater than 3 in Genecards ( Supplementary Table S1 ).

Before conducting differential analysis on GSE102485, we filtered for mRNAs with expression counts greater than one in at least the number of replicates. Subsequently, we analyzed gene expression using the DESeq2 package. DEGs were identified using the criteria of an adjusted P-value <0.05 and an absolute |log2 Fold Change| ≥ 2. Next, we employ the R software packages “heatmap” and “ggplot2” to create visual representations, including heatmap and volcano plots. We performed an intersection analysis and created a Venn diagram to visually represent the overlap among these ER stress-related genes and DEGs in GSE102485.

We conducted GO enrichment analysis and KEGG pathway analysis on the ER stress-related DEGs using the R package “ClusterProfiler” (16). We considered adjusted p-values of < 0.05 to be statistically significant in our analysis. Furthermore, we employed the Metascape database (https://metascape.org/) to further explore the functional mechanisms. The criteria set were a minimum overlap of 3, p ≤ 0.01, and a minimum enrichment of 1.5.

We used the STRING database (https://string-db.org/) to explore the interactions among the ER stress-related DEGs (17). A threshold of a combined score ≥ 0.4 was set to identify significant interactions between these genes. The Cytoscape plugin MCODE was employed to filter out important modules of core genes in the PPI network with the following criteria: degree cutoff = 2, node score cutoff = 0.2, K-core = 2, maximum depth = 100. To identify hub genes, we utilized the CytoHubba plug-in (version 0.1) within Cytoscape software. CytoHubba identifies hub genes based on the consensus of multiple algorithms (MCC, MNC, Degree and EPC).

We utilized the CIBERSOFT algorithm from the R package “IOBR” (18) to analyze the infiltration of multiple immune cell subtypes in PDR and control samples from the GSE102485 dataset with the transcriptome of neovascular membranes. Linear regression analysis was conducted to analyze the correlation between the expression of ER stress-related hub genes and immune cells. The results were visualized using the R package “ggplot2”.

Human umbilical vein endothelial cells (HUVECs) (ATCC, Cat. CRL-1730) were cultured in DMEM containing 10% fetal bovine serum (FBS) and 1% antibiotic-antimycotic under standard conditions (5% CO2, 37°C). The HUVECs in the logarithmic growth phase were categorized into two groups for the experiment: the control group and the high glucose (HG) group. HUVECs cultured in a medium with 5.5 mmol/L glucose were assigned to the control group, whereas those cultured in a medium containing 30 mmol/L glucose were assigned to the HG group.

Protein extraction from HUVECs was performed using RIPA lysis buffer containing protease inhibitors and phosphatase inhibitors (Yesen, China). The protein concentration was determined using the BCA protein assay kit (ZJ102, Epizyme, China). Subsequently, proteins were separated on 10% SDS-PAGE gels and transferred onto PVDF membranes (R9A84148, Millipore). The membranes were then blocked in QuickBlock™ Western Blocking Buffer (P0252, Beyotime) for 30 minutes and incubated overnight at 4°C with primary antibodies against TRAM1 (18243-1-AP, Proteintech) and TXNIP (12705-1-AP, Proteintech). The following day, the membranes were incubated with HRP-conjugated secondary antibodies at room temperature for 1 hour. Protein bands were visualized using enhanced chemiluminescence (SQ201, EpiZyme, China), and protein densitometry was quantified using ImageJ software (version 6.0; Media Cybernetics, Inc.). β-actin (66031-1-lg, Proteintech) was used as the internal reference.

RNA was extracted from HUVECs using the EZ-press-RNA purification kit (EZBioscience, USA) according to the manufacturer’s plan. Then, cDNA was reverse-transcribed from total RNA using the reverse transcription kit (Takara, Japan) and RT-PCR was conducted using TB Green^® Premix Ex Taq™ II kit (Takara, Japan). Primer sequences are detailed in Supplementary Table S2 . β-actin was chosen as the reference gene for normalizing mRNA expression levels, and quantitative analysis was performed using the 2^-ΔΔCT method.

The Connectivity Map (cMAP) website was utilized to investigate small molecule drugs with the potential to inhibit the formation and progression of PDR (19). We submitted the 51 ER stress-related DEGs to the cMAP website, focusing specifically on the 46 upregulated genes, to identify potential small molecule drugs that could inhibit the formation and development of PDR. The score in the results list returned by cMAP represents the percentage by which the reference gene set is more similar to the current perturbation compared to the similarity of the query to the current perturbation. Drugs with negative scores and high absolute values are considered potential treatments because they can inhibit the expression of ER stress characteristic genes.

An independent Student’s t-test was employed to compare the two groups based on statistically significant differences of normally distributed variables. For non-normally distributed variables, the Wilcoxon rank-sum test was utilized. Statistical analysis was performed using R (version 4.2.0), with a significance threshold set at p < 0.05.

The study was designed according to the flow chart outlined in Figure 1 . Overall, 2451 DEGs were recognized in GSE102485, with 1815 genes displaying significant upregulation and 636 genes exhibiting significant downregulation. The volcano plots illustrating the DEGs are presented in Figure 2A . Subsequently, employing a filtering criterion of a correlation score > 3, we identified a total of 328 ER stress-related genes from the GeneCards database. By generating Venn diagrams, a total of 51 ER stress-related DEGs were identified, among which 46 were upregulated and 5 were downregulated as ER stress-related DEGs ( Figures 2B, C ).

To further explore the potential biological functions of these ER stress-related DEGs at the biological level, we conducted GO and KEGG analyses. As depicted in Figures 3A–C , Biological Process (BP) terms were predominantly enriched in “response to endoplasmic reticulum stress”, “intrinsic apoptotic signaling pathway in response to endoplasmic reticulum stress”, and “intrinsic apoptotic signaling pathway”. In the context of Cellular Component (CC) ontology, significant enrichment was noted in both the “outer membrane” and the “mitochondrial outer membrane”. Turning to Molecular Function (MF) analysis, the predominant enrichment was identified in functions related to “cytokine receptor binding” and “protein phosphatase 2A binding”. Moreover, illustrated in Figure 3D , the KEGG analysis showcased enrichment in pathways such as the “AGE-RAGE signaling pathway in diabetic complications”, “NOD-like receptor signaling pathway” and the “TNF signaling pathway”. As illustrated in Figure 3E , the bar chart depicts the results of Metascape enrichment analysis for the provided gene list. The enrichment analysis of ER stress-related DEGs revealed significant associations with biological processes such as lipid and atherosclerosis pathways, immune responses including NOD-like receptor signaling and infectious diseases like influenza. Additionally, ER stress-related DEGs were implicated in cellular responses to various stimuli, including abiotic and mechanical stress, as well as pathways related to ER stress and apoptotic signaling.

To further elucidate the potential relationships among the proteins encoded by these ER stress-related DEGs and to identify hub genes, a PPI network analysis was conducted using STRING. The PPI network comprised 51 nodes and 283 edges, with a highly significant enrichment (PPI enrichment p-value < 1.0e-16) as shown in Figure 4A . Subsequently, module analysis using the MCODE plugin revealed the most significant module: Module 1, which included 18 nodes and 140 edges, with a cluster score of 16.471 ( Figure 4B ).

In the quest to pinpoint the hub genes among the ER stress-related DEGs, multiple topological analysis algorithms, including MCC, MNC, Degree, and EPC, were employed. The results from the top 10 genes obtained from each algorithm were cross-referenced, leading to the identification of 6 hub genes: BCL2, CCL2, IL-1β, TLR4, TNF, and TP53 ( Figure 4C ). The location of the 6 hub genes on chromosomes is shown in Figure 4D .

We used the CIBERSORT algorithm to assess the proportions of different infiltrating immune cell types between the PDR group and the control group. The bar chart presented in Figure 5A illustrates the proportions of 22 immune cell types across the 25 samples. Compared to the control group, there was an increase in eosinophil infiltration while memory B cells and T follicular helper cell infiltration decreased ( Figure 5B ). Concerning the correlation between hub gene expression and immune cell infiltration, there was a notable negative correlation between the expression of BCL2, CCL2, IL-1β, TLR4, and TP53, and the levels of infiltration of memory B cells and T follicular helper cells ( Figure 5C ).

Typically, there is a degree of consistency observed between the outcomes of the training set and those of the validation set. After normalizing the raw data from the validation set GSE60436, we found that in the comparison between the PDR group and the control group, the expression differences of CCL2, IL-1β, TLR4, TNF, and TP53 were consistent with our results from the GSE102485 dataset, and these differences were statistically significant with a p-value < 0.05 ( Figure 6A ).

In addition, we cultured HUVECs in a HG environment (30mM) to simulate the DR model in vitro. Firstly, we examined the expression of the ER stress markers TRAM1 and TXNIP proteins using Western blotting, which are key executive factors of ER stress. The results revealed a significant upregulation of TRAM1 and TXNIP protein levels in HUVECs after incubation with high glucose ( Figure 6B ), indicating the occurrence of ER stress. Subsequently, we employed qRT-PCR analysis to validate the expression of ER stress-related hub genes. The results revealed a significant upregulation in the expression of BCL2, CCL2, IL-1β, TLR4, TNF, and TP53 in HUVECs after 48 hours under HG conditions compared to the low-glucose environment ( Supplementary Table S3 ). These findings align with the results obtained from bioinformatics analysis ( Figure 6C ).

To predict potential small molecule drugs that may inhibit ER stress in PDR, we uploaded the upregulated ER stress-DEGs to the cMAP online tool. We identified 8 drugs with the highest negative scores (diazepam, FG-7142, benzanthrone, AR-A014418, rucaparib, phenamil, quercetagetin, and parbendazole), indicating that they may inhibit the expression of ER stress markers ( Figure 7A , Supplementary Table S4 ). Furthermore, we present the chemical structures of these eight small molecular compounds ( Figure 7B ).