Screening and identification of genes related to ferroptosis in keratoconus

Corneal keratoconus (KC) is a dilated (ectatic) corneal disease characterized by a central thinning of the cornea, which causes protrusion into a conical shape that seriously affects vision. However, due to the complex etiology of keratoconus, its entire mechanism remains unclear and there is no mechanism-directed treatment method. Ferroptosis is a novel programmed cell death mechanism related to lipid peroxidation, stress, and amino acid metabolism, which plays a crucial role in various diseases. This study aimed to explore the relationship between keratoconus and ferroptosis, to provide new insights into the mechanism of keratoconus development, and potential treatment options based on further elucidation of this mechanism. The corresponding mRNA microarray expression matrix data of KC patients were obtained from GEO database (GSE204791). Weighted co-expression network analysis (WGCNA) and support vector machine recursive feature elimination (SVM-RFE) were selected to screen hub genes, which were overlapped with ferroptosis genes (FRGs) from FerrDb. GO and GSEA were performed to analyze differential pathways, ssGSEA was used to determine immune status, and then, feasible drugs were predicted by gene-drug network. Additionally, we predicted the miRNA and IncRNA of hub genes to identify the underlying mechanism of disease so as to predict treatment for the disease. The epithelial transcriptome from keratoconus tissue mRNA microarray data (GSE204791) was extracted for the main analysis, including eight epithelial cells and eight epithelial control cells. The differential genes that were overlapped by WGCAN, SVM-RFE and FRGs were mainly related to oxidative stress, immune regulation, cellular inflammation, and metal ion transport. Through further analysis, aldo–keto reductase family 1 member C3 (AKR1C3) was selected, and negatively correlated with mature CD56 natural killer (NK) cells and macrophages. Then, gene-drug interaction network analysis and miRNA prediction were performed through the website. It was concluded that four immune-related drugs (INDOMETHACIN, DAUNORUBICIN, DOXORUBICIN, DOCETAXEL) and a miRNA (has-miR-184) were screened to predict potential drugs and targets for disease treatment. To our knowledge, this was the first report of KC being associated with ferroptosis and prompted search for differential genes to predict drug targets of gene immunotherapy. Our findings provided insight and a solid basis for the analysis and treatment of KC.


