A Balancing Act in Corneal Epithelial Repair: A MAPK-JUN/EGR1/TFAP2A Network Regulates Ferroptotic Cell Fate

Introduction

The corneal epithelium is essential for maintaining corneal transparency and ocular surface integrity, both of which are critical for vision.^1–3 Precise epithelial regeneration following injury is indispensable for restoring visual function and preventing secondary complications such as infection, ulceration, and opacity.^4,5 However, due to continuous exposure to external environmental stressors, the corneal epithelium is highly susceptible to physical, chemical, and microbial damage.⁶ The efficiency of its repair directly determines corneal structural restoration and visual prognosis. Therefore, elucidating the molecular mechanisms that govern epithelial cell fate and tissue repair under stress conditions is fundamental for understanding corneal homeostasis and developing novel therapeutic strategies to promote epithelial regeneration.

Recent advances in high-throughput transcriptomics have enabled system-level characterization of molecular dynamics during tissue repair. Zebrafish, sharing a highly conserved corneal architecture and approximately 70–80% genomic homology with humans, has emerged as a powerful model for investigating ocular development and regeneration.^7,8 Notably, stress response and regeneration-associated pathways exhibit strong functional conservation between zebrafish and mammals, underscoring the model’s translational value in exploring corneal repair mechanisms.^9,10

Mounting evidence indicates that tissue repair is tightly coupled with signaling networks regulating cellular stress adaptation, survival, and migration.^11–13 Ferroptosis, a regulated form of cell death characterized by iron-dependent lipid peroxidation, has recently been implicated in various contexts of tissue injury and regeneration.^14,15 Impaired repair is often accompanied by the accumulation of ferroptotic markers, while inhibition of ferroptosis can promote regenerative recovery.¹⁶ Given that corneal injury typically involves oxidative stress, disrupted iron metabolism, and weakened antioxidant defenses, ferroptosis may represent a critical determinant of corneal epithelial integrity under stress.^17–20

Although studies have suggested that ferroptosis may cooperate with stress signaling pathways to regulate corneal epithelial repair, the specific mechanisms of the spatiotemporal dynamics and functional connection between the two remain unclear. Does ferroptosis participate in repair as a downstream effector module of a specific signaling pathway? Is its activation cell type or context specific? These key scientific issues have yet to be systematically analyzed.

To address these gaps, this study integrates transcriptomic analysis of a zebrafish corneal repair model with in vitro validation in human corneal epithelial cells. We aim to delineate the role of ferroptosis in epithelial regeneration and to uncover its molecular connection with key regulatory pathways. Our findings provide a new conceptual framework linking stress-responsive signaling to ferroptotic regulation, offering mechanistic insight and potential therapeutic targets for enhancing corneal epithelial repair.

Materials and Methods

Research Framework and Data Curation

The present study was conducted following the integrated analytical workflow outlined in Figure 1. Publicly available RNA sequencing data were obtained from the NCBI Gene Expression Omnibus under accession number GSE193784. The dataset comprised eight samples, including four uninjured corneal epithelial controls and four epithelial injury repair models. Raw expression data were comprehensively preprocessed, which included background correction using the RMA (Robust Multi-array Average) algorithm, quantile normalization, and log2 transformation prior to downstream analysis.

Figure 1 Overview of the analytical pipeline for transcriptomic data processing, hub gene identification, and experimental validation.

DEGs Screening and Key Module Identification

Screening of DEGs aims to accurately capture the transcriptome differences between corneal epithelium under physiological steady state and corneal epithelium during injury repair, thereby targeting gene sets with significant fluctuations in transcriptional profiles. This study performed differential gene expression analysis using the limma package (v3.52.4) in R, applying an empirical Bayes moderated t-test with Benjamini-Hochberg false discovery rate (FDR) correction. Genes with an absolute log₂ fold change |log₂FC| > 1 and an adjusted p-value < 0.05 were defined as significantly differentially expressed. The results were visualized using (1) volcano plots generated with ggplot2 (v3.4.0), and (2) hierarchically clustered heatmaps created with pheatmap (v1.0.12) using Euclidean distance and complete linkage.

To further capture gene groups that are not significantly differentially expressed but have cooperative fluctuation characteristics during the repair process, we performed WGCNA using the WGCNA package (v1.72–1) in R.²¹ We determined an optimal soft-thresholding, using the pickSoftThreshold function to achieve a scale-free topology. The adjacency matrix was transformed into a topological overlap matrix (TOM) to quantify network connectivity. Gene modules were identified using average linkage hierarchical clustering and a dynamic hybrid tree-cutting algorithm, with parameters set to a minimum module size of 30, deepSplit = 2, and a merge cut height of 0.25. The first principal component (module eigengene) of each module was extracted, and module-trait associations were assessed using Pearson correlation, controlling for potential confounders via linear regression. The DEGs obtained by differential expression analysis were intersected with genes in key related modules for subsequent research.

Functional Enrichment Analysis

The overlapping genes identified from the intersection of DEGs and the WGCNA module were subjected to GO and KEGG enrichment analyses to interpret their biological functions.

PPI Network Construction and Screening of Key Gene Sets

Subsequently, a PPI network was constructed using the STRING database (v11.5) with a confidence score threshold of > 0.4.²² The network was visualized and analyzed in Cytoscape (v3.9.1), and the cytoHubba plugin was used to compute multiple centrality measures (including Degree, Closeness, and Betweenness). Based on a comprehensive assessment of these topological features, the high-confidence hub genes were identified as central nodes in the network. Combined with indicators such as node connectivity, the key nodes in the network are sorted to obtain key gene sets.

Screening for Ferroptosis-Associated Hub Genes

To investigate the role of ferroptosis in corneal epithelial repair, ferroptosis-associated genes—including drivers, suppressors, and markers—were systematically curated from the FerrDb v2 database (http://zhounan.org/ferrdb/current/). A Venn diagram approach was used to intersect these ferroptosis-related genes with the previously identified network hub genes, leading to the identification of candidate genes that may functionally link ferroptosis to the corneal epithelial wound healing process.

Cell Culture and Ferroptosis Induction

The HCE-T line was procured from Procell Life Science & Technology Co., Ltd. (Wuhan, China) and originally sourced from the RIKEN BioResource Research Center (Catalog number: RCB2280). This cell line was derived from the corneal tissue of a 49-year-old female donor and immortalized by transformation with an SV40-adenovirus 12 (SV40-Ad12) hybrid vector. Cells were maintained in DMEM/F12 medium (Hyclone, Logan, USA), supplemented with 5% fetal bovine serum (Gibco), 5 mg/mL recombinant insulin (Beyotime, Shanghai, China), 10 ng/mL human epidermal growth factor (Thermo Fisher Scientific, Waltham, USA), 0.5% dimethyl sulfoxide (DMSO), and 1% penicillin–streptomycin (Beyotime) under standard conditions (37°C, 5% CO2, humidified atmosphere). When cells reached 70–80% confluence, cells were washed with PBS and treated with the ferroptosis-specific inducer 2.5 μmol/L Erastin for 24 hours to establish an in vitro ferroptosis model for subsequent experimental investigations.

Validation of Hub Gene Expression via RT-qPCR

RNA extraction was performed on treated cell samples utilizing a standardized commercial RNA isolation (GOONIE, China, Cat#400-100) according to the manufacturer’s protocol. The purified RNA was subsequently reverse-transcribed into complementary DNA (cDNA) using a cDNA synthesis kit (TIANYA BIO, China, Cat#P1504). Quantitative amplification reactions were carried out on a LightCycler 480 real-time PCR platform (Roche Diagnostics, Switzerland) with SYBR Green chemistry for fluorescence detection. For data normalization, the housekeeping gene GAPDH served as an endogenous reference, and relative quantification of target gene expression was determined through the comparative threshold cycle (2⁻ΔΔCt) analytical approach. The complete set of oligonucleotide primers designed for specific amplification is systematically presented in Table 1.

Table 1 Oligonucleotide Primer Sequences for Quantitative PCR Validation

Statistical Analysis

Statistical analyses were performed using R software and GraphPad Prism (v9.0). Data distributions were first assessed for normality (Shapiro–Wilk test) and homogeneity of variance prior to hypothesis testing. When parametric assumptions were satisfied, differences between groups were analyzed using two-way ANOVA followed by appropriate multiple-comparisons tests. If normality or variance assumptions were not met, corresponding nonparametric tests were applied. A two-sided p value < 0.05 was considered statistically significant. Quantitative data are presented as mean ± standard deviation (SD), with at least three independent biological replicates per group.

Results

Identification of DEGs and Key Co-Expression Modules

Following data normalization and quality control, a total of 14,914 genes were included in the analysis. Using the predefined screening thresholds of |log2 fold change| > 1 and adjusted p-value < 0.05, 252 genes were identified as significantly differentially expressed in the corneal injury repair model compared to the uninjured control. Among these, 170 genes were up-regulated and 82 were down-regulated. The overall distribution of these expression changes is visualized in the volcano plot (Figure 2A). Furthermore, to specifically highlight the most pronounced expression changes, a clustered heatmap was employed to visualize the expression patterns of the top 100 most significantly altered genes (comprising 50 up- and 50 down-regulated genes) across all samples, with all data being log₂-transformed and standardized (Figure 2B).

Figure 2 Continued.

Figure 2 Identification of DEGs and Key Co-Expression Modules. (A) Volcano plot of differentially expressed genes (DEGs). The plot visualizes the global distribution of gene expression changes. The vertical and horizontal blue dashed lines denote the filtering thresholds for |log2 fold change| > 1 and adjusted p-value < 0.05, respectively, which define the boundaries for statistical significance. Genes represented by black dots indicate those with non-significant differential expression that failed to meet these criteria. Statistically significant upregulated genes are highlighted in red (n = 170), while downregulated genes are shown in green (n = 82). A total of 252 DEGs were identified for further analysis. (B) Heatmap of the top 100 most variably expressed genes. The expression patterns of the top 50 upregulated and 50 downregulated genes across all samples are shown. Expression values were log2-transformed and normalized, with red and blue indicating relatively high and low expression levels, respectively. The clustering clearly segregates the injured and uninjured control samples. (C) Sample clustering dendrogram. Hierarchical clustering of all samples based on gene expression profiles shows no outliers, confirming data integrity for network construction. (D) Network topology analysis for various soft-thresholding powers. The Scale Independence plot (left) and Mean connectivity plot (right) show the fit index and connection density against candidate powers (β). In the Scale independence panel, the red horizontal line represents the fit threshold (R2 = 0.72) for model selection. The numeric labels (1–20) denote specific βvalues for each data point, indicating network stabilization at higher powers. An optimal power of β= 16 was selected where the Scale independence curve flattens above the threshold while maintaining biologically relevant Mean connectivity. (E) Cluster dendrogram of co-expressed genes. Genes were clustered based on topological overlap, with color bands below indicating module assignment. A total of 42 distinct co-expression modules were identified. (F) Module-trait relationships heatmap. This heatmap illustrates the correlation between module eigengenes (rows) and the corneal epithelial repair phenotype (columns). Each cell contains the Pearson correlation coefficient (r) and its associated P-value in parentheses. The color gradient from blue to red represents negative to positive correlations, respectively. Notably, the antiquewhite1 module demonstrates the strongest positive correlation with the repair phenotype (r = 0.94, P = 6e−04). (G) Scatterplot of Gene Significance (GS) versus Module Membership (MM) in the antiquewhite1 module. The scatterplot depicts the highly significant positive correlation was observed (Pearson correlation r = 0.84, P < 1e−200), indicating that genes most central to this module are also the most highly associated with the repair phenotype. (H) Venn diagram of gene intersection. The diagram shows the overlap between differentially expressed genes (DEGs, 252 genes) and genes from the WGCNA-derived antiquewhite1 module (1,630 genes), yielding 99 overlapping genes for subsequent analysis. The statistical significance of this overlap was evaluated using a hypergeometric test (P < 0.001).

Subsequently, a WGCNA was performed on the 14,914 genes to identify modules associated with corneal epithelial repair. Sample clustering analysis revealed no outliers, indicating good data integrity for subsequent analysis (Figure 2C). The optimal soft-thresholding power was determined as β = 16 to achieve a scale-free network topology while maintaining adequate connectivity (Figure 2D). This parameter setting resulted in the identification of 42 distinct co-expression modules through hierarchical clustering (Figure 2E). Analysis of module-trait relationships indicated that the antiquewhite1 module exhibited the strongest positive correlation with corneal epithelial repair (Figure 2F and G). Consequently, the 1630 genes comprising this module were defined as being associated with corneal epithelial repair. Finally, by intersecting the DEGs with the corneal epithelial repair-related genes, a set of 99 overlapping genes was obtained for further investigation (Figure 2H).

Functional Enrichment Analysis of Overlapping Genes

To decipher their biological roles, this gene set was performed on the 99 overlapping genes obtained above. GO analysis revealed significant enrichments in key categories, including regulation of the MAPK cascade (Biological Process, BP), transcriptional regulatory complex (Cellular Component, CC), and ion channel regulator activity (Molecular Function, MF) (Figure 3A). Concurrently, KEGG pathway analysis confirmed the MAPK signaling pathway as the most significantly enriched pathway (Figure 3B).

Figure 3 Functional enrichment analysis of overlapping genes. (A) GO enrichment analysis. The bar plot displays the significantly enriched GO terms in three categories: BP, CC, and MF. The most significantly enriched terms include regulation of the MAPK cascade (BP), transcriptional regulatory complex (CC), and ion channel regulator activity (MF). (B) KEGG pathway enrichment analysis. The bubble plot shows the significantly enriched KEGG pathways. The MAPK signaling pathway was identified as the most significantly enriched pathway. Bubble size represents the number of genes in each pathway, and color indicates the significance level (-log10 adjusted p-value).

PPI Network Analysis

Subsequently, a PPI network was constructed from these 99 genes (Figure 4A). Topological analysis of this network using multiple centrality measures via the cytoHubba plugin identified 15 high-confidence hub genes that occupied the most critical positions within the network architecture (Figure 4B).

Figure 4 PPI Network Analysis. (A) The PPI network was constructed from the 99 overlapping genes using the STRING database with a confidence score threshold > 0.4. Nodes represent proteins and edges represent functional associations. (B) Identification of hub genes in the PPI network. Topological analysis using the cytoHubba plugin in Cytoscape identified 15 high-confidence hub genes based on multiple centrality measures (degree, closeness, and betweenness). Node color and size represent the hub gene score, with darker colors and larger sizes indicating higher centrality in the network.

Screening and Functional Annotation of Ferroptosis-Associated Hub Genes

We performed an intersection analysis between the 15 identified hub genes and all ferroptosis-related genes (including drivers, suppressors, and markers) cataloged in the FerrDb v2 database (Figure 5). This analysis identified three genes that possess both network hub properties and a close association with ferroptosis: EGR1 was classified as a ferroptosis driver, while JUN and TFAP2A were annotated as ferroptosis suppressors.

Figure 5 Venn diagram illustrating the overlap between the 15 identified hub genes and ferroptosis-related genes obtained in FerrDb v2. The statistical significance of the overlap was evaluated using a hypergeometric test (P <0.05).

Experimental Validation of Hub Gene Expression Patterns

To experimentally verify the involvement of these identified genes in ferroptosis, we examined their expression patterns in an Erastin-induced ferroptosis model using HCE-T cells. RT-qPCR analysis confirmed significant dysregulation of all three hub genes (n = 3). Consistent with their functional annotations in FerrDb v2, EGR1 (a driver) demonstrated significant upregulation (p < 0.001), while both JUN and TFAP2A (suppressors) showed marked downregulation (p < 0.001 and p < 0.05, respectively) compared to controls (Figure 6).

Figure 6 RT-qPCR analysis of JUN, EGR1, and TFAP2A expression in Erastin-induced ferroptotic HCE-T cells. EGR1 was significantly upregulated, while JUN and TFAP2A were downregulated compared to control cells (n = 3). Data represent mean ± SD. Significance thresholds: *P < 0.05, ***P < 0.001.

Cross-Species Comparison of Transcriptomic Features

To assess whether the regulatory network identified in this study is conserved across species, we systematically compared the key signaling pathways and hub genes involved in the zebrafish model with recently published mammalian single-cell and multi-omics datasets. The comparative analysis revealed a high degree of functional consistency across species. First, our findings indicated a significant enrichment of MAPK signaling cascades and apoptosis-related pathways during corneal epithelial repair in zebrafish, which closely paralleled the transcriptional activation patterns observed in a cynomolgus monkey corneal wound healing model at day 1 post-injury.²³ Furthermore, the core ferroptosis-related hub genes and highly proliferative gene clusters identified in our network—particularly the MAPK-centered regulatory nodes—showed substantial overlap with the molecular signatures driving cellular plasticity and lineage reprogramming in murine models.²⁴ These features were also consistent with the hyperproliferative transcriptional profiles of CKS2⁺/STMN1⁺ transit-amplifying cells defined in the human corneal single-cell atlas.²⁵ Collectively, these cross-species comparative results support that iron-dependent metabolic stress responses, along with MAPK-driven stress repair regulatory programs, constitute core components of corneal epithelial injury and repair mechanisms.

Discussion

This study integrates transcriptomic data from a zebrafish corneal epithelial repair model with ferroptosis-related gene databases and in vitro validation, revealing a pivotal role for ferroptosis in corneal wound healing. We identified a regulatory network centered on the MAPK pathway, composed of transcription factors JUN, EGR1, and TFAP2A, which fine-tunes the balance between cell death and survival to drive effective epithelial regeneration.

Rather than viewing ferroptosis merely as a pathological consequence of lipid peroxidation, our findings position this iron-dependent cell death as an intrinsic, highly regulated component of the regenerative program. Corneal repair relies on the timely clearance of damaged cells and the efficient migration of healthy epithelium.^26–28 Our results support the concept that ferroptosis functions as a tightly regulated “molecular debridement” process facilitating tissue renewal.

Based on our bioinformatic analyses and in vitro Erastin-induced ferroptosis models, we propose a biphasic temporal model for ferroptosis during corneal repair: a “repair initiation phase” and a “repair execution phase”. Rapid upregulation of the immediate response gene EGR1 in the early stages of injury may enhance cellular sensitivity to ferroptosis by inhibiting key protective factors such as GPX4,^29–31 while downregulation of TFAP2A relieves the basal inhibition of ferroptosis, and the two synergistically trigger a local, controllable ferroptosis response.^32,33 This process is similar to “molecular debridement” by clearing damaged cells to create space for healthy cells to migrate. At the same time, the released DAMPs signal promotes inflammation and proliferation responses, thereby initiating the repair process.^34,35 As repair progresses, excessive ferroptosis becomes detrimental; JUN upregulation establishes an anti-ferroptotic barrier by activating antioxidant genes (eg., SLC7A11),³⁶ thereby preserving cell viability at the migratory front. The sustained expression of JUN and attenuation of pro-ferroptotic signals mark the transition from “clearance” to “reconstruction”.^37,38

The high degree of consistency observed at the level of key genes and signaling pathways between the zebrafish model and higher mammalian systems substantially enhances the translational value of our findings. According to our results, corneal repair involves profound transcriptomic reprogramming within relevant cell populations, a phenomenon closely aligned with observations by Feret et al,²⁴ in murine injury models, wherein the regenerative state accompanied by elevated metabolic demands is defined as “cellular plasticity” Functionally, this regenerative plasticity in the human cornea is primarily manifested within the highly proliferative CKS2⁺/STMN1⁺ transit-amplifying cells (TACs), which constitute the principal proliferative source for epithelial renewal and wound coverage.²⁵ Furthermore, studies by Zhou et al²³ have demonstrated that the acute phase of corneal epithelial regeneration is characterized by enhanced cellular plasticity concomitant with increased metabolic flux, intrinsically elevating cellular susceptibility to ferroptosis-related damage. Therefore, the MAPK–JUN/EGR1/TFAP2A regulatory axis should not be interpreted solely as a classical pro-proliferative signaling pathway, but rather defined as a dual-functional “stress–regeneration integrative hub.” While driving cell cycle activation and proliferative expansion, this axis concurrently functions as an evolutionarily conserved survival regulatory module that buffers metabolic stress by inhibiting lipid peroxidation and maintaining redox homeostasis, thereby sustaining regenerative potential and ensuring efficient wound closure.

Considering the corneal epithelium’s unique immune privilege, neural regulation, and tear film microenvironment likely shape the activity of the JUN/EGR1/TFAP2A network. Future studies will require the use of knockout models and spatial transcriptomics to validate causal relationships and dissect how local microenvironment govern ferroptotic outcomes.

While integrating previously isolated molecules into a coherent MAPK-centered ferroptosis network, we acknowledge certain limitations. Notably, the identification of these regulatory nodes may be subject to selection and appraisal bias, as the initial screening framework and database-driven appraisal were conducted within a focused scope. To mitigate this, future studies employing broader multi-omic integration and unbiased high-throughput functional screens are necessary to further validate these causal relationships.

Despite these constraints and the lack of direct in vivo mammalian validation, our proposed “precise clearance for regeneration” model redefines ferroptosis as a dual, context-dependent process in corneal repair. From a translational perspective, this suggests that therapeutic interventions should shift from global inhibition toward temporal modulation: early, transient activation may facilitate “molecular debridement,” while late-stage inhibition supports epithelial regeneration. This framework provides a promising, time-resolved strategy for targeting specific downstream protective axes.

Conclusion

This study elucidates a MAPK-driven transcriptional network centered on JUN, EGR1, and TFAP2A, positioning ferroptosis as a key “cell fate arbiter” in corneal epithelial repair. These findings advance our understanding of tissue regeneration and lay a conceptual and mechanistic foundation for time-resolved therapeutic modulation of ferroptosis in corneal healing.