Multiple Hub Genes as Diagnostic and Therapeutic Targets Inducing Inflammation and Angiogenesis in Retinopathy of Prematurity

Introduction

Retinopathy of prematurity (ROP) is a significant contributor to blindness that can be treated in preterm infants worldwide.¹ It results from immature retinal blood vessels and is one of the main complications facing the survival of premature infants. The pathological process of ROP mainly affects the fundus of the eye, and the degree of the disease can be divided into five different stages, from mild vascular dilation and thinning to severe neovascularization and retinal detachment. Accordingly, it is greatly important to clarify the etiopathogenesis and identify key biomarkers of ROP for personalized and effective treatment.

The etiopathogenesis of ROP is not yet clear. One of the characteristics of ROP is pathological blood vessels that develop in a disorderly fashion into the vitreous, where they can cause retinal traction, detachment, and bleeding, ultimately leading to blindness, which is prompted by growing local vascular endothelial growth factor (VEGF) levels from peripheral avascular retina.² Recent evidences suggested that inflammation could make contributions to a continuous rise in the risk of ROP.³ Recent researches have indicated macrophages, monocytes and microglia and so on play important roles in the development of intraocular neovascularization.⁴ Studies have explored how chemokines are involved in angiogenesis, growth control and hematopoiesis.⁵ For example, interleukin-8 (IL-8), a chemokine, is involved in inflammation and pathological blood vessel formation in the eye.⁶ Nonetheless, our comprehension of ROP remains quite restricted. Currently, there are main methods for treating ROP, including laser photocoagulation, cryotherapy, scleral buckling or vitrectomy, and anti-VEGF drugs. Laser photocoagulation aims at destroying the pathological neovascularization. However, laser photocoagulation treatment comes with side effects like subjective pain, ciliary and exudative retinal detachment, and macular edema.⁷ Cryotherapy is the conventional treatment for ROP in the past, which employs a probe chilled to subzero temperatures, and which necessitates general anesthesia.⁷ Cryotherapy has not been widely used to treat ROP since the late 1980s because it induces more inflammation than laser therapy.⁷ Scleral buckling or vitrectomy is used to alleviate vitreoretinal traction and can reattach the retina to preserve vision and prevent blindness but can be associated with significant ocular side effects such as myopia.⁷ Four drugs, administered intravitreally as anti-VEGF agents, have been applied for ROP treatment: ranibizumab, bevacizumab, aflibercept and conbercept. However, anti-VEGF drugs make the possibility of harmful side effects due to undifferentiated suppression of every form of angiogenesis, like increased arterial pressure, endothelial cell detachment and thrombotic events.⁸ It can thus be seen that exploring more complete etiopathogenesis and more important therapeutic targets is of great significance.

In our research, three hub genes Plxnd1, Esm1 and Cd248 were screened out by WGCNA and various machine learning methods based on the dataset in the Gene Expression Omnibus (GEO),⁹ and validated through nomograms, ROC curves and qRT-PCR. The relationships between these hub genes, immune reaction and ROP were explored by means of immune cell infiltration analysis. The ceRNA network revealed regulatory interactions between these hub genes and their related miRNAs and lncRNAs. Via drug–gene interaction networks, the potential therapeutic drugs were screened for these hub genes. Overall, the 3 hub genes were correlated with immune response and angiogenesis during the occurrence and development of ROP. Our research is the first to explore the hub genes of ROP through bioinformatics and artificial intelligence methodologies (including Least Absolute Shrinkage and Selection Operator (LASSO), Support Vector Machine (SVM), and Random Forest (RF)) and elucidate diagnostic and therapeutic value of the 3 hub genes, which will offer novel insights into for the nosogenesis and treatment of ROP.

Materials and Methods

Data Collection and Preprocessing

Analysis process of this research is shown in Figure 1. Because clinical samples of infant retina are unavailable, mouse model of oxygen-induced retinopathy (OIR) has been widely employed for ROP-related researches. The classical mouse model of OIR has proven useful in imitating the processes in prematurity retinopathy and abnormal angiogenesis,^10–14 and can be used to explore the pathogenesis of ROP. In this classical OIR model, mice were exposed to 75% oxygen for 5 consecutive days at postnatal day (P) 7. At P12, mice were taken out of the oxygen chamber and brought back to normal air conditions. At P 17, abundant new vascularization occurred, and proliferative retinopathy developed reliably (and quantifiably) over 17 days.¹⁵ Accordingly, we selected and obtained, namely, GSE130400 from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), which was based on classical mouse OIR models.¹⁶ Furthermore, this dataset had expression data of whole genome and the more experimental samples compared with other datasets, including 13 OIR mouse samples (biological replicates) and 12 control mouse samples (biological replicates). GSE130400 was processed and analyzed using the dplyr package.

Figure 1 The workflow of our research.

Differentially Expressed Analysis