Materials and methods
Microarray data acquisition.We searched the keywords ("keratoconus" AND "Homo sapiens" [porgn:__ txid9606]) and "Expression profiling by array" using the GEO database (https:// www.ncbi.nlm.nih.gov/ geo/) to obtain mRNA Microarray Gene Expression Matrix of Keratoconus.We obtained GEO profiles (GSE204791) which included epithelial cells (8 EN, 8 EKC) and stromal cells (7 SN, 8 SKC) of keratoconus.The GSE204791 dataset matrix is a public database based on platform document GPL21185 (Agilent-072363 SurePrint G3 Human GE v3 8 × 60 K Microarray 039494 [Probe Name Version]), the patients involved in the database had previously given ethical approval, so our study did not involve ethical issues and other conflicts of interest.The workflow is shown in Fig. 1.

Data stability verification.
The dist function was used to calculate the distance between each data category; the R package (ggplot) was executed for visual analysis of the data to make a clustering dendrogram.In addition, we performed dimensionality reduction on the data, and used the ggplot (3.3.6)package to apply visualization to the PCA analysis of the data.
Gene set enrichment analysis.The gene set enrichment analysis (GSEA) was used to evaluate the distribution trend of the defined gene set in the gene table ranked relative to the phenotype, to judge its contribution to the phenotype.Gene expression data in GSE204791 was analyzed with the use of the c2.cp (C2.cp.v2023.1hs.Symbols.gmt) collection, and the Molecular Signature Database (MSigDB) served as a reference gene set.Gene set data were included in the GSEA software (version 4.3.2),and item significant enrichment criteria: p value < 0.05.

Construction of weighted gene co-expression network analysis. Weighted Gene Co-Expression
Network Analysis was performed using R package (WGCNA).The "limma" package normalized the data and used Pearson to cluster the samples and average linkage analysis, delete small fluctuations and missing values; the adjacency matrix was transformed into a TOM matrix and genes were then clustered.The screened groups were represented by different colors, and the gene set (p < 0.05, correlation ≥ 0.5) were overlapped with the ferroptosis database (http:// www.zhoun an.org/ ferrdb/ curre nt/); there were 728 genes involved in ferroptosis, including driver genes, suppressor genes and marker genes.

Gene enrichment analysis.
To obtain the function and biological process of ferroptosis-related differential mRNA, the data were annotated by gene ontology (GO) using the "go plot", "cluster profiler", and "circlize" package in R version (4.3.2), and then the ggplot2 package was used for visualization.

Support vector machine-recursive feature elimination (SVM-RFE). Support vector machine
(SVM) is a monitoring machine learning technique widely used in classification and regression analysis, and the best genes were selected from the metadata cohort using the RFE algorithm.In conclusion, to identify gene sets Gene expression analysis.The gene expression matrix was extracted from the transcriptome samples of GSE204791 using R (4.3.2).And then, the ggpubr package was used to visualize the gene expression levels of AKR1C3 and SLC7A11.A significant difference was considered when p < 0.05.The samples for data analysis are the transcriptome samples based on GSE204791, and the dataset matrix is a public database based on the platform file GPL21185.Therefore, the use and analysis of this data do not involve medical ethical review and other conflicts of interest.

Gene-drug analysis network.
The gene-drug network was predicted using the website (https:// dgidb.genome.wustl.edu/).The predicted results were based on drugs that target genes reported in previous studies, and a total of 10 different drugs were predicted, including 4 immunomodulatory drugs.

Results
Significant signaling pathways were screened by GSEA.Data were obtained from keratoconus epithelial cells (KEC) and epithelial basal cells (EN) in GSE204791, and SKC and SN were used for partial validation.The results of the box-plot (Supplementary Fig. S1A) showed that the black lines on the medial axis of all samples were roughly at the same level, and, therefore, the standardization degree of the samples was satisfactory.The PCA results also confirmed the replicability of the data (Supplementary Fig. S1B).Similarly, we obtained the same results in SKC and SN (Supplementary Fig. S2A and B).In sample processing, we deleted the stromal cell sample GSM6193943, because the sample deviated greatly and did not have satisfactory repeatability characteristics.Next, we performed a cluster tree analysis on the 14 samples using the class-averaging method of hierarchical clustering-Euclidean clustering, and the clustering dendrogram confirmed the clustering relationship between each group (Supplementary Fig. S1C).Our purpose was to illustrate that GSE204791 sample www.nature.com/scientificreports/data were reliable and could be used for further analysis.There were 25,964 genes displayed in the KEC and EN samples, and 362 genes with significant differences (p value < 0.05, |Log 2 FC (Fold Change)|> 0.58), of which 156 were significantly up-regulated (indicated in red), and 206 were significantly down-regulated (indicated in blue) (Supplementary Fig. S1D).

GSEA analysis.
The gene expression data in GSE204791 were analyzed holistically by GSEA using the C2 (CP) gene sets (MSigDB).As shown in Fig. 2A, the most abundant item identified by GSEA was Cytokine Cytokine Receptor Interaction (p.adj < 0.001), and the enrichment results showed down-regulation.In addition, we found that the gene sets were also gathered in Signaling By interleukin (p.adj = 0.002) (Fig. 2B).Partial GSEA results are provided in Supplementary Table 1.Based on GSEA analysis, the ferroptosis-related gene set was also enriched significantly in the EKC (p.adj = 0.007), and mostly downregulated (Fig. 2C).Similarly, we found the same results in SKC and SN (Supplementary Fig. S3A and B); the ferroptosis-related gene set was enriched significantly, and also mostly downregulated (Supplementary Figure S3C).
WGCAN screened the hub genes.The hub genes were screened for EKC and EN.We normalized the data and Pearson's correlation coefficient to cluster the samples.WGCNA analysis generally included genes with a false discovery rate < 0.05 and log 2 FC (Fold Change) ≥ 0.5.In our study, the genes with low variability and samples with missing values were identified and removed; all samples were checked and no outlier samples were deleted (Supplementary Fig. S4A).We used the PickSoft-Threshold function in the WGCNA package for analysis.When the soft threshold was 19, the topological relationship of the scale-free network with R2 equal to 0.9 was satisfied, and it had the best connectivity (Supplementary Fig. S4B and C).We constructed gene networks and identified modules using the one-step network building function of the WGCNA R package.The adjacency matrix was converted into TOM matrix, then the number of genes in the dynamic clipping module was set to 60 and the depth segmentation was set to 2 (which implies a medium sensitivity) to construct the WGCNA network (Supplementary Fig. S4D).Finally, 10 modules were successfully constructed, which are represented in 10 different colors (Supplementary Fig. S4E).It is noteworthy that the ME brown (r = 0.82, p = 1e−04), ME black (r = − 0.57, p = 0.02), and ME magenta (r = − 0.65, p = 0.007) modules were significantly correlated with the disease state and may have played an important role in keratoconus (Supplementary Fig. S4F); they were used for further analysis.Furthermore, we also confirmed there was a highly significant correlation between Module membership (Mm) and Gene significance (GS) in these modules (black module: Cor = 0.51, p = 1.4e−45 (Fig. 3A); brown module: Cor = 0.7, p = 1.2e−148 (Fig. 3B); magenta module: Cor = 0.5, p = 9.8-11e (Fig. 3C)).We overlaid these 4 modules with FRGs and found that 8 FRGs (LCE2C, NOS2, AKR1C3, CAV1, MMD, LINC00618, SNCA, SLC7A11) were closely related to the disease state of keratoconus (Fig. 3D).

Relationship between core genes and immune microenvironment.
In previous studies, from the perspective of data integrity analysis, KC was affected by inflammatory immune regulation, especially the Signaling by Interleukins.We performed Single Sample Gene Set Enrichment Analysis (ssGSEA) on the data by www.nature.com/scientificreports/analyzing the association between samples and immune cells as a whole.Immune scoring of the samples was performed, and then visual analysis of the immune cell subsets of EKC and EN were performed to explore the different immune cell subsets in the samples.The results showed that keratoconus epithelial cells significantly affected immune subtype cells such as B cells, T cells, and natural killer cells (Fig. 7A,B).The immune cells regulated by AKR1C3 included CD56dim natural killer cells and macrophages (Fig. 7C).Similarly, immunophenotypic differences were observed between SN and SKC samples (Supplementary Fig. S7).These results suggested that KC may be affected by immune regulation.
Gene-drug therapy and mRNA target prediction.The Drug and Gene Interaction website was used to predict drugs that interacted with AKR1C3 (https:// dgidb.genome.wustl.edu/).The use of these drugs was based on previously reported studies that targeted AKR1C3 for gene therapy.After systematic screening, 10 drugs were identified (DAUNORUBICN, EXEMESTANE, TESTOSTERONE, DOXORUBICIN, BACCHARIN, INDOMETHACIN, DOCETAXEL, CHEMBL1682202, CHEMBL1682201, CHEMBL1682200).Among them, INDOMETHACIN, DAUNORUBICIN, DOXORUBICIN and DOCETAXEL exerted therapeutic effects by immunotherapeutically targeting AKR1C3 (Fig. 8A).Next, we predicted miRNA binding to AKR1C3 through miRNA databases (miRanda, miRDB, TargetScan), and the prediction results needed to satisfy the three databases at the same time.Five miRNAs satisfied the conditions (hsa-miR-3136-5p, hsa-miR-4282, hsa-miR-184, hsa-miR-3176, hsa-miR-379-5p) (Fig. 8B).We found that hsa-miR-184 was directly correlated with AKR1C3, miR-184 was the most abundant miRNA expressed in the cornea, which is related to the regulation of protein levels in the cornea.Next, we predicted the IncRNAs that bind to has-miR-184: RP11-830F9.6 and MLLT4-AS1.Therefore, we speculated that hsa-miR-184 may be an effective target for the treatment of KC and exert its effect through immunization.

Discussion
KC is characterized by corneal thinning, central protrusion, and irregular astigmatism, resulting in serious compromise of visual acuity.The most widely accepted etiologies of keratoconus have been ascribed to genetic and environmental factors; it has been stated that KC is well known as an independent symptom that does not involve any systemic or ocular disease 11 , and that KC is a noninflammatory corneal disease 29 .Treatments of keratoconus have principally relied on corneal transplantation 30 and corneal cross-linking 31 .However, a growing body of evidence suggests that the progression of KC involves the up-regulation of proinflammatory mediators and the involvement of Toll-like 2/4 receptors, which are influenced by the immune system 32,33 .The ocular surface of patients with KC has been clearly shown to demonstrate immune-inflammatory features 13 .Moreover, large amounts of superoxide, hydrogen peroxide and hydroxyl radicals have been reported on the ocular surface of keratoconus patients 34,35 .The oxidative stress response induced by environmental changes in KC is the direct cause of increased ROS levels and decreased antioxidant levels in keratoconus cells, which ultimately are responsible for extracellular matrix degradation and subsequent corneal stromal thinning 16,36 .
Ferroptosis is a unique, iron-dependent, non-apoptotic form of cell death.The accumulation of iron-dependent lipid peroxidation is mainly triggered by the Fenton reaction and lipoxygenases.Ferroptosis can affect the inflammatory response of cells through immunogenicity, and inhibition of ferroptosis can effectively suppress the inflammatory response; ferroptosis inhibitors have been shown to play important roles in several diseases through their anti-inflammatory effects 37 .Evidence of inflammation in ferroptosis has been reported, and the infiltration of neutrophils and the expression of pro-inflammatory cytokines were inhibited by Fer-1 38 .It is logical, therefore, that ferroptosis may play an important role in KC, but prior research supporting this concept has been lacking.In the present study, we elucidated the relationship between KC, ferroptosis and the immune microenvironment for the first time, which could reasonably provide new clues for the molecular mechanism of KC.In this study, we used the data from GSE204791 to identify mRNAs associated with ferroptosis, extracted KEC and EN data for analysis, using the SKC and SN for study verification.Then we found that cell-receptor interaction, inflammation, and ferroptosis genes were significantly enriched when GSEA was used to analyze RNA-seq data.Notably, certain immune and inflammatory factors have previously been detected on the ocular surface of SK patients, as well as products of mtDNA damage in the KC cornea 39 .In our study, excessive ROS and RNS, disruption of antioxidant machinery, and mitochondrial dysfunction were found in KCS, suggesting that increased oxidative stress played an important role in the pathogenesis of KCS 40 .Therefore, we postulated that keratoconus and ferroptosis were significantly associated.With further pursuit using WGCNA analysis, genes were classified, and 8 genes related to ferroptosis were screened, and GO analysis demonstrated that most of them (NOS2, AKR1C3, CAV1, SNCA, SLC7A11) were closely related to oxidative stress and iron ion transport.Then, hub genes were screened by SVM-RFE, and 9 genes (AQP5, ATF2, ARNTL, AKR1C3, FXN, MGST1, ARNT, LCN2, SLC11A2, SMPD1) related to oxidative stress and metal ion transport were analyzed by GO.Finally, the two gene sets were intersected and AKR1C3 came to the fore.
AKR1C3, a member of the aldosterone-reductase (AKR) family, is a hormone activity regulator with various cellular functions, including the regulation of prostaglandin, steroid hormone, and retinoid metabolism 41 .AKR1C3 has been reported to have oxidoreductase activity and AKR (NADPH) activity 42 , and is involved in the regulation of cellular redox homeostasis and immunity of tumor cells 43 .As an inhibitor of ferroptosis, the reduction of AKR1C3 protein expression can stimulate the occurrence of ferroptosis 44 .The cystine/glutamate antiporter SLC7A11 (commonly known as xCT) can import cystine for glutathione biosynthesis and antioxidant defense, and high expression of SLC7A11 can effectively suppress ferroptosis, which has been clearly verified in many cancers 45,46 .However, studies have shown that overexpression of AKR1C3 can promote the increase of SLC7A11 in HCC cells 44 .In our analysis, SLC7A11 expression was significantly reduced, which may have been caused by the reduction of AKR1C3, which in turn inhibited the xCT system and promoted ferroptosis.Genetics has demonstrated that maintaining GSH synthesis or promoting the activity of the system xCT protects cells from death triggered by oxidative stress conditions 47,48 .Similarly, the accumulation of reactive oxygen species (ROS) is an important mechanism that causes oxidative stress death of endothelial cells, and ROS can promote cell death through ferroptosis 49 ; this seems to provide strong support for ROS accumulation in keratoconus.Similarly, evidence suggests that ferroptosis plays an important role in immune cell-induced cell death 50 .For example, NK cells have reportedly been involved in ferroptosis, lipid peroxidation and increased expression of proteins related to oxidative damage in tumor cells 51 , and macrophages showed resistance to ferroptosis 52 .
Prior reports have indicated that miR-184 is the most abundant miRNA that is expressed in corneal and lens epithelial cells, which is associated with the horizontal adjustment of proteins in the cornea 53,54 .In the present study, AKR1C3 was found to be responsible for the synthesis of miR-184, and the expression of AKR1C3 was decreased in keratoconus.Therefore, the synthesis of miR-184, which is responsible for the regulation of protein level, was also impaired, likely an important reason affecting KC.We also predicted IncRNA (MLLT4-AS1 and RP11-830F9.6), associated with miR-184 synthesis, which is essential for the molecular function of miR-184.The ultimate aim was to predict the corresponding immune drugs, which could serve in interfering with keratoconus at the genetic level.Of the 10 drugs predicted to target AKR1C3, 4 (INDOMETHACIN, DAUNORUBICIN, DOXORUBICIN, DOCETAXEL) were selected, principally because they target AKR1C3.Unfortunately, all except Indomethacin are confined to the treatment of tumors.Nevertheless, drug development aimed at the treatment of genes is still in its relative infancy.
In conclusion, our study revealed that the pathogenesis of keratoconus may be due at least in part to ferroptosis on the ocular surface, and the core gene AKR1C3 was identified by bioinformatics analysis, showed a connection to the expression of miR-184 and other miRNAs, as well as the over expression of inflammatory factors and immune subtype cells.

Figure 1 .
Figure 1.The workflow chart in this study.

Figure 5 .
Figure 5. GO analysis of ferroptosis genes.(A) GO analysis of differentially expressed genes screened by WGCNA.(B) GO analysis of differentially expressed genes screened by SVM-RFE.(C) Venn diagram shows intersection genes.

Figure 6 .
Figure 6.The differential gene expression and Receiver operating characteristic analysis (ROC).(A) The mRNA expression of AKR1C3.(B) The mRNA expression of SLC7A11.(C) The AUC analysis of AKR1C3.Receiver operating characteristic, ROC; AUC, area under the ROC curve.