A Single-Cell Transcriptomic Atlas Elucidates a Microglial Gene Signature Linking Ferroptosis to Mitochondrial Dysfunction in Epilepsy

Introduction

Epilepsy, a chronic brain disease affecting over 70 million people worldwide, is characterized by spontaneous and recurrent seizures. The pathogenesis of epilepsy results from the synergistic interaction of multiple mechanisms, including disrupted excitatory-inhibitory balance, neuroinflammation, ion channel abnormalities, and mitochondrial dysfunction.^1,2 Despite significant advances in treatment, approximately 30% of patients develop drug-resistant epilepsy, for whom surgical interventions have limited clinical applicability.³ Consequently, identifying novel therapeutic targets through comprehensive biological analyses of epileptogenesis has become an urgent research priority.

Ferroptosis, a regulated form of cell death characterized by excessive intracellular iron accumulation, has been increasingly recognized as a critical contributor to the progression of various neurological disorders.^4,5 Notably, mitochondria undergo distinct morphological and functional transformations in ferroptosis. Morphologically, these changes involve decreased bilayer membrane density, reduced volume, and diminished or absent mitochondrial cristae. Functionally, these modifications result in depolarized mitochondrial membrane potential (ΔΨm) and disrupted energy metabolism.^6,7 Concurrently, mitochondrial dysfunction exacerbates oxidative stress, disrupts oxidative phosphorylation (OXPHOS), lowers neuronal excitability thresholds, and amplifies neuroinflammation, further promoting seizure generation.^8,9 However, the specific molecular mechanisms linking ferroptosis and mitochondrial dysfunction in epilepsy, particularly across different cell types, remain poorly understood.

Traditional studies on epilepsy are mostly based on the brain tissue level, and it is difficult to reveal the heterogeneity of individual cell types and their specific roles in epileptogenesis. Conventional anti-seizure medications predominantly target broad-spectrum ion channels or neurotransmitter systems rather than specific pathogenic genes. Gene therapy represents a transformative therapeutic paradigm for epilepsy due to its capacity for precise targeting of pathogenic mechanisms. Existing limited studies have demonstrated significant alterations in ion channel function and synaptic transmission (involving KCNQ2/3, SCN1A, and STXBP1) in epilepsy pathogenesis.^10,11 Thus, there is an urgent need to elucidate more molecular bases of epilepsy at single-cell resolution, pinpointing key cellular subtypes and signaling pathways to provide a theoretical foundation for developing new therapeutic approaches.

Single-cell transcriptomic studies explore cellular heterogeneity, revealing pathogenic subpopulations and targetable pathways in epilepsy. Machine learning, a transformative approach for analyzing biological data, has been applied to systematically identify critical genes, significantly improving the efficiency and reliability of diagnostic biomarker discovery.^12–14 We combined snRNA-seq with bulk RNA-seq to explore distinct subpopulations within epileptic microenvironment. By integrating differential expression profiling, machine learning-based gene screening, intercellular communication networks, and pathway enrichment, we systematically delineated the functional roles of ferroptosis and mitochondrial dysfunction-associated molecules, particularly within microglia, during epileptogenesis. Drug-target predictions and molecular docking approaches suggested innovative treatment options. Using a pilocarpine rat model, we observed evidence of ferroptosis and mitochondrial dysfunction in epileptic rats, with immunohistochemical validation consistent with the bioinformatics findings. This research provides critical theoretical insights and potential therapeutic targets relevant to molecular subtyping and precision therapy in epilepsy.

Methods

Data Acquisition

