WGCNA combined with GSVA to explore biomarkers of refractory neocortical epilepsy.

About two-thirds of epilepsy patients relapse within five years after surgery. It is significant to note that the limitations of current treatments stem from a lack of understanding of molecular mechanisms. In this study, Weighted Gene Co-expression Network Analysis (WGCNA) and Gene set variation analysis (GSVA) methods were used to analyze the total RNA data from 20 surgical removal samples (epileptogenic zone and irritative zone, EZ and IZ) of 10 Chinese patients with refractory neocortical epilepsy downloaded from the original microarray dataset (GSE31718) of the National Center for Biological Information -Gene Expression Omnibus database (NCBI-GEO). The late stages of the estrogen response pathway, the IL6-JAK-STAT3-signal pathway and G2 checkpoints are correlated with the EZ, whereas the early stages of the estrogen response pathway and TGF-β signal are more strongly expressed in the IZ. The allogeneic rejection, apical surface and the TGF-β signal are relevant to the high seizure frequency, the unfolded protein response and MYC-target are mostly expressed in patients with low-frequency seizures. Genes with high gene significance(GS) values that were correlated with seizure frequency include OSR2, CABP4, CAPSL, CYP4F8, and FRK in the pink module, and SH3GLB2, CHAC1 and DDX23 in the yellow module. The occurrence of EZ and IZ act on different biological mechanisms. The upregulated genes associated with seizure frequency include OSR2, CABP4, CAPSL, CYP4F8, and FRK, and the downregulated genes include SH3GLB2, CHAC1 and DDX23. The evidence of key genes and differential pathways obtained by WGCNA and GSVA may be biomarkers for novel preventive and pharmacological interventions in clinical practice.

Introduction

Epilepsy, one of the most common disorders in neurology, affects more than 70 million people in the world (Thijs RD et al., 2019). Prolonged clustered seizures increase the risk of social-psychological dysfunction, cognitive decline, as well as sudden death without expectation. The investigation indicated the prevalence of active epilepsy was about 6.38/1000, and the lifetime prevailing was about 7.60/1000 (Huang Q et al., 2020).

Neocortical epilepsy, is a specific type of epilepsy, mainly due to parts of the cerebral cortex not migrating to the correct location during the unborn and post-birth growth of the patient, causing other material to grow in the originally correct location. Neocortical epilepsy has a greater possibility of developing refractory epilepsy and a higher frequency of abnormal neuronal discharges than normal epilepsy disorders. Surgical resection as a treatment for refractory neocortical epilepsy. However, patients treated with surgery are still at risk for recurrence. In a retrospective comparison of patients undergoing surgery between 1999 and 2015, two-thirds relapsed within five years after surgery. MRI scans after recurrence showed progression of original lesions in 13% of patients and new lesions in 13% of patients, compared with early postoperative MRI (Petrik S et al., 2021) The situation indicated above may be caused by the lack of understanding of the molecular mechanisms of epilepsy.

As we all know, recognizing the location of the epileptogenic zone (zones of seizure initiation, EZ) and irritative zone (interictal epileptic spikes, IZ) has always been a key step in preoperative evaluation. And the relationship between the EZ and IZ has aroused debate. The researchers recorded the distribution of peaks between interphase seizures and topography of EZ using stereoelectroencephalography(SEEG)in patients with focal neocortical epilepsies. There is a 44% inconsistency between EZ and IZ (Bartolomei F et al., 2016). However, the mechanism of the difference between EZ and IZ has rarely been reported.

Gene set variation analysis (GSVA) evaluates pathway expression between different clinical characteristics. The study enumerated the differential pathways between the irritative and epileptic zones, and between high and low-frequency seizures.

Weighted Gene Co-expression Network Analysis (WGCNA) makes use of the concept that a disease should be seen as a network of interferences rather than the change of individual genes. As Gaiteri et al. have pointed out, a gene may alter clinical characteristics that differential expression analysis cannot detect (Gaiteri, C. et al., 2014). The gene net is separately in distinct modules according to expression resemblance, the modules can be associated with clinical patterns, and key gene identification. In this study, WGCNA was used to identify key genes associated with refractory neocortical epilepsy.