Differentially expressed analysis (|log₂FC| ≥ 1, adjusted p-value < 0.05) was conducted via comparing oxygen-induced retinopathy (OIR) tissues to normal tissues in the R computing environment employing the Linear Model Fitting (LimFit) package and the Empirical Bayes (eBayes) of the limma package. The data analysis was conducted using GraphPad Prism, and a volcano plot was created to display differentially expressed genes (DEGs).

WGCNA Analysis

WGCNA is conducted to divide gene modules through calculating correlation between gene expression levels and clinical features. We can identify candidate biomarkers or therapeutic targets through WGCNA. In this study, WGCNA was used to obtain co-expression modules and genes connected with OIR. A scale-free analysis was carried out to identify whether interactions between genes followed scale-free distribution. We chose a number close to an initial flat value of the scale-free network as a soft-threshold, based on which, the correlation matrix was converted into the adjacency matrix and the topological overlap matrix (TOM). We identified modules with hierarchical clustering (minModuleSize = 30) using dynamic tree-cutting method. Similarity between modules was measured and modules with a similarity greater than 0.75 were merged before extracting the final modules. The correlation between traits and different modules was calculated, and modules with a correlation exceeding 0.5 were identified as key modules.

From a Venn diagram generated in the online platform called OmicStudio (https://www.omicstudio.cn/), key genes existed in the intersection of key modules and the DEGs.

Functional Enrichment Analysis

Functional enrichment analysis is aimed at looking for biological functions of specific genes and biological processes that the genes participate in. The biological functions of key genes were assessed in terms of Gene Ontology (GO)¹⁷ function enrichment analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG)¹⁷ pathway analysis and reactome enrichment analysis. We uploaded 681 key genes to the online platform DAVID (https://david.ncifcrf.gov/) in order for the annotated information of GO and KEGG. Subsequently, a heatmap was generated using https://www.bioinformatics.com.cn, an online tool for data analysis and visualization, accessed last on 10 Nov 2023. Metascape (http://metascape.org) helped to make sure key signaling pathways. The results were visualized using networks and heat maps.

Gene set enrichment analysis (GSEA) version 4.3.2 software was employing for analyzing regulation of the associated pathways from the Canonical Pathways Gene Set (m2.cp.v2023.1.Mm.symbols.gmt) and the Gene Ontology Gene Set (m5.go.v2023.1.Mm.symbols.gmt). The random combination was set for 1000 times.

Identification by Machine Learning

LASSO, SVM and RF are utilized to sift the biomarkers related to clinical features, of which the former belongs to linear regression and the latter two belong to classification algorithm. In this study, OIR-related biomarkers were independently screened using three machine learning methods. In short, LASSO was executed using “glmnet” package in R to sift the gene signatures with the least classification error at the best lambda. SVM was conducted by “msvmRFE” and “e1071” package in R at 10-fold cross-validation to find the variables that achieve maximum accuracy; RF was carried out with “randomForest” package in R, and then gene features along with their contributions were analyzed through “ggplot2” R package. Ultimately, the overlapped genes screened via LASSO, SVM, and RF were regarded as reliable OIR-related biomarkers in ROP. To more accurately forecast patient incidence risk, nomograms were established utilizing hub genes. The Receiver Operating Characteristic (ROC) curves were plotted by the “pROC” package of R and calculated for performance evaluation to identify hub genes. Furthermore, we conducted the ROC curves of known ROP-related genes Hif1α and Gstp1, and compared the areas under the curve (AUC) values of hub genes and mixed genes.

Cell Culture

Human umbilical vein endothelial cells (HUVECs) are ones of cells most commonly used in researches related to ROP and angiogenesis.^18–20 HUVECs were purchased from the Wuhan Sunncell Biotechnology Company (SUNNCELL; Wuhan Sunncell Biotech, Inc; Hubei; China). These cells were maintained in HUVEC-SV40 cell culture medium (SUNNCELL; Wuhan Sunncell Biotech, Inc.; Hubei, China) with 5% fetal bovine serum (FBS), 1% P/S, 8 mM/L L-Glutamine, 100 mg/mL Sodium Pyruvate and endothelial cell growth factor at 37°C with 5% CO₂ for 24 h, which were taken for the normal control (NC) group. According to OIR cell lines,²¹ HUVECs were cultured in HUVEC-SV40 cell culture medium with 200 µg/mL CoCl₂ at 37°C with 5% CO₂ for 24 hours, which were regarded as the ROP group. We detected the mRNA level of Hif1α as a significant insight to confirm whether the OIR model was constructed successfully. Both cells of these two groups were continuously cultivated in HUVEC-SV40 cell culture medium at 37°C with 5% CO₂ for 24 h and used for the mRNA levels of hub genes.

Validation of Hub Gene by qRT-PCR

Hif1α and hub genes expression levels of ROP and NC groups were validated by qRT-PCR. The methods were displayed in earlier released articles.²² These RNAs were extracted by TRIzol^® (Invitrogen; Thermo Fisher Scientific, Inc; CA; USA) and reverse transcribed by cDNA Reverse Transcription Kit (EZBioscience, Roseville, MN, USA). Quantitative analysis was performed using SYBR Green I Master Mix (EZBioscience, Roseville, MN, USA) and a LightCycler® 480 Real-time PCR system (Roche, Basel, Switzerland). The qPCR primers were provided by BioTNT (Shanghai, China). GAPDH was used as an internal standard. Forward primer of Hif1α was as followed: GAACGTCGAAAAGAAAAGTCTCG; reverse primer was as followed: CCTTATCAAGATGCGAACTCACA. Forward primer of Plxnd1 was as followed: TGAGTCTGTTGTACGCTGTGA; reverse primer was as followed: GCCCCTTTAGTTGGAGGCT. Forward primer of Esm1 was as followed: ACAGCAGTGAGTGCAAAAGCA; reverse primer was as followed: GCGGTAGCAAGTTTCTCCCC. Forward primer of Cd248 was as followed: ATCGCAGCCAACTATCCAGAT; reverse primer was as followed: TTCCAGGCAAATGAGTGGTGG. The 2^–ΔΔct method was utilized for calculating relative gene expression levels.

Immune Cell Infiltration

Based on gene expression profiles, the CIBERSORT algorithm estimated the distribution of various immune cell types and the “mice” datasets from “Inference of immune cell composition on the expression profiles of mouse tissue”.²³ We analyzed the immune cells with a significance threshold of p < 0.05. The data of ROP was assayed to infer the relative proportions of 25 kinds of immune infiltrating cells, and Spearman correlation analysis was conducted on hub gene expression and immune cell content.

Exploration of ceRNA Network Analysis

CeRNA network analysis is used to explore the complex regulatory relationships between mRNAs, lncRNAs and miRNAs, which helps to understand how these RNA molecules interact with each other to regulate gene expression. TargetScan and miRWalk¹⁷ databases were used to identify possible miRNAs that target the hub genes, followed by the miRNet¹⁷ database to identify possible lncRNAs targeting the miRNAs. Data visualization was performed using Cytoscape.

Predication of Drug–Gene Interaction

The identification of interactions between existing or potential drugs and hub genes was conducted through Comparative Toxicogenomics Database (CTD; http://ctdbase.org/).¹⁷ Selected hub genes, regarded as potential pharmaceutical targets for ROP treatment, were imported into CTD to explore existing drugs and small organic compounds. Additionally, data visualization was performed through the Cytoscape software.

Statistical Analysis

The diagnostic value of the predictive model was assessed via univariate and multivariate logistic regression analyses. Two-sided p values were used in all statistical tests, and p (or FDR) < 0.05 was considered to indicate statistical significance. All analyses were performed using R software (version 4.0.2).

Results

Identification of DEGs Potentially Associated with ROP

We carried out bioinformatics analysis based on the annotated file of a dataset (GSE 130400). There were 25 retina samples that involved 13 samples from OIR groups and the others from the control groups. The “limma” R package was used to screen the differential expression OIR-related genes in GSE130400 with a p < 0.05 and a |log₂FC| ≥ 1. A total of 798 DEGs, including 190 down-regulated and 608 up-regulated genes, were identified between OIR and control samples and displayed in the volcano plot (Figure 2A). The heat map of differently expressed OIR-related genes was shown in Figure 2B.

Figure 2 Identification of DEGs potentially associated with ROP. (A) The volcano map of DEGs. Blue and red dots represent significantly lower and higher DEGs. (B) Heatmap of DEGs. Yellow and green in the group bar respectively represent the OIR and the control groups. Red and blue denote gene expressed at high and low levels in the samples, respectively.

WGCNA Analysis

WGCNA was applied to analyze genes of GSE 130400, and the co-expression modules and genes connected with OIR were found. The OIR-correlated coefficient was employed to cluster the samples, and a sample clustering tree was obtained (Figure 3A). We chose the soft-threshold 9 to achieve a scale-free topology (Figure 3B). The adjacency matrix was converted into a TOM matrix, which was utilized to show the similarity between nodes by considering the weighted correlation. Finally, 40 modules were identified through dynamic tree clipping (Figure 3C).

Figure 3 Construction and module analysis of WGCNA. (A) The samples clustering tree. (B) Network topology analysis under various soft threshold power settings. (C) Clustering tree of genes with different similarities based on topological overlap, and assigned module colors. (D) The heat map of different modules. The relationship between modules and ROP. (E) The scatter plot between the module membership in the dark olive green module and gene significance for cases. (F) The scatter plot between the module membership in the brown module and gene significance for cases. (G) The scatter plot between the module membership in the light steel blue1 module and gene significance for cases. (H) The Venn diagram consisting of two regular cycles, which severally represent 798 DEGs and genes of three modules in WGCNA.

Compared with other modules, we found three key modules (light steel blue1 module (cor = 0.92, p = 6e⁻¹¹), brown module (cor = −0.53, p = 0.007), and dark olive green module (cor = 0.5, p = 0.01), respectively) (Figure 3D). The former two modules showed positive correlation, while the latter exhibited a negative effect. There were 1147 (light steel blue1 module), 183 (dark olive green module), and 2381 (brown module) genes which were, respectively, found in the three modules (Figure 3E–G).

From the Venn diagram (Figure 3H), 681 key genes existed in the intersection of 798 DEGs and 3711 genes from three modules in WGCNA.

Enrichment Analysis

There were 681 key genes uploaded to DAVID and Metascape for GO and KEGG analyses. Based on p-value, GO functional enrichment analysis, respectively, revealed the top ten terms based on the gene count across biological process (BP), cellular component (CC) and molecular function (MF) categories (Figure 4A). Key genes were markedly enriched in positive regulation of cell proliferation, angiogenesis and positive regulation of apoptotic process and so on in the BP category. Enriched CC terms contained cytoplasm, cytosol and extracellular space terms and so on. In the MF category, genes were mainly enriched in protein binding, identical protein binding and actin binding activity and so on. KEGG pathway enrichment results showed that the mostly enriched pathways included PI3K-Akt signaling pathway, human papillomavirus infection and Kaposi sarcoma-associated herpesvirus infection and so on (Figure 4B). Metascape enrichment analysis results displayed that key genes were mainly enriched in functions and pathways like blood vessel development, positive regulation of cell migration, and extracellular matrix organization and so on (Figure 4C–E).

Figure 4 Functional enrichment analysis. (A) The bubble plots of GO enrichment analysis results. (B) The bubble plots of KEGG enrichment analysis results. (C) The heat map describing the top 20 most enriched items. (D) The network diagram showing related items. (E) The network diagram based on the p-value of different pathways. (F–I) GSEA analysis including reactome developmental biology (F), reactome nervous system development (G), reactome Rho GTPase cycle (H), reactome signaling by Rho GTPases miro GTPases and Rhobtb3 (I).

GSEA indicated that 681 key genes were mainly enriched in reactome developmental biology (Figure 4F), reactome nervous system development (Figure 4G), reactome Rho GTPases cycle (Figure 4H), reactome signaling by Rho GTPases miro GTPases and Rhobtb3 (Figure 4I).

Identification and Validation of Hub Genes

For further selection of the hub genes with a significantly characteristic value of classifying the ROP groups and control groups, three different algorithms (LASSO, SVM-RFE, and RF) were used on the 681 key genes. Using a tenfold cross-validation framework, we identified the optimal λ value of 0.03926209 in the LASSO analysis by minimizing the partial likelihood deviation (Figure 5A). Afterward, the results of LASSO model revealed that 25 genes were associated with the occurrence of ROP (Figure 5B). Meanwhile, the feature vectors generated by SVM were removed to find the best variables and 9 genes were identified (Figure 5C). RF algorithm was also used to rank the 681 key genes in the lights of the variable importance of each gene, and the genes with the MeanDecreaseGini > 2 were selected, which found 30 genes (Figure 5D).

Figure 5 Exploration and validation of hub genes by machine learning. (A) LASSO regression of the 25 genes. (B) Cross-validation for tuning the parameter selection in the LASSO regression. (C) The root mean square error was calculated from 10-fold CV and verified the results of SVM-RFE. (D) The mean decrease Gini coefficients of genes in the RF classifier. (E) Venn diagram by identifying genes common to results of LASSO, SVM and RF. (F–H) Nomogram model diagram based on the expression levels of three hub genes in groups of ROP and control. (I) The ROC curves. Red, blue, green, black and purple separately represent three hub genes Esm1, Plxnd1 and Cd248, and two ROP-related genes Hif1α and Gstp1. (J) The ROC curves. Red represents the combination of Plxnd1, Esm1 and Cd248. Blue represents the combination of Hif1α and Gstp1. (K) The bar chart of Hif1α mRNA levels in ROP and NC groups, **p < 0.01 indicated statistical significance compared with the control. (L–N) The bar charts of Esm1, Plxnd1 and Cd248 mRNA levels in ROP and NC groups, ***p < 0.001, ****p < 0.0001 indicated statistical significance compared with the control.

From the Venn diagram (Figure 5E), we found three hub genes including Plxnd1, Esm1 and Cd248 that existed at the overlapping position of LASSO, SVM and RF. The nomogram model diagram showed that each hub gene was a high-risk factor of developing ROP (Figure 5F–H) and Esm1 had greater significance in the risk of developing ROP, while compared to the other two hub genes. ROC analysis indicated that all hub genes had area under the curve (AUC) values > 0.9, with Esm1 having the largest AUC value (AUC, 0.987) and Plxnd1 the smallest (AUC, 0.91), and HIF1α (AUC, 0.609) and Gstp1 (AUC, 0.59) had smaller values compared with three hub genes (Figure 5I). In addition, the combined value of three hub genes (AUC, 1) was the highest and remarkably more than the combined value of two ROP-related genes (AUC, 0.654) (Figure 5J). qRT-PCR showed that the expression of hub genes in the ROP group was significantly higher than that in the NC group (Figure 5K–N).

Immune Infiltration Analysis

We found that the function and pathway analysis of genes intersected in WGCNA and 798 DEGs in ROP were tightly associated with inflammatory and immune processes. CIBERSORT was performed to analyze the abundances of immune cells of all samples in different groups, and the results were shown in Figure 6A. At the side of case tissue, ROP had a higher abundance of M0 macrophage, and M2 macrophage was significantly lower (Figure 6B). Additionally, the correlation analysis demonstrated that Plxnd1 showed significantly positive correlation to Th2 cells (R = 0.437, p = 0.029) and NK resting (R = 0.405, p = 0.044), and that Esm1 was positively associated with Th2 cells (R = 0.668, p < 0.001), NK resting (R = 0.622, p < 0.001) and M0 macrophage (R = 0.411, p = 0.041), and finally, that Cd248 had greatly positive association with Th2 cells (R = 0.622, p < 0.001) and M0 macrophage (R = 0.477, p = 0.016) (Figure 6C, F and G). The scatter plots between hub genes and immune cells were shown in Figure 6D, E and H–L.

Figure 6 Immune infiltration analysis and the correlation between the hub gene and the immune cells. (A) Heat map displaying the ratio of 25 immune cells in ROP and control. (B) Boxplot displaying the immune cell proportions between ROP and control groups. (C, F and G) The lollipop plots showing the correlation between immune cells and hub genes. (D, E and H–L) Scatter plots of hub genes and related immune cells, including Plxnd1 and Th2 cells (D), Plxnd1 and NK resting (E), Esm1 and Th2 cells (H) Esm1 and NK resting (I), Esm1 and M0 macrophage (J), Cd248 and Th2 cells (K), Cd248 and M0 macrophage (L).

Further investigation of the connection between hub genes and immune checkpoint genes indicated Plxnd1 had a strongly positive correlation (R = 0.608, p < 0.001) (Figure 7A) Esm1 (R = 0.438, p = 0.029) and Cd248 (R = 0.489, p = 0.013) had lowly positive correlation with H19 immune checkpoint gene (Figure 7B and C). Furthermore, the association between hub genes and immune chemokine genes showed that Plxnd1 had positive relevance to H19 (R = 0.608, p < 0.001), C2 (R = 0.460, p < 0.002), F7 (R = 0.449, p = 0.024) and F5 (R = 0.403, p = 0.046), Esm1 was positively related to F5 (R = 0.598, p = 0.002), F7 (R = 0.546, p = 0.005) and H19 (R = 0.438, p = 0.029), and Cd248 was connected with H19 (R = 0.489, p = 0.013) (Figure 7D–F). The outcomes demonstrated a notable positive association between the expression of hub genes and multiple immune characteristics.

Figure 7 The correlation of hub genes and immune checkpoints and chemokines. (A–C) The lollipop plots showing the correlation between immune checkpoint genes and hub genes. (D–F) The lollipop plots showing the correlation between immune chemokine genes and hub genes.

The mRNA–miRNA–lncRNA ceRNA Network of ROP

To establish the ceRNA network, we searched for miRNA–mRNA targets and obtained 33 miRNAs of Esm1, 1 miRNA of Cd248 and 3 miRNAs of Plxnd1, results of which were in the intersection of TargetScan and miRWalk database, which was shown in Figure 8A. Then, the potential lncRNA that targeted the 36 miRNAs was predicted through miRNet online databases, and 17 lncRNAs targeting 5 miRNAs (including mmu-miR-679-3p, mmu-miR-293-5p, mmu-miR-3059-5p, mmu-miR-669d-3p and mmu-miR-1306-5p) were obtained.

Figure 8 CeRNA network analysis of hub genes and drug–gene interaction analysis. (A) The linkage map of 3 hub genes and their corresponding miRNAs and lncRNAs, in which blue represents miRNA, green represents lncRNA, and yellow represents hub gene. (B) The linkage map of 3 hub genes and potential target drugs, in which yellow represents hub gene, green represents common target drugs for three hub genes, blue represents potential drugs targeting double hub genes, and dark blue represents prospective drugs targeting simple hub gene.

Drug–Gene Interaction of Hub Genes

Possible therapeutic drugs for targeting three hub genes were sifted from the CTD database. The drug–gene interaction network of the above genes was exhibited in Figure 8B, which involved in 71 potential therapeutic drugs for Esm1, 29 latent drugs for Cd248 and 44 prospective drugs for Plxnd1. There were 14 common potential target drugs for these hub genes including sodium arsenite, CGP 52608, bis(4-hydroxyphenyl) sulfone, Doxorubicin, ect.

Discussion

ROP is a disease that affects the development of the eyes of premature infants. It develops in the first few months after birth and is usually associated with low birth weight and underdeveloped retina in premature infants. VEGF is the angiogenic factor that has been researched the most. In the vitreous of ROP infants, there is a significant rise in VEGF levels.^24,25 The factor is generally regarded as the key element of angiogenesis in ROP, whether physiological or pathological, and in various other neovascular conditions.²⁶ In addition, inflammation response also is thought to have a major influence on the progression of ROP.²⁷ Four drugs, administered intravitreally as anti-VEGF agents, have been applied for ROP treatment: ranibizumab, bevacizumab, aflibercept and conbercept. However, the effect of using medication alone is not ideal, and even with surgery, there will be various complications such as infection and conjunctival congestion. Moreover, these drugs have many adverse effects, such as a high incidence of complications and myopia²⁸ and systemic VEGF suppression.^29,30 Therefore, it is necessary to explore better therapeutic drugs. However, pathogenesis of genes in ROP has not been fully elucidated, which poses great difficulties for drug research and development. Thus, our research is aimed at screening out hub genes associated with ROP and find out potential therapeutic drugs targeting them. In this study, based on the expression profiles of ROP, we carried out differently expressed analysis and WGCNA for key genes. Then, hub genes were explored by multiple machine learning methods. Furthermore, hub genes were validated by the nomograms, ROC curves and qRT-PCR. Immune infiltration analysis deeply explored the relationships between hub genes, immune response and ROP occurrence. Finally, through the ceRNA and drug–gene interaction networks, the gene regulation and potential therapeutic drugs were screened for ROP hub genes.

It is difficult to obtain clinical retinal tissues for ROP researches, so the mouse model is usually used for the mechanism study of ROP. It has been elaborated in many influential researches in recent years^31,32 that OIR mouse model is a well-established tool for studying ROP and abnormal angiogenesis. In the OIR mouse model, P12-P17 is characterized by pathological angiogenesis at the junction of avascular and vascularized retina.³³ At P 17, the most abundant pathological vascularization occurred.¹⁵ As a result, we analyzed the samples collected at P17 (including the 13 OIR samples and the 12 control samples) in GSE130400, and screened 681 key genes via differently expression analysis and the weighted gene co-expression network. Then, three hub genes (Esm1, Plxnd1 and Cd248) were obtained by machine learning methods (LASSO, SVM-RFE and RF). WGCNA exhibited that they all existed in the light steel blue1 module which had a close relationship with ROP (module eigengene > 0.9). Module Membership of hub genes were more than 0.9, which suggested that the hub genes played key roles in the characteristics of the light steel blue1 module. The gene significance of each hub gene was greater than 0.9, displaying that three hub genes had implications for the ROP’s development. In addition, based on the results of the expression profiles, these hub genes showed an important upregulation in the ROP group compared to the control group (|log₂FC| > 4) and expression differences of them were greatly higher than most other genes in differently expressed genes. qRT-PCR results demonstrated that hub genes expression levels were higher in ROP groups compared with NC groups, which showed three hub genes had relatively great effects with ROP. Furthermore, we validated the risk that these hub genes induced ROP. Esm1, Plxnd1 and Cd248 had areas under the curve (AUC) values >0.9 that were more than Hif1α and Gstp1, displaying ideal predictive performance for ROP. The nomogram models, which relied on the expression levels of hub genes, showed good performance in significantly differentiating the morbidity of ROP, implying that Esm1, Plxnd1 and Cd248 were key factors in ROP, particularly Esm1.

In the enrichment analysis based on 681 key genes, we identified angiogenesis, blood vessel development, cell migration, positive regulation of endothelial cell proliferation and reactome Rho GTPases cycle and so on. The retinal vascular diseases are marked by aberrant angiogenesis accompanying the dysregulated proliferation and migration of endothelial cells (ECs).^34,35 Therefore, angiogenesis, blood vessel development, positive regulation of endothelial cell proliferation, cell migration and regulation are involved in the occurrence and development of ROP. Furthermore, members of the Rho GTPase family, such as RhoA, Cdc42 and Rac1, are involved in signal transduction and cellular behaviors during retinal blood vessel growth and vessel regression.³⁶ They affect retinal vessel development by regulating mechanisms such as cytoskeletal reshaping, cell migration, adhesion and contraction.³⁷ Thus, reactome Rho GTPases cycle has close association with ROP. In general, the results of enrichment analysis indicated that 681 key genes mainly took part in angiogenesis, blood vessel development, cell migration and so on.

The enrichment results suggested that Esm1, Plxnd1 and Cd248 went hand in hand with ROP. Plxnd1 is part of the Plexin transmembrane protein family and mediates semaphorin signalling. PLXNs contain an extracellular SEMA domain, one membrane spanning region and a cytoplasmic domain that has a GTPase activating protein (GAP) domain and a Rho GTPase binding domain insert.^38,39 When the extracellular domain binds to a SEMA and the intracellular domain binds to a Rho GTPase concurrently, the GAP activity of Plxnd1 is activated, involved in reactome Rho GTPases cycle to modulate integrin—mediated cell adhesion, cytoskeletal dynamics and cell migration. Esm1 is a dermatan sulfate proteoglycan that is expressed and secreted mostly by endothelial cells.⁴⁰ This gene holds a place of importance in angiogenesis. Esm1 enhances the promigratory and proliferative effects of VEGF on endothelial cells.⁴¹ Esm1 positively modulates VEGF-A165 signaling by binding to fibronectin and preventing the sequestration of VEGF-A165 by fibronectin, which increases the bioavailability of VEGF-A165 to regulate angiogenesis.⁴² Esm1 can also interact with multiple angiogenic molecules like HGF/SF and FGF to form a positive feedback loop, which further induces the proliferation and tubular structure formation of vascular endothelial cells.⁴³ As a result, Esm1 plays a crucial role in angiogenesis, blood vessel development, blood vessel morphogenesis, and routing angiogenesis. Cd248 is a transmembrane glycoprotein that belongs to the C-type lectin-like receptor family. Cd248 can strengthen the interaction of ITGB1 and extracellular matrix proteins like FN and CYR61, thus activating FAK-paxillin pathway to facilitate cell migration and metastasis.⁴⁴ Thus, Cd248 fulfills the function of cell migration. All the results showed that hub genes Plxnd1, Esm1 and Cd248 played important roles of the biological processes pertinent to ROP.

Some factors make contributions to the pathogenesis of ROP; nevertheless, inflammatory processes stand as major one.⁴⁵ To thoroughly investigate the changes regarding immune cells in ROP, we conducted an immune infiltration analysis. In our study, ROP tissue had a greater percentage of M0 macrophage, but relatively smaller ones of M2 macrophage. Macrophages in a resting state (M0), originating from bone marrow, are typically seen as the precursors of polarized macrophages. The mainstream is that M0 merely represents a resting state of macrophages, without a specific function before their polarization, and M1 macrophages and M2 macrophages are derived from M0.⁴⁶ M1 macrophages are seen as phagocytic and pro-inflammatory, secreting pro-inflammatory cytokines, such as TNFα and IL-6, which can lead to damage and inflammatory response of retinal endothelial cells. Some substances secreted by M1 macrophage, like rhFABP4, could promote vascular endothelial growth factor α (VEGFα) protein expression to significantly enhance tube formation.⁴⁷ M2 macrophages release cytokines that promote the proliferation of contiguous cells and tissue repair. M2 macrophages have anti-inflammatory effects, and the decrease in the proportion of M2 macrophages may facilitate pro-inflammatory of M1 macrophages. The previous research showed that M2 macrophage polarization could be mitigated by IL-17, which increased retinal neovascularization in ROP.²³ Our results verified that immune response is closely related to ROP, further revealing that the reduction of M2 macrophages may promote inflammation and lead to angiogenesis in ROP.

Then, we explored the relationship between three hub genes and immune cells. In our study, hub genes Plxnd1, Esm1 and Cd248 had positive correlation with Th2 cells, M0 macrophage and NK resting cells. The Th2 subset is one of the subtypes of helper T cells, which has mutual regulation with Th1 subset, and generates cytokines like IL-4 and IL-5 that facilitate B cell proliferation and differentiation, which is related to humoral-type immune responses.^48–50 What’s more, Th2 cells can release pro-angiogenic factors⁵¹ to promote vessel formations. NK resting cell is one kind of NK cell state. When NK cell is in a resting state, its proangiogenic activity is enhanced.⁵² Our research suggested that Plxnd1, Esm1 and Cd248 had positive association with the upregulation in Th2 cells, M0 macrophage and NK resting cells, which might promote inflammation and angiogenesis to take part in the development of ROP. It was confirmed that the autocrine Sema4A-PlexinD1 axis was a negative regulator of human Th1 differentiation and played a key role in Th2 skewing.⁵³ In addition, changes in activation and resting states of NK cells are associated with the expression of Cd248.⁵⁴ However, the mechanisms by which hub genes regulate Th2 cells, M0 macrophage and NK resting cells are still unclear and need to be studied.

Furthermore, we explored the relationship between three hub genes and immune checkpoints and chemokines. We found that Plxnd1, Esm1 and Cd248 positively correlated with lncRNA H19 (H19), CXCR5 (F5), CCL27 (F7) and CCL2 (C2). It was demonstrated that H19 might take part in several cellular processes such as inflammation and apoptosis via enhancement of the expression of X-box binding protein 1.⁵⁵ Inflammation is a major element giving an impetus to ROP, and H19 may be involved in ROP through pro-inflammatory effects. In the bargain, H19 regulated the production of phosphorylated ERK1/2 proteins to repress the MAPK–ERK1/2 signaling pathway, which adjusts fibroblast proliferation to assist angiogenesis.^56,57 CXCR5 is the sole receptor of CXCL13.^58,59 CXCL13/CXCR5 triggers chemotaxis, proliferation, antibody secretion, and intracellular signal transduction in activated B cells.^60,61 They also promote immunosuppressive cytokine TGF-β1 production, which may impair NK cell cytotoxic functions and cause the increased proportion of NK resting cells to promote vessel creation.⁵² In addition, CXCL13/CXCR5 drives VEGF expression to induce angiogenesis.⁶⁰ CCL27 has been connection with the recruitment and infiltration of immune cells via ligation with CCR10. Besides, CCR10/CCL27-CCL28 results in endothelial cells migration and angiogenesis.⁶² CCL2, also referred to as MCP-1, is a pro-inflammatory chemokine encoded by the Ccl2 gene, and it plays a role in recruiting monocyte-derived macrophages to inflamed sites.⁶³ In our results, Plxnd1, Esm1 and Cd248 had certain relevance with the up-regulated expression of H19, F5, F7 and C2, which may make contributions to the occurrence and development of ROP. It is suggested that CXCR5 is expressed on the surface of B cells.⁶⁴ Plxnd1 influences activated B cell migration toward CXCL13 chemokine,⁶⁵ implying that Plxnd1 is an important foundation of combining CXCL13 and CXCR5. Furthermore, the overexpression of H19 increases the level of Esm1 mRNA and protein,⁶⁶ which hints Esm1 has high relationship with H19. Nevertheless, mechanisms of Plxnd1, Esm1 and Cd248 regulating H19, F5, F7 and C2 in ROP are not yet clear and remain to be studied.

In sum, we carried out a comprehensive and systematic bioinformatics analysis and discovered 3 hub genes, Plxnd1, Esm1 and Cd248, which are mainly involved in the angiogenesis, blood vessel development, cell migration, positive regulation of endothelial cell proliferation and reactome Rho GTPases cycle. Plxnd1 hastens cell migration by reactome Rho GTPases cycle. Esm1 promotes angiogenesis by increasing the bioavailability of VEGF-A165 and induces the proliferation and tubular structure formation of vascular endothelial cells through interacting with molecules of vascularity, like HGF/SF. Cd248 can stimulate cell migration via activating FAK-paxillin pathway. Furthermore, our study found that Plxnd1, Esm1 and Cd248 had close contact with inflammation. Plxnd1, combining with Sema4A, can furtherance an increase in the proportion of Th2 cells to redound to angiogenesis. Cd248 can regulate angiogenesis by affecting the number of NK resting cells. Therefore, hub genes promoted vascularization by immune cells Th2 cells and NK resting cells to participate in ROP. Additionally, three hub genes promote the upregulation of checkpoints and chemokines H19, F5, F7 and C2 to stimulate inflammation and play key roles in furtherance of cell proliferation and migration in ROP. Hence, Plxnd1, Esm1 and Cd248 might be possible key biomarkers in the progression of ROP. Furthermore, our research demonstrated that mRNAs of Plxnd1 and Cd248 could be inhibited by miRNA mmu-miR-7050-3p, and meanwhile Esm1 could be regulated by Meg3/Malat1-mmu-miR-679-3p/mmu-miR-669d-3p axis. Research showed that lncRNA Meg3 could regulate apoptosis and inflammation via NF-κ B signaling pathway.⁶⁷ Malat1 could promote endothelial cells proliferation and induced retina angiogenesis in the neonates.⁶⁸ Therefore, the result of the ceRNA network further revealed that Plxnd1, Esm1 and Cd248 played significant roles in the development of ROP. In summary, our research revealed three new hub genes closely related to ROP, which would be of great significance for further the elucidation of ROP mechanism, the study of new diagnostic methods and the development of novel drugs.

Finally, we identified 14 potential therapeutic drugs targeting three hub genes through CTD, and 4 of them were approved for clinical treatment. CTD provides high-quality and valuable drugs based on existing literature and studies and has been widely used for a great deal of drug prediction researches based on targets.^69–71 Among them, sodium arsenite was previously widely used for the treatment of certain types of leukemia, such as Acute Promyelocytic Leukemia. Tretinoin is utilized to treat acne vulgaris and photodamage. Add to that, Doxorubicin was utilized against malignancies including solid tumours, transplantable leukemias and lymphomas,⁷² the mechanism of which was to reduce the expression of VEGF and restrain the proliferation and migration of endothelial cells.^73,74 Esm1 functioned as an important factor in the activity of VEGF, which will influence the proliferation and migration of endothelial cells.⁴² It meant that Doxorubicin was likely to inhibit the effect of Esm1 to hinder the function of VEGF and proliferation and migration of endothelial cells and further to impede the process of ROP. In conclusion, we found out these potential therapeutic drugs for ROP treatment based on the three hub genes. Although there was still a lack of researches of further function validation for 3 hub genes and drug treatments for ROP based on three hub genes, our study provided a novel viewpoint for the nosogenesis and treatment progress of ROP.

Conclusion

In this study, we obtained three hub genes Plxnd1, Esm1 and Cd248 related to ROP by WGCNA analysis, machine learning methods, and then validated hub genes by the nomograms, ROC curves and qRT-PCR. Plxnd1, Esm1 and Cd248 took part in endothelial cell proliferation and migration and angiogenesis, and evoked inflammatory responses through upregulation of immune cells TTh2 cells, M0 macrophage and NK resting cells and a higher expression of chemokines H19, F5, F7 and C2 to make effects in ROP progression. Obtainment of clinical retina is inaccessible, so data source of this study was based on animal models. However, this research laid down the foundation for further clarifying pathogenesis of ROP. In summary, this work revealed that Plxnd1, Esm1 and Cd248 could be diagnostic markers and potential therapeutic targets for ROP. However, further research is needed to explore the molecular mechanism by which these three hub genes promote ROP and to develop drugs targeting these three genes.