We downloaded the epilepsy-related datasets GSE190452 and GSE256068 from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The single-cell GSE190452 dataset includes 4 epilepsy samples and 4 control samples. Bulk transcriptome data were obtained from the GSE256068 dataset, comprising 135 epilepsy samples and 27 control samples. Gene sets analyzed in this study included 484 FRGs from the FerrDb database (http://www.zhounan.org/ferrdb/)¹⁵ and 1136 MRGs from the MitoCarta3.0 database (http://www.broadinstitute.org/mitocarta/).¹⁶

Animals

Twenty-seven six-week-old male Wistar rats (200 ± 20 g), provided by Shandong Jinan Pengyue Laboratory Animal Breeding Co., Ltd., were housed in a specific pathogen-free grade laboratory. The environmental conditions were carefully controlled, with temperature maintained at 22 ± 2 °C, humidity at 55 ± 5%, and a 12-h light/dark cycle enforced.

Single-Cell Data Quality Control and Processing

Single-cell data were first processed using the Seurat package.¹⁷ Cells with a high proportion of mitochondrial gene expression typically exhibit low RNA expression, a hallmark of cellular damage and an indicator of programmed cell death pathways.¹⁸ We employed the median absolute deviation (MAD) method to exclude cells with mitochondrial gene expression values exceeding 3 MADs above the median.¹⁹ Then, we used DoubletFinder (version 2.0.4)²⁰ to detect and remove doublets within each sample separately.

The filtered data were globally standardized using LogNormalize, and cell cycle scores were calculated using CellCycleScoring. Hypervariable genes were detected by FindVariableFeatures. Concurrently, the ScaleData function was employed to regress out technical variations due to mitochondrial gene content, ribosomal gene content, and cell cycle phase differences. Principal component analysis (PCA) was applied to reduce the linear dimensionality of the expression matrix. Harmony was used to remove batch effects, followed by nonlinear dimensionality reduction via uniform manifold approximation and projection (UMAP). Cell annotation was conducted by integrating marker genes from the CellMarker database (http://117.50.127.228/CellMarker/),²¹ the PanglaoDB database (https://panglaodb.se/),²² literature reviews, and automated annotation with SingleR software.²³

Machine Learning

Least Absolute Shrinkage and Selection Operator (LASSO) logistic regression was implemented using the glmnet package.²⁴ In this study, Random Forest (RF) averaged its outcomes to rank feature importance via %IncMSE, defined as the percentage increase in prediction error after feature permutation. The top 10 diagnostic gene candidates were retained for downstream analysis. By recursively removing the features contributing the least to classification performance, the Support Vector Machine–recursive feature elimination (SVM-RFE) algorithm selected the variables most critical for distinguishing samples.

Immune Infiltration

Single-sample gene set enrichment analysis (ssGSEA) was applied to the gene expression profiles to estimate the relative abundance of 29 distinct immune cell types.²⁵ What’s more, the correlations between the expression levels of target genes and the estimated proportions of various immune cells were systematically evaluated.

Pathway Dynamics and Transcriptional Regulatory Network Analyses

GSEA was conducted to identify differentially activated pathways between patient cohorts stratified according to key gene expression thresholds.²⁶ The background gene set was downloaded from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/) (version 7.0). This dataset served as the annotated background for conducting pathway subtype analyses. Differences in pathway expression between groups were assessed, and enriched gene sets (adjusted p < 0.05) were ranked by consistency score. By comprehensively scoring the gene sets of interest, gene set variation analysis (GSVA) transforms gene-level alterations into pathway-level changes, thereby enabling the inference of biological functions underlying the samples.²⁵ The R package RcisTarget was employed to predict transcription factors.²⁷ The normalized enrichment score (NES) for each motif is influenced by the total number of motifs present in the database. Besides the motifs originally annotated in the source dataset, we created supplementary annotation files based on motif similarity and gene sequence information.

Drug Target Prediction and Protein-Ligand Interaction Analyses

Drug prediction provides crucial support for target gene research and serves as an important link between basic research and clinical applications. We employed the Drug-Gene Interaction Database (DGIdb) to identify drugs or compounds that exhibit molecular interactions with key target genes. To further elucidate the interactions between candidate drugs and their target genes, molecular docking analyses were conducted to characterize interaction mechanisms. The structure of the protein encoded by each target gene was retrieved from the AlphaFold database (https://alphafold.com/), and the corresponding drug structures were obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). Molecular docking was performed using AutoDock software with nine independent runs. The docking result with the lowest binding energy, excluding the first reference conformation, was selected for presentation. Based on the molecular docking results, molecular dynamics (MD) simulations were performed using GROMACS 2023 software.²⁸

Lithium Chloride-Pilocarpine Model

After a one-week acclimatization period, adult rats were randomly assigned to two groups: (1) Control group, in which rats received an intraperitoneal (i.p.) injection of saline (10 mL/kg); (2) Epilepsy group, in which rats were pre-injected with lithium chloride (127 mg/kg, i.p., MilliporeSigma, Shanghai, China) 16–20 h prior. They were then administered scopolamine (1 mg/kg, i.p., MCE, Shanghai, China), followed 30 min later by pilocarpine hydrochloride (30 mg/kg, i.p., Targetmol, Shanghai, China). Seizure behaviors over the subsequent 120 min were recorded and graded according to the Racine scale as previously described:²⁹ I, chewing and facial movements; II, neck muscle spasms and head nodding; III, bilateral forelimb clonus (with or without upright posture); IV, tail hyperextension and rearing with loss of balance; V, generalized tonic-clonic seizures with limb rigidity. Rats exhibiting stage IV or V seizures were considered to have successfully developed status epilepticus. Three rats in the epilepsy group died and were excluded from the analysis. The rats (n = 12 per group) were deeply anesthetized via inhalation of 3% isoflurane in an induction chamber. Euthanasia was then performed by transcardial perfusion, with the depth of anesthesia confirmed by the absence of a pedal withdrawal reflex.

Histopathological Examination

Fourteen days after seizures, rat brains were extracted, immersed in 4% paraformaldehyde overnight, and embedded in paraffin. Coronal sections (4 μm) were prepared using a rotary microtome and mounted on slides. As previously described,³⁰ hematoxylin and eosin (H&E) staining was performed to assess hippocampal architecture. In addition, Prussian blue staining was employed to evaluate iron deposits in the hippocampus and cortex. Following the instructions of the respective kits, stained sections were visualized using a light microscope (Leica, Wetzlar, Germany).

Transmission Electron Microscopy (TEM)

Fourteen days after seizures, fresh hippocampal tissue (approximately 1 mm ³ in volume) was fixed in 2.5% glutaraldehyde for 24 h. After dehydration using graded acetone, the hippocampal samples were embedded in epoxy resin and sectioned into 60–80 nm ultra-thin slices. Following staining with uranyl acetate for 10 min and lead citrate for an additional 10 min, mitochondrial ultrastructure features were examined using TEM (Hitachi, Tokyo, Japan).

JC-1 Fluorescence Staining

Mitochondria were isolated from fresh rat hippocampal tissues using a tissue mitochondria isolation kit (Beyotime, Shanghai, China). Mitochondrial extract was rapidly stored in storage buffer at 4 °C and quantified for protein content. Detection of JC-1 polymers and monomers was conducted using the enhanced ΔΨm assay kit (Beyotime, Shanghai, China) and visualized under an inverted fluorescence microscope (Leica, Wetzlar, Germany).

IHC

Brain tissues were sectioned into 4 μm slices as described previously. Following dewaxing and rehydration, the sections underwent water-bath heating for antigen retrieval. Next, they were incubated overnight at 4 °C with an anti-PARP14 antibody (1:100, Affinity Biosciences, Jiangsu, China) and then treated with an anti-rabbit secondary antibody (1:100, Affinity Biosciences, Jiangsu, China). After staining with diaminobenzidine solution (ZSGB-BIO, Beijing, China) and counterstaining with hematoxylin (Biosharp, Anhui, China), all mounted sections were photographed using a microscope (Leica, Wetzlar, Germany).

Statistical Analysis

The R software package (version 4.2.2) was used to conduct all statistical analyses. Comparisons between groups were performed using the Student’s unpaired t-test. Statistical significance was defined as a p-value < 0.05.

Results

The snRNA-Seq Profiling of Epilepsy

To ensure data quality, the Seurat package was used to filter out captured outliers and cells expressing fewer than 200 genes. Subsequently, the DoubletFinder package was utilized to identify and exclude doublets, leading to a total of 41,218 retained cells. The filtering results were visualized by violin and scatter plots (Figure S1a and b). Next, 2000 hypervariable genes were identified (Figure S1c), followed by data normalization, homogenization, PCA, and Harmony integration (Figure S1d–f).

Finally, 13 cellular subgroups were identified via UMAP analysis (Figure 1a). Among these, 7 cell types were classified: oligodendrocyte precursor cell (OPC), oligodendrocyte, inhibitory neurons, excitatory neurons, microglia, astrocyte, and endothelial cell (Figure 1b). Additionally, a bubble plot illustrating the classical markers of each cell type (Figure 1c) and bar plots showing the cell proportions for each group (Figure 1d) were generated.

Figure 1 Single-cell annotation. (a) Uniform manifold approximation and projection (UMAP) visualization revealing 13 cellular subgroups. (b) Annotation of the 13 subgroups into 7 cell categories including oligodendrocyte precursor cell (OPC), oligodendrocyte, inhibitory neurons, excitatory neurons, microglia, astrocyte, and endothelial cell. (c) Bubble plot showing classical markers of the 7 cell types. (d) Bar charts depicting the proportional composition of each cell type in both groups.

FRGs Analysis in Microglia

FRGs sets in single-cell data were evaluated using five distinct methods: AUCell, UCell, singscore, ssGSEA, and AddModuleScore.³¹ These scores were subsequently integrated into a comprehensive composite metric. The results demonstrated that ferroptosis scores were higher in microglia, with significant upregulation observed in the disease group (Figure 2a and b). To compare the disease group with controls, we applied the FindMarkers function to identify 606 differentially expressed genes (DEGs) in microglia, applying a stringent threshold of logFC > 0.25 and adjusted p-value < 0.05 (Figure 2c). By intersecting these DEGs with 484 FRGs, we obtained 20 overlapping genes involved in ferroptosis (Figure 2d). These genes included CYBB, FTL, FTH1, TMSB4X, TGFBR1, TLR4, BACH1, RBMS1, MEF2C, PARP14, SAT1, HSPA5, ALOX5, PARP8, ZEB1, PTPN18, GJA1, HMOX1, ATM, and JUN.

Figure 2 Heterogeneity in ferroptosis-related genes (FRGs) expression. (a) Bubble plot showing FRGs expression scores across the 7 cell types. (b) Violin plot showing differences in FRGs scores (****P < 0.0001). (c) Volcano plot representing differentially expressed genes (DEGs) identified. (d) Venn diagram illustrating the overlap between DEGs and FRGs.

Correlation and Machine Learning Analysis of Mitochondrial Genes in Bulk RNA-Seq Identifies Epilepsy Biomarkers

A total of 1136 MRGs were obtained, and mitochondrial pathway activity in the samples was subsequently quantified using ssGSEA. Correlation analysis was then performed to assess the association between ssGSEA scores and the expression of 20 overlapping FRGs, with a threshold of p < 0.05 (Figure 3a–L). Finally, 12 genes significantly related to mitochondrial function were identified, including ALOX5, ATM, BACH1, CYBB, GJA1, HMOX1, JUN, PARP8, PARP14, TGFBR1, TLR4, and ZEB1.

Figure 3 Correlation analysis between intersecting genes and mitochondrial scores. (a–l) Scatter plots showing the correlations between mitochondrial pathway activity, quantified by single-sample gene set enrichment analysis (ssGSEA) scores, and the expression levels of ALOX5, ATM, BACH1, CYBB, GJA1, HMOX1, JUN, PARP8, PARP14, TGFBR1, TLR4, and ZEB1.

We applied three machine learning methods to identify key epilepsy-associated genes from the initial 12 genes. LASSO regression identified 7 genes as potential signature biomarkers for epilepsy (Figure 4a and b). Complementarily, the RF algorithm selected the top 10 most relevant candidates (Figure 4c). For feature screening via SVM-RFE, the highest classification accuracy was achieved with 9 feature genes (Figure 4d). The intersection of the three algorithms revealed six consensus genes (ATM, TGFBR1, ZEB1, PARP14, GJA1, and ALOX5) as the most promising signatures for further investigation (Figure 4e).

Figure 4 Construction of machine learning models. (a) The trajectories of Least Absolute Shrinkage and Selection Operator (LASSO) regression coefficients and the optimal gene combination at the minimum lambda value. (b) Minimum lambda determined via ten-fold cross-validation in LASSO regression. (c) Top 10 signature genes identified by the Random Forest (RF) algorithm. (d) Top 9 signature genes selected by the Support Vector Machine-recursive feature elimination (SVM-RFE) algorithm, achieving peak classification accuracy. (e) Venn diagram showing six consensus genes identified by all three algorithms.

Immune Infiltration Analysis

We characterized the immune infiltration landscape by quantifying immune cell proportions and their intercorrelations (Figure 5a and b). Compared with controls, disease samples exhibited significantly elevated APC co-stimulation and Th1 cell levels but reduced CD8+ T cell abundance (Figure 5c). Furthermore, the relationships between hub genes and immune cell abundances were explored. Specifically, ATM, TGFBR1, PARP14, GJA1, and ALOX5 showed significant negative correlations with CD8+ T cells while positively correlating with multiple immune features, such as APC co-inhibition, CCR, checkpoint, and HLA. ZEB1, TGFBR1, and GJA1 displayed strong negative correlations with T cell co-stimulation (Figure 5d). In the immune microenvironment, systematic evaluation of hub genes and various immune factors, including chemokines, receptors, major histocompatibility complexes (MHC), immunoinhibitors, and immunostimulators, revealed that these key genes significantly modulate immune cell infiltration (Figure 6a–e).

Figure 5 Immune infiltration landscape. (a) Compositional analysis of 23 types of immune-infiltrating cells. (b) Pearson’s correlation coefficients among 23 immune-infiltrating cells. (c) Comparative violin plots of differentially abundant immune cells. (d) Heatmap illustrating correlation coefficients between candidate genes and various immune cell populations. *P < 0.05, **P < 0.01, ***P < 0.001, ns = no significance.

Figure 6 Hub gene-immune factor correlations. (a–e) Significant associations between critical genes and chemokines (a), receptors (b), major histocompatibility complexes (MHC) (c), immunoinhibitors (d), and immunostimulators (e).*P < 0.05, **P < 0.01.

Functional Enrichment and Transcriptional Regulation of Epilepsy Hub Genes

GSEA showed that ATM is enriched in cytokine-cytokine receptor interaction, antigen processing and presentation, and Toll-like receptor signaling pathway (Figure 7a). TGFBR1 is enriched in NF-κB signaling pathway, B cell receptor signaling pathway, and Th17 cell differentiation (Figure 7b). ZEB1 is enriched in Th1 and Th2 cell differentiation, JAK-STAT signaling pathway, and NOD-like receptor signaling pathway (Figure 7c). PARP14 is enriched in NF-κB signaling pathway, TNF signaling pathway, and IL-17 signaling pathway (Figure 7d). GJA1 is enriched in ferroptosis, JAK-STAT signaling pathway, and PPAR signaling pathway (Figure 7e). ALOX5 is enriched in ferroptosis, Toll-like receptor signaling pathway, and NOD-like receptor signaling pathway (Figure 7f).

Figure 7 Gene set enrichment analysis (GSEA) of hub genes associated with epilepsy. (a–f) Significantly enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways identified for ATM (a), TGFBR1 (b), ZEB1 (c), PARP14 (d), GJA1 (e), and ALOX5 (f), respectively.

GSVA revealed that ATM is enriched in mitotic spindle and interferon-alpha response (Figure 8a). TGFBR1 is enriched in allograft rejection and androgen response (Figure 8b). ZEB1 is enriched in IL2-STAT5 signaling and KRAS signaling up (Figure 8c). PARP14 and ALOX5 are both enriched in interferon-gamma response and interferon-alpha response (Figure 8d and f). GJA1 is enriched in Notch signaling and TGFB signaling (Figure 8e).

Figure 8 Gene set variation analysis (GSVA) of hub gene pathway activity. (a–f) Pathways enriched in high-expression and low-expression clusters for ATM (a), TGFBR1 (b), ZEB1 (c), PARP14 (d), GJA1 (e), and ALOX5 (f), using the Hallmark gene set as the background.

Using the six key genes (ATM, TGFBR1, ZEB1, PARP14, GJA1, and ALOX5) as input, we performed motif enrichment analysis. The cisbp__M1271 motif emerged as the top-scoring regulatory element (NES: 5.72). Figure 9a and b displays all significantly enriched motifs alongside their associated transcription factors, revealing the transcriptional regulatory network during epileptogenesis.

Figure 9 Transcriptional regulatory network of epilepsy-associated genes. (a) Network illustrating the connections between six key genes (ATM, TGFBR1, ZEB1, PARP14, GJA1, and ALOX5) and their corresponding transcription factors. (b) Enriched binding motifs and their cognate transcription factors, with cisbp__M1271 (NES: 5.72) identified as the most significant regulatory element.

Single-Cell Expression Profiling of Key Genes and Immunometabolic Pathways

Using the Seurat R package, we analyzed single-cell expression patterns of key epilepsy-associated genes via the DotPlot and FeaturePlot functions, and found that PARP14, ATM, TGFBR1, and ALOX5 were highly expressed in microglia (Figure 10a and b). We employed AUCell to quantify the activity of genes related to immunometabolic pathways and used bubble plot for visualization to confirm expression differences of genes involved in these pathways. AUCell quantification demonstrated significant enrichment of these genes in immunometabolic pathways, particularly Notch signaling, mitotic spindle, and PI3K-Akt-mTOR signaling (Figure 10c).

Figure 10 Single-cell expression and immunometabolic pathway activity of key genes. (a) Featureplot showing the normalized expression levels of key genes across 13 distinct cell clusters. (b) Dotplot depicting cell-type-specific expression distributions, with the highest expression in microglial populations. (c) Bubble plot illustrating differential activity of immunometabolic pathways associated with key genes.

Drug Target Prediction Analysis

We utilized DGIdb to predict potential antiepileptic drugs targeting the key genes (ATM, TGFBR1, ZEB1, PARP14, GJA1, and ALOX5). This analysis may accelerate the development of antiepileptic drugs and lay a foundation for future translational validation. As shown in Table S1, 78 drugs were predicted for ATM, including TUVUSERTIB and ATR INHIBITOR ART0380. Seven drugs were predicted for TGFBR1, including VACTOSERTIB and METELIMUMAB. Three drugs were predicted for ZEB1, including SALINOMYCIN and OXALIPLATIN. Eleven drugs were predicted for GJA1, such as BLEOMYCIN and SELENIUM. Thirty-four drugs were predicted for ALOX5, including ZILEUTON and ATRELEUTON.

Molecular Docking and MD Simulation Reveal Stable Binding of GJA1 with Bleomycin

Based on interaction scores and drug approval status, GJA1 and bleomycin were selected as the exploratory drug target and ligand in this study. The docking binding energy is presented in Table 1. Structural docking results are visualized in Figure 11a. These results show that GJA1 exhibits strong binding affinity to bleomycin, as reflected by the binding energy (−6.0 kcal/mol).

Table 1 Binding Energy of Bleomycin to GJA1

Figure 11 Molecular docking and molecular dynamics (MD) simulation results of GJA1 with bleomycin. (a) Docking structure of the GJA1-bleomycin complex visualized using PyMOL software. (b) Root-mean-square deviation (RMSD) plot elucidating fluctuations of the complex. (c) Radius of gyration (Rg) plot evaluating the structural compactness of the complex. (d) Root-mean-square fluctuation (RMSF) plot depicting residue-specific flexibility and structural stability of the complex. (e) Hydrogen-bond plot identifying stabilizing intermolecular interactions within the complex. (f) Solvent-accessible surface area (SASA) plot reflecting intermolecular interactions and solvation effects of the complex.

The GJA1-bleomycin complex revealed significant fluctuations in the ligand binding site between 40,000 and 60,000 ps (Figure 11b and f), indicative of binding site rearrangement. Following this transition period (>60,000 ps), the complex achieved structural stability, as evidenced by three convergent metrics: (1) convergence of root-mean-square deviation (RMSD) values (Figure 11b); (2) stabilization of the radius of gyration (Rg) (Figure 11c); and (3) plateauing of the buried solvent-accessible surface area (SASA) (Figure 11f). These collectively demonstrate the establishment of a new stable binding pose. Residue-specific flexibility analysis via root-mean-square fluctuation (RMSF) measurements showed limited fluctuations (Figure 11d), with most residues exhibiting RMSF values below 2 Å. Stabilized hydrogen-bond occupancy indicates the persistence of intermolecular hydrogen bonds between the protein and ligand (Figure 11e). Collectively, these structural and dynamic analyses confirm that the GJA1-bleomycin complex forms a stable binding interface.

Morphological and Functional Characterization in Pilocarpine-Induced Epileptic Rats

H&E staining was performed for histological analysis to assess pilocarpine-induced hippocampal neuronal damage. Compared to the respective control group, the hippocampal dentate gyrus (DG), cornu ammonis 1 (CA1), and cornu ammonis 3 (CA3) regions of epileptic rats exhibited significant neuronal loss, condensed nuclei, swollen cytoplasm, widened intercellular gaps, and irregular cellular morphology (as indicated by arrows, Figure 12a). To elucidate the involvement of ferroptosis in epileptogenesis, Prussian blue staining was employed to detect free iron accumulation in the tissue. In the epilepsy group, distinct Prussian blue-positive hemosiderin deposits were observed, indicating pathological iron deposition (Figure 12b).

Figure 12 Histopathological examination and mitochondrial assessment. (a) Representative hematoxylin and eosin (H&E) staining images. Green arrows indicate neurons with normal morphology. Red arrows point to representative pathological features, including neuronal loss and nuclear condensation. (b) Prussian blue staining of the hippocampus and cerebral cortex. Red arrows show iron particle deposition in epilepsy rats. (c) Representative transmission electron microscopy (TEM) images (Arrowheads indicate mitochondria). (d) Representative images using the JC-1 fluorescent probe. (e and f) Quantitative analyses of fluorescence intensities, showing JC-1 aggregates (e) and JC-1 monomers (f) Data are presented as mean ± SD. n = 6. ****P < 0.0001.

Epileptic seizures cause significant structural and functional mitochondrial impairments in the hippocampus. Representative images of TEM revealed intact nuclear membranes and clearly defined mitochondrial cristae in controls. As shown in Figure 12c, epileptic rats exhibited pathological features, including increased mitochondrial membrane density, cristae fragmentation, irregular swelling, and vacuolar degeneration (indicated by arrows). Notably, autophagosomes—double-membrane structures encapsulating damaged mitochondria—were markedly more abundant in the epileptic hippocampus. To further evaluate mitochondrial dysfunction associated with ferroptosis, JC-1 staining was conducted to evaluate changes in ΔΨm. The epilepsy group showed increased levels of JC-1 monomers (green) and decreased levels of JC-1 aggregates (red) compared with the control group, indicating a collapse of ΔΨm during pilocarpine-induced ferroptosis (Figure 12d–f). This finding aligns with ultrastructural mitochondrial defects observed by TEM, collectively suggesting ferroptosis-driven bioenergetic failure.

To further elucidate the expression of candidate genes, IHC staining was performed to analyze different hippocampal regions in both groups. PARP14 expression was significantly higher in pilocarpine-induced epileptic rats compared to controls, validating our previous findings (Figure 13a–d).

Figure 13 IHC analysis of PARP14 expression. (a) Representative IHC staining images of PARP14. Quantitative analyses of PARP14 protein expression in the hippocampal (b) dentate gyrus (DG), (c) cornu ammonis 1 (CA1), and (d) cornu ammonis 3 (CA3) subregions. Data are presented as mean ± SD. n = 6. *P < 0.05, ***P < 0.001, ****P < 0.0001.

Discussion

Our study focuses on the roles of ferroptosis and mitochondrial dysfunction in epileptogenesis. We employed single-cell transcriptome sequencing combined with machine learning to identify pivotal candidate genes. These findings establish a molecular framework that could inform the future development of diagnostic tools and provide new directions for the discovery of antiepileptic drugs.

Through integrative analysis of multiple single-cell datasets, we identified seven distinct cellular subpopulations to investigate alterations in cell-type proportions following epilepsy. Leveraging FerrDb, we employed five computational methods to assess ferroptosis activity at single-cell resolution. Intriguingly, microglia exhibited markedly elevated ferroptosis scores, with a significant increase in epilepsy samples. As the primary immune cells in the brain, microglia accumulate iron deposits that interact synergistically with inflammatory responses, thereby exacerbating damage to the central nervous system.³² Microglia, among the most iron-retentive cells in the brain, exacerbate ferroptosis-driven neurodegeneration via the vesicle trafficking gene SEC24B.³³ As versatile sentinels of the central nervous system, microglia represent promising therapeutic targets; our investigation into their regulatory mechanisms in epilepsy thus provides novel insights for developing targeted anti-seizure interventions.

Building upon the identified microglial gene signatures at the intersection of epilepsy and ferroptosis, we integrated mitochondrial functional correlation analysis with three machine learning algorithms, ultimately identifying six hub genes (ATM, TGFBR1, ZEB1, PARP14, GJA1, and ALOX5). This study uncovers novel genes that link ferroptosis and mitochondrial dysfunction in epilepsy, which have not been previously described. ATM (ataxia telangiectasia mutated) acts as the primary sensor of DNA double-strand break damage. Activation of ATM, mediated by DNA damage, promotes acetylation of p53 and NF-κB, which exacerbates neuroinflammation and accelerates cellular senescence, thereby indirectly contributing to epileptogenesis.^34,35 ATM deficiency increases hippocampal neuronal excitability through dysregulated expression of the potassium-chloride cotransporter KCC2 and impaired gamma-aminobutyric acidergic (GABAergic) synaptic inhibition.³⁶ ATM promotes ferroptosis by phosphorylating the nuclear receptor coactivator 4, which increases intracellular labile iron levels.³⁷ Moreover, its deficiency disrupts mitochondrial function, manifesting as impaired mitophagy, aberrant mitochondrial accumulation, dysregulated energy metabolism, and elevated reactive oxygen species (ROS) levels.³⁸ TGFBR1 (transforming growth factor-beta receptor 1) orchestrates diverse pathophysiological processes, such as cell proliferation, differentiation, and apoptosis, by mediating signaling cascades of the TGFB superfamily.³⁹ In temporal lobe epilepsy (TLE) patients, upregulated TGFB1 expression in brain tissue activates the TGFBR1-Smad2/3 signaling pathway, thereby driving astrocyte proliferation and glial scar formation.⁴⁰ Additionally, microglia modulate myelin growth and integrity via the TGFB1-TGFBR1 axis in neurodegenerative diseases.⁴¹ ZEB1 (zinc finger E-box binding homeobox 1), a master regulator of epithelial-mesenchymal transition, modulates ferroptosis and mitochondrial function through heat shock protein family A (Hsp70) member 5 (HSPA5)-mediated pathways.⁴² Previous studies have identified ZEB1 as a pathogenic gene in Fuchs’ corneal dystrophy (FCD), a condition linked to epileptogenic lesions.⁴³ GJA1 (gap junction protein alpha 1) forms functional gap junctions between astrocytes, propagating pro-inflammatory signals that exacerbate epileptogenesis.⁴⁴ It regulates mitochondrial homeostasis in cardiomyocytes by coordinating ion transport, energy metabolism, and stress defense mechanisms. In vitro silencing of GJA1 markedly impairs mitochondrial energetics, as evidenced by depolarized ΔΨm and diminished adenosine triphosphate (ATP) synthesis.^45,46 ALOX5 (arachidonate 5-lipoxygenase), a key enzyme in arachidonic acid metabolism, drives ferroptosis by mediating lipid peroxidation. Mitochondrial DNA stress and dysfunction exacerbate ALOX5-dependent ferroptosis via activation of the cyclic GMP-AMP synthase (cGAS)/stimulator of interferon genes (STING) pathway.⁴⁷ Additionally, ALOX5 promotes glutamate hyperexcretion by activating the Ca2⁺/calmodulin-dependent protein kinase II (CAMKII)-mediated phosphorylation of synapsin I, thereby exacerbating seizure initiation.⁴⁸

This study provides the first evidence implicating PARP14 in epileptogenesis, offering a novel mechanistic perspective on epilepsy. As a member of the poly ADP-ribose polymerase (PARP) family, PARP14 modulates cellular stress responses and immune signaling pathways. Notably, activation of the PARP family is known to influence lipid metabolism.⁴⁹ NAD⁺ serves as an essential coenzyme for mitochondrial function and participates in the tricarboxylic acid (TCA) cycle. PARP1 and PARP2 impair mitochondrial homeostasis through NAD⁺ depletion.⁵⁰ PARP14, as a member of the PARP family, may similarly disrupt mitochondrial integrity via ADP-ribosylation activity-mediated NAD⁺ consumption. Emerging evidence demonstrates that PARP1 stabilizes ALOX5 via post-translational protein modification in neutrophils, thereby accelerating lung cancer progression.⁵¹ By analogy, PARP14 may drive epilepsy-associated neuroinflammation by promoting ALOX5-mediated lipid peroxidation. Previous studies have identified PARP14 as a shared ferroptosis marker gene in diabetic retinopathy and atherosclerosis, with its expression correlating with immune infiltration in patients with these conditions.⁵² Our study provides the first evidence that PARP14 contributes to epileptogenesis via ferroptosis and mitochondrial dysfunction, uncovering a previously unrecognized mechanistic link. We investigated the regulatory role of PARP14 in pilocarpine-induced rat models of epilepsy. IHC analysis revealed significant upregulation of PARP14 expression during epileptogenesis, consistent with our findings. A limitation of our study is its reliance on rat models, which may not fully recapitulate human epileptogenesis. In subsequent studies, we will selectively knockout PARP14 in microglia and employ adeno-associated virus (AAV)-mediated PARP14 overexpression in the hippocampus to further investigate its roles in epilepsy, with particular emphasis on its regulation of ferroptosis and mitochondrial dysfunction.

Studies have revealed that epileptic foci harbor a hyperinflammatory immune microenvironment, characterized by aberrant microglial activation, peripheral immune cell infiltration, and enhanced microglia–T cell crosstalk that exacerbates neuroinflammation.^53–55 A previous study demonstrated elevated Th1 cell levels in both TLE patients and animal models, where these cells modulate microglial activation—a finding consistent with our current immune infiltration analysis.⁵⁶ Biallelic loss-of-function mutations in TGFB1 cause central nervous system disorders characterized by epilepsy, cerebral atrophy, and posterior leukoencephalopathy.^39,57 Consistent with these prior findings, our study identified positive correlations between the expression of genes (ATM, TGFBR1, ZEB1, PARP14, GJA1, and ALOX5) and TGFB1 expression, reinforcing their potential synergistic roles in epileptogenic pathways. Their correlations with specific immune modulators also allow us to hypothesize how this gene signature might reshape concrete microglial effector functions, moving beyond general immune infiltration. First, the positive correlation with CSF1R (colony stimulating factor 1 receptor) suggests a potential modulation of the core signaling pathway that governs microglial general phagocytic capacity for cellular debris.⁵⁸ Second, the significant associations with the neuronal fractalkine (CX3CL1)-microglial fractalkine receptor (CX3CR1) signaling—displaying a decoupled pattern (negative correlation with the ligand CX3CL1 and positive correlation with the receptor CX3CR1)—point toward a specific dysregulation of the neuron-to-microglia communication vital for synapse-specific pruning.⁵⁹ Third, the strong correlations with a suite of MHC molecules, particularly MHC class II components, indicate an enhanced antigen-presenting and immunologically active phagocytic phenotype. ⁶⁰ Collectively, these associations suggest a microglial reprogramming that is poised to exacerbate neuroinflammation and synaptic instability, thereby reinforcing their synergistic roles in epileptogenesis.

Through systematic analysis of signaling pathways, we elucidated the biological mechanisms underlying epilepsy progression. GSEA revealed that ATM is enriched in cytokine–cytokine receptor interactions. Evidence has shown that seizures induce BBB disruption, which activates microglia and elevates pro-inflammatory cytokines, ultimately triggering inflammasome activation via Toll-like receptor signaling.^55,61 ALOX5 is enriched in pathways including Toll-like receptor signaling and ferroptosis. Toll-like receptor 4 (TLR4) is activated by the endogenous ligand high-mobility group box 1 (HMGB1), which promotes neuronal hyperexcitability through NF-κB pathway activation, thereby driving epileptogenesis.^62,63 TGFBR1 and PARP14 are involved in the NF-κB signaling pathway, a key driver of epilepsy. Activated NF-κB drives astrocyte proliferation and polarizes microglia toward the pro-inflammatory M1 phenotype, releasing inflammatory mediators and enzymes that exacerbate neuroinflammation and oxidative stress.⁶⁴ ZEB1 and GJA1 are enriched in the JAK-STAT signaling pathway, which regulates immune function and inflammatory responses. Clossen and his colleagues demonstrated that JAK-STAT inhibitors prolong seizure latency, reduce seizure frequency, and ameliorate synaptic plasticity impairments in epileptic models.⁶⁵ In addition, GSVA revealed that elevated expression of ATM, PARP14, and ALOX5 is linked to interferon-alpha signaling. Studies indicate that interferon-alpha signaling may contribute to epileptogenesis through multiple mechanisms, including neuroinflammation, synaptic plasticity dysregulation, and BBB disruption.^66,67 ZEB1 is highly expressed in the IL2-STAT5 signaling pathway. TGFBR1 and GJA1 both participate in the TGFB/Smad signaling pathway, coordinately regulating immunosuppressive processes. Integrative analysis of differential gene expression and functional enrichment reveals that these key genes contribute to epilepsy pathogenesis through interconnected pathways, including neuroinflammation, immune dysregulation, and ferroptosis.

Through DGIdb, we computationally predicted a panel of drugs potentially interacting with the key genes. A subsequent literature review identified several pharmacologic agents with potential therapeutic efficacy for epilepsy treatment. SULFASALAZINE inhibits the cystine/glutamate antiporter, thereby decreasing glutamate release and blocking glioma-associated seizures.⁶⁸ MELATONIN, an indoleamine primarily secreted by the pineal gland, has been demonstrated to attenuate epileptiform activity through its well-established neuroprotective, antioxidant, and anti-cytotoxic properties.⁶⁹ MONTELUKAST modulates cysteinyl leukotrienes—implicated in BBB disruption—effectively suppressing epileptogenesis and reducing seizure frequency in vivo.⁷⁰ SALINOMYCIN induces genotoxic DNA damage and exhibits synergistic effects with radiotherapy, significantly delaying glioblastoma progression.⁷¹ Glioblastoma and epilepsy are interconnected through complex bidirectional pathophysiological relationships involving tumor microenvironment alterations, aberrant neuronal excitability, and crosstalk between molecular pathways.⁷² FLUNARIZINE DIHYDROCHLORIDE, a selective calcium channel blocker, has been used for the prophylactic management of epilepsy, with particular efficacy in patients with comorbid migraine.⁷³ VALPROIC ACID, one of the most widely used antiseizure medications since its approval in 1978, exerts antiepileptic effects through multiple mechanisms, including GABAergic potentiation, ion channel modulation, and epigenetic regulation.⁷⁴ EVEROLIMUS, an mTOR inhibitor derived from rapamycin, reduces neuronal hyperexcitability by prolonging the opening of Ca²⁺ and K⁺ channels, thereby decreasing seizure frequency.⁷⁵ METHYLENE BLUE ANHYDROUS, a phenothiazine derivative, exhibits neuroprotective effects against acute seizures by scavenging ROS, modulating glutamatergic neurotransmission, and restoring mitochondrial bioenergetics.⁷⁶

Molecular docking and MD simulation are employed to analyze protein-ligand structures, investigating complex stability, dynamics, geometry, and binding affinity. These approaches provide temporal trajectories of atomic motions, offering novel hypotheses for structural interactions and drug development. This study represents the first application of temporal trajectory analysis to characterize the interaction between GJA1 and bleomycin. Our findings demonstrate that GJA1 forms a stable binding conformation with bleomycin, as evidenced by an RMSD < 2 Å, stable hydrogen-bond occupancy, and consistent structural stability and compactness. These structural data support a potential mechanism whereby bleomycin binding could modulate intercellular communication between neurons and astrocytes. Bleomycin, a glycopeptide antibiotic, exhibits dose-dependent pharmacological effects in vivo and is clinically employed as a first-line therapeutic agent for Hodgkin’s lymphoma and testicular carcinomas.⁷⁷ At low concentrations, bleomycin has been shown to activate the Kelch-like ECH-associated protein 1 (KEAP1)/nuclear factor erythroid 2-related factor 2 (Nrf2) signaling axis, resulting in the subsequent nuclear accumulation of Nrf2 and enhanced expression of Phase II detoxification enzymes such as glutamate-cysteine ligase and heme oxygenase-1.⁷⁸ At elevated concentrations, bleomycin forms redox-active complexes with intracellular ferrous iron, thereby inducing ferroptosis in neoplastic cells. Bleomycin-mediated activation of GJA1 in pulmonary epithelial cells promotes fibrogenesis.^77,79 Integrated GSEA and GSVA analyses reveal that GJA1 serves as a pivotal signaling hub within epileptic networks by orchestrating crosstalk among multiple pathways, including JAK-STAT, PPAR, ferroptosis, Notch, and TGFB signaling pathways. Intriguingly, bleomycin may modulate epileptogenesis through analogous mechanisms in neuronal systems. Future investigations should incorporate retrospective clinical analyses of seizure incidence in bleomycin-treated oncology patient cohorts to establish evidence-based dosing parameters that optimize therapeutic efficacy while minimizing neurotoxicity.

This study has certain limitations that warrant further improvement. The snRNA-seq samples included in this study are relatively limited, which may introduce bias and restrict the generalizability of our analysis. Additionally, while we focused on a microglia-enriched signature, the potential roles and interactions of the identified pathways in other critical cell types, including neurons and astrocytes, within the epileptic neurovascular unit were not fully explored. The functional implications of the key genes are inferred from bioinformatic correlations and pathway analysis. Although we observed consistent trends across public datasets, validation in diverse models (including clinical specimens) and across broader neuroanatomical regions will be essential to confirm the translational relevance of this gene signature. Furthermore, subsequent studies should conduct a series of functional verification experiments, including overexpression and knockout assays of the six key genes, to further confirm the reliability of our conclusions.

Conclusion

We identified a six-gene signature (ATM, TGFBR1, ZEB1, PARP14, GJA1, and ALOX5) centered in microglia that is associated with ferroptosis and mitochondrial dysfunction in epilepsy. Subsequent analyses explored their correlations with immune infiltration, signaling pathways, and transcriptional regulatory networks. In vitro experimental validation further demonstrated that PARP14 is upregulated during seizures, a previously unrecognized molecular insight for epilepsy research. Pharmacological prediction analysis based on these key genes provided a promising direction for mitigating seizures, with particular emphasis on bleomycin’s potential role in modulating GJA1. These findings prioritize these genes as high-value candidates for future functional and therapeutic exploration in epilepsy.