Gene Ontology(GO) and Kyoto Encyclopedia of Genes and Genomes(KEGG) analysis are gene annotation tools used in bioinformatics, including structural annotation and functional annotation. GO describes the genes in terms of molecular function (MF), biological process (BP) and cellular components(CC), KEGG lists the signaling pathways in which genes are located, which reveal potential mechanisms at the molecular level of the disease.

Method

Patients and specimens

The search term " neocortical epilepsy " was used to access the Gene Expression Omnibus (GEO) database at the National Center for Biological Information (NCBI, https:// www.ncbi.nlm.nih.gov/geo/) to search for published epilepsy gene dataset. The original microarray dataset (GSE31718) was obtained from Molecular Biology and Human Genetics, Tzu Chi University, Taiwan, China. The dataset includes the total RNA data from 20 surgical removal samples (the EZ and IZ identified by clinical evaluation by EEG, SPEC, MRI, and v-EEG) of 10 patients with refractory neocortical epilepsy in varying clinical characteristics. (Table 1) (Hsin Y et al., 2015).

Table 1

Distribution of clinical characteristics of 10 patients with neocortical epilepsy.

Phenotypic characteristics		Number
Gender	Male	18
	Female	2
Resection zone	Epileptogenic zone	10
	Irritative zone	10
Age-onset	Early(≤4 yrs)	4
	Late(>4 yrs)	14
Frequency	Low(≤10 times a month)	8
	High(>10 times a month)	10
Brain hemisphere	Left	8
	Right	12
Brain region	temporal lobe	8
	frontal lobe	8
	frontal lobe and temporal lobe	2
	temporal lobe and parietal lobe	2

Note.The rightmost column shows the number of samples, epileptogenic zone and irritative zone of 10 patients, for a total of 20 samples

Weighted gene co-expression network analysis (WGCNA)

The study of the epileptogenic and irritative zones surgically obtained from refractory neocortical epilepsy patients described the illness as a disturbing net of the gene-gene network. The network was constructed using the WGCNA package (version 1.71; http://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/.

) in R environment (version R*64 4.1.1). Firstly, the samples clustering method was applied to identify outliers and match them with features. No outliers as Fig. 1 A showed, and all 20 samples could be analyzed. And then, the 6873 genes (top 75% absolute median deviation) from the total RNA samples in the database were subjected to network construction. The Pearson correlation coefficient was used for obtaining gene co-expression similarity measures, topological overlap matrix (TOM) was used to filter weak connections during network construction. The soft threshold was set as 16, which is more in line with the scale-free distribution and provides an appropriate average connection for the subsequent establishment of co-expression modules. (Fig. 1B) Finally, a dynamic tree-cutting algorithm was used to select branches of the tree graph and realized the module identification.

Fig. 1

（A）Sample clustering: screening for the presence of outlier samples. （B）The influence of network topology on soft threshold power is analyzed.

The left panel shows a scale-free fitting index (Y-axis) as a function of soft threshold power (X-axis). The panel on the right shows average connectivity (degree, y-axis) as a functional soft-threshold power (x-axis). The power was set as 16 for further analysis. （C）Gene clustering analysis of module identification: Twelve co-expression modules were constructed according to topological overlap and assigning module colors. The gene dendrograms and relationships between the gene modules are at the top and bottom of the images. （D）Heatmap plot of the adjacencies in the eigengene network: positive or negative module correlations are shown in red or blue, respectively. (E) TOM between genes is depicted in a heat map plot, and the depth of red represents the correlation intensity between pairs of genes on a linear scale.

Besides, module-trait correlations were constructed by WGCNA considering demographic, clinical, and molecular phenotypic. The gene significance (GS), is the value with the definition as the association between specific traits and gene expression within modules. Module Membership (MM), a measurement of connections within a module (Langfelder P and Horvath S.2008), represents the correlation between a single module and its genes. Both were used for the selection of key genetic markers. And the correlations mentioned above were investigated using Pearson’s correlation. The modules presenting high scores(r > |0. 50| and P < 0.05) with the single treat could be selective for progressive study to identify the high GS value genes. All online crossover analysis was performed and the WGCNA algorithm could be screened with R *(Langfelder P and Horvath S.2008).

Gene set variation analysis (GSVA)

We employed the R package “GSVA 1.36.2” to score gene sets of biological processes for each sample. GSVA evaluates potential variation in pathway activity in each sample and compares pathway expression between the two groups of patients by pre-inputting a selected gene set. (Hänzelmann S et al., 2013). The import of the GSVA algorithm was a gene expressive matrix of log2 microarray expressive scores as well as a pre-defining gene set or pre-existing gene set database (MSig). The GSVA score was non-parametric, using K-like random walk statistics and negative values of a specific specimen as well as a gene set. [Ning S et al., 2020].

By applying the above methods, here are the problems to be investigated: (i) if in patients different clinical patterns are connected with transcriptional modules; (ii) if patients with different resection zone and seizure frequencies express varying pathways; (iii) if the key genes with high GS value in the trait-relevant module, were associated with epilepsy-correlated molecular mechanisms.

GO and KEGG

GO analysis was processed by accessing the Panther Classification System.

(http://pantherdb.org/) The key genes were input in the Gene list box to submit, with selecting Homo sapiens for organisation type and functional classification for analysis mode. Eventually, the information of key genes were enumerated in terms of MF, BP and CC.

The above key genes were respectively used as search terms to access the KEGG pathway at the GenomeNet Database Resources (www.genome.jp) to search for signaling pathways of the key genes.

Result

Co-expression module

Total RNA samples were obtained from 20 surgical removal of 10 patients with refractory neocortical epilepsy, and 6873 genes (top 75% absolute median deviation) from the total RNA samples were subjected to network construction by WGCNA. The final 12 transcriptional modules represented by different colors were identified (Fig. 1 C). The resulting eigengene network proved two meta-modules, as displayed in the graded clustering dendrography, here named I and II. (Fig. 3E) The adjacencies of modules in the eigengene network were shown in Fig. 1D. The correlation intensity between pairs of genes were depicted in Fig. 1E.

Fig. 3

(A) (B) Module-trait associations, red and blue, show positive as well as negative correlations, separately, with darker colors in the section and larger values outside the brackets indicating higher correlations. (C) (D) Scatterplot of individual module-single trait associations, the whole gene scores were plotted in MM (x-axis) vs. GS graphic (y-axis), the horizontal axis (MM) stands for mutual relation among gene and the module, as well as vertical axis (GS) represents the absolute value of the correlation between the gene and the phenotypic feature. (E) The traits listed that are significantly related to the module, (P < 0.05) marked with red and green stand for positive and negative correlations, separately; the bold underlined letters stand for high correlations (r > |0.50|). ERL presents estrogen-response-late, ERE presents estrogen-response-early. AR, AS, UPR, and MYC are abbreviations of allograft rejection, apical surface, unfolded protein response, and MYC-targets-V1, respectively.

GSVA

Here is pathway activity variation related to the resection zone(the EZ and IZ) in the Fig. 2A. The late stages of the estrogen response pathway, the IL6-JAK-STAT3-signal pathway and G2 checkpoints are correlated with the EZ, whereas the early stages of the estrogen response pathway and TGF-β signal are more strongly expressed in the IZ. And as shown in Fig. 2B, the allogeneic rejection, apical surface, the TGF-β signal, et.al are relevant to the high seizure frequency(HSF), on the contrary, the unfolded protein response, MYC-target, et.al are mostly expressed in patients with low-frequency seizures(LSF).

Fig. 2

A: Pathway differences between the EZ and IZ; B: Pathway differences between HSF and LSF. The vertical axis stands for the signal pathways, and the horizontal axis stands for the degree of expression differences in different subgroups (moderated-statistic). The pathways marked in blue and green (absolute value of GSVA>1) were statistically significant for specific clinical features.

Module-trait correlation

Two module-trait correlation analyses were performed. The modules presenting significant correlation values (r > |0. 50| and P < 0.05) were considered to be statistically correlated with the treat. One for clinical traits, (Fig. 3a) age at onset shows a negative correlation with the four modules in meta-module I and positive relation with three modules in meta-module II, respectively. The yellow module in meta-module I is negatively related to seizure frequency, whereas the pink module in meta-module II is positively related to that trait. Based on the above, the meta module II module is positively correlated with severe disease manifestations (early onset of disease and high frequency of attack), whereas in meta module I, these relationships are completely reversed.

Another for pathways, (Fig. 3B) the late stages of estrogen response is negatively correlated with the modules in meta-module I and positively relate to modules in meta-module II, and the early stages of estrogen response are quite inverted. Besides, the allogeneic rejection, the HSF-relevant pathway, is negatively correlated with the yellow module and positively correlated with the pink module. Inversely, the unfolded protein response and MYC-target-V1, which are highly expressed in patients with low seizure frequency, are negatively correlated with the pink module and positively correlated with the yellow module. All the relationships between modules and traits were displayed in Fig. 3E.

Screening results for key genes

The pink and yellow modules showed a statistically meaningful correlation with seizure frequency and frequency-relevant differential pathways simultaneously. Thus, the genes significantly correlated to seizure frequency and related pathways were screened in the pink and yellow modules. All gene values from the pink and yellow modules were depicted in a MM (x-axis) vs. GS graphic (y-axis). (Fig. 3C and D) The key genes with high GS values (the absolute value≥0.6) in both seizure frequency trait and relevant pathway trait include OSR2, CABP4, CAPSL, CYP4F8, and FRK in the pink module, and SH3GLB2, CHAC1, and DDX23 in the yellow module. The higher the absolute values of GS, the more biological importance is the gene. (Table 2) The related gene annotation by GO and KEGG analysis was concluded in Table 3.

Table 2

Gene-trait correlations for key genes with high GS value.

module	key genes	Frequency-GS	Pathway	Pathway-GS
pink	OSR2	0.65	allograft rejection	0.78
	CABP4	0.65	apical surface	0.68
	CAPSL	0.60		0.74
	CYP4F8	0.75	unfolded protein response	-0.67
	CABP4	0.65		-0.68
	FRK	0.63		-0.66
	CAPSL	0.60		-0.73
yellow	SH3GLB2	-0.64		0.82
	CHAC1	-0.63		0.81
	DDX23	-0.62		0.86

Note. The third and the rightmost column showed the GS value associated with the seizure frequency and the relevant pathway respectively.

Table 3

The key gene annotation.

module	key genes	GO	KEGG
pink	OSR2	transcription by RNA polymerase IIregulation of transcription by RNA polymerase II	SNARE interactions in vesicular transport
	CABP4	–	–
	CAPSL	calmodulin-related	–
	CYP4F8	Oxygenase	Arachidonic acid metabolism， Metabolic pathways
	FRK	innate immune responsetransmembrane receptor protein tyrosine kinase signaling pathway	MAPK signaling pathwaypathogen interaction
yellow	SH3GLB2	membrane organization，cellular developmental processmulti-organism cellular processvesicle targeting	Endocytosis
	CHAC1	sulfur compound catabolic processglutathione metabolic processorganonitrogen compound catabolic process	–
	DDX23	mRNA splicing, ATP hydrolysis activity	Spliceosome

Note. The GO column listed the molecular function (MF), biological process (BP) and cellular components(CC) of key genes, the KEGG column listed the signaling pathways where the genes are located. The table with “-” indicated that none of corresponding gene annotation found in

the analysis.

Discussion

Among the 30–40% of patients who have seizures that persist after taking medication, more are opting for epilepsy surgery. (Rugg-Gunn F et al., 2020) Nevertheless, two-thirds can relapse within five years after surgery (Petrik S et al., 2021). Therefore, surgical excision alone is not sufficient for curing refractory epilepsy. We should draw more attention to exploring the molecular mechanisms of epilepsy. Always, the inconsistency between EZ and IZ arouses debate. In the study, none of the modules presented a vital relation(r > |0.50|,P < 0.05) with the resection zone phenotype(EZ and IZ) in the module-trait correlation analysis, but the GSVA screened the pathway expression between EZ and IZ.

The IL6-JAK-STAT3-signal pathway promotes the expression of various acute-phase proteins and triggers a range of inflammatory and immune responses (Wang SW and Sun YM, 2014). The previous research strives to illustrate mechanisms of potential TGF-β signal-induced epileptogenesis of mice. Levy. N et al. indicated TGF-β signal upregulated IL-6 in mice brains during the period of epileptogenesis (Levy N et al., 2015). These two inflammatory responses are the most relevant pathways in the EZ and IZ, respectively. In recent decades, many researchers have reported a correlation between the inflammatory response and epilepsy, (Vezzani A et al., 2011) The result mentioned above is also consistent with previous reports. Moreover, it may note that the inflammatory response in EZ and IZ may arise from different mechanisms.

Of particular interest, early and late estrogen responses were highly expressed in the irritative and epileptic regions, respectively. Studies have shown that estrogen enhances epilepsy susceptibility through intracellular receptor binding and ion channel regulation. Binding to intracellular receptors can regulate gene expression and thus interfere with protein synthesis pathways to affect neuronal excitability [Stevens SJ and Harden CL 2011]. In addition, another approach is to activate Ca ion channels and promote the excitability of pyramidal cells [Navis A andd Harden C 2016] by directly affecting ion channels, including increasing the density of NMDA receptors on dendritic spines of hippocampal CAl pyramidal cells and Purkinje neurons in the cerebellum. Or lead to down-regulation of glutamate decarboxylase, which mediates the conversion of glutamate to GABA, and reduced the production of GABA (Verrotti A et al., 2012). The estrogen responses provide a new direction for exploring the relationship between the EZ and IZ, but the detail of that needs subsequent studies.

There are some modules presenting a significant correlation with seizure frequency, so we compared pathway expression between the high seizure frequency and low seizure frequency(HSF and LSF). The allogeneic rejection, apical surface, and the TGF-β signal, et.al are relevant to HSF. The inhibition of allograft rejection has been shown to reduce neuronal apoptosis following seizures (Liu. J et al., 2018), which suggests that immune response is involved in epileptogenesis. Besides, it is worthy to note that the TGF-β signal is also relevant to the IZ. The production of the irritative zone and high frequency of seizures may derive from similar pathogenesis.

On the contrary, the unfolded protein response (UPR) are mostly expressed in patients with low-frequency seizures. UPR is a protective mechanism elicited by misfolded proteins accumulating in the endoplasmic reticulum(ER) after body injury (Howell, SH, 2017). A large number of evidence indicates that endoplasmic reticulum stress may be involved in various pathological processes of epilepsy (Fu J et al., 2020). The UPR involves two heat shock proteins (HSP), HSPH4 and HSPA5. Nowakowska, Marta et al. found that persistent epilepsy (SE) induced up-regulated HSPA5 in the piriform cortex as well as down-regulated HSPA5 and HSPH4 in the hippocampus (Nowakowska M et al., 2020). Taken together suggests that UPR, the LSF-relevant pathway, may be initiated to exhibit protective effects against the ER stress in the epileptic pathologistic process and the related heat shock proteins may be the pharmacological regulatory targets of epilepsy.

In addition, the module-trait correlation analysis showed that the pink module is positively correlated with seizure frequency, while the yellow module is negatively correlated with seizure frequency. Meanwhile, a strong relationship also exists between the two modules and frequency-correlated pathways. Therefore, the study selected the key genes with high GS values that were highly connected with the seizure frequency at the clinical and pathway level in the pink and yellow modules. The genes and the signal pathways listed may respectively provide the upstream and downstream markers of neocortical refractory epilepsy.

Ultimately, the upregulated genes associated with seizure frequency include OSR2, CABP4, CAPSL, CYP4F8, and FRK in the pink module, and the downregulated genes include SH3GLB2, CHAC1 and DDX23 in the yellow module. Among these genes, CABP4, as well as DDX23, have been reported to be involved in epileptic processes, being potential biomarker candidates. CaBP4 (Ca2 +-binding protein 4) encodes the regulation of voltage-gated Ca(v)1 Ca(2 +) channels, which may be a sort of apical surface. The mutation in the CABP4 was predicted to affect protein function and may be connected to autosomal dominant nocturnal frontal lobe epilepsy (ADNFLE) (Chen ZH et al., 2017). Missense alterations in DDX23 (Dead-Box Helicase 23) can manifest as seizures in some patients, gathered through GeneMatcher (Burns W et al., 2021). As shown in Table 3, DDX23 is connected with mRNA splicing and ATP hydrolysis activity.

The gene CYP4F8 with the highest GS value encodes for proteins belonging to the cytochrome P450 (CYP) 4 family enzymes (Liu J et al., 2021), which is distributed in the ER. And the relationship between ER stress and UPR has been mentioned above.

It participates in regulating several biological processes, including cholesterol, lipids and steroid synthesis, and arachidonic acid (AA) metabolism (Edson, K and Rettie, A.2013). CYP4F8 is associated with arachidonic acid metabolism and metabolic pathways according to KEGG analysis, in line with previous reports. Certain pro-soluble lipids produced by AA-activated lipoxygenases (LOXs) can control the inflammatory response, (Serhan CN, 2017; Serhan CN et al., 2008) and the pro-soluble response in the hippocampus is out of balance during epilepsy (Frigerio, F et al., 2018]. Based on these findings, the upregulation of CYP4F8 plays an important role in high seizure frequency as its overexpression was also identified in our study.

As for other upregulated genes, CAPSL encode calcyphosin-like protein and is involved in calcium modulation reactions. OSR2, odd-skipped related transcription factor 2, participates in the regulation of transcription by RNA polymerase. It is closely associated with immune infiltration, (Ma T et al., 2022) as expected, is relatively higher expressed in the allograft rejection. Also, OSR2 downregulation can diminish TGF-β signaling (Anh LPH et al., 2022). OSR2 upregulation may relate to the inflammation response of epilepsy. FRK, UPR-relevant gene, is correlated to transmembrane receptor protein tyrosine kinase signaling pathway.

Other genes that are negatively associated with seizure frequency include SH3GLB2 and CHAC1. SH3GLB2, which encodes the endophilin B2, is responsible to regulate the trafficking of endocytic vesicles and autophagosomes to late endosomes or lysosomes.(Serfass JM et al., 2017) This is consistent with GO and KEGG analysis. CHAC1, Cation transport regulator-like protein 1, was found in a co- regulated group of genes enriched for components of the ATF4 (activating transcription factor 4) arm of the unfolded protein response, (Mungrue IN et al., 2009) which has the significant correlation with the seizure frequency.

In conclusion, by adopting the GSVA and the WGCNA approach and integrating clinical and transcriptomic data, we unveiled the inconsistency in biological pathway expression of EZ and IZ, and the differences in signaling pathways and targets may explain the inconsistency in brain locations between them in the clinical pattern, we also found several key genes and signal pathways significantly associated with the seizure frequency. Above all may provide a novel direction for the potential molecular mechanisms for refractory epilepsy phenotypes and can be targeted for new preventive as well as drug-based therapy interventions.

Conclusion

The occurrence of EZ and IZ act on the different biological mechanism. Of particular interest, early and late stages of estrogen responses are highly expressed in the IZ and EZ, respectively. The upregulated genes associated with seizure frequency include OSR2, CABP4, CAPSL, CYP4F8, and FRK, and the downregulated genes include SH3GLB2, CHAC1 and DDX23.

Study Limitations

There are still limitations in this study, including the small sample size and the lack of experimental verifications. Besides, the EZ and IZ in the original data were not located by SEEG. Thus, further study will require more samples localized by SEEG for validation.

CRediT authorship contribution statement

Rui Zhang: Conceptualization, Methodology, Investigation, Data analysis, Writing – original draft. Yan Chen: Data curation, Formal analysis. Jia He: Data curation, Formal analysis. Hai-yan Gou: Investigation. Yu-lan Zhu: Investigation. Yan-mei Zhu: Conceptualization, Supervision, Writing – review & editing.