Profiling Analysis of Circular RNA and mRNA in Human Temporal Lobe Epilepsy with Hippocampal Sclerosis ILAE Type 1

Hippocampal sclerosis (HS) is the most common surgical pathology associated with temporal lobe epilepsy (TLE). However, the cause of TLE with or without HS remains unknown. Our current study aimed to illustrate the essential molecular mechanism that is potentially involved in the pathogenesis of TLE-HS and to shed light on the transcriptional changes associated with hippocampal sclerosis. Compared to no-HS group, 341 mRNA transcripts and 131 circRNA transcripts were differentially expressed in ILAE type 1 group. The raw sequencing data have been deposited into sequence-read archive (SRA) database under accession number PRJNA699348.Gene Ontology analysis demonstrated that the dysregulated genes were associated with the biological processes of vesicle-mediated transport. Enrichment analysis demonstrated that dysregulated genes were involved mainly in the MAPK signal pathway. Subsequently, A total of 441 known or predicted interactions were formed among DEGs, and the most important module was detected in the PPI network using the MCODE plug-in. There were mainly four functional modules enriched: ER to Golgi transport vesicle membrane, Basal transcription factors, GABA-gated chloride ion channel activity, CENP-A containing nucleosome assembly. A circRNA-mRNA co-expression network was constructed including 5 circRNAs(hsa_circ_0025349, hsa_circ_0002405, hsa_circ_0004805, hsa_circ_0032254, and hsa_circ_0032875) and three mRNAs (FYN, SELENBP1, and GRIPAP1) based on the normalized mRNA signal intensities. This is the first to report the circRNAs and mRNAs expression profile of surgically resected hippocampal tissues from TLE patients of ILAE-1 and no-HS, and these results may provide new insight into the transcriptional changes associated with this pathology.


Introduction
Circular RNAs (circRNAs) are a unique type of non-coding RNAs. Unlike normal linear RNAs, their 3'-end of an exon is spliced to the 5'-end of an upstream exon resulting in a circular RNA molecular, which makes them more stable and conserved. CircRNAs are known to be able to interact with many molecules. A single can bind to several microRNAs (miRNAs) and suppress downstream target genes, resulting in gene suppression. This inhibition function is called a "sponge effect". Recently, several circRNAs have been reported to be translated into protein whether exogenous Gene bank: The raw sequencing data have been deposited into sequence read archive (SRA) database under accession number PRJNA699348.
1 3 circRNAs or endogenous circRNAs. In addition, circRNAs can directly bind to RNA-binding proteins and affect the expression and function of their target genes (Zang et al. 2020). With the rapid development of high-throughput sequencing, thousands of circRNAs are expressed at high levels in the brain (Meng et al. 2019). In recent years, several studies focused on differentially expressed circRNAs and their function in epilepsy patients and animal models, suggesting circRNAs as potentially novel biomarkers or therapeutic targets in this disease.
Epilepsy is one of the most common chronic neurological diseases, affecting more than 65 million people in the world. Epilepsy is a condition characterized by recurrent epileptic seizures, and causes impairment in cognition, memory, quality of life, and social psychology. Over the past three decades, the use of more than 12 third-generation antiepileptic drugs has provided doctors and patients with more options for treating multiple seizures. However, almost onethird of patients with epilepsy have poor control of seizures after drug treatment (Wiebe and Jette 2012). Hippocampal sclerosis (HS) is the most common type of histopathologic abnormality in patients with drug-resistant temporal lobe epilepsy (Blümcke et al. 2013). The International League Against Epilepsy (ILAE) put forward a classification of HS according to histopathologic abnormality and three subtypes were identified. The ILAE HS type 1 is the most common type of HS with 60-80% of all patients with TLE-HS (Blumcke et al. 2013). Our previous research found that the ILAE HS type 1 has a favorable prognosis (Na et al. 2015). However, the specific molecular mechanism of different prognosis in HS type 1 patients is still not clear.
In this study, we examined the expression profile of circR-NAs in HS ILAE type 1 and no-HS type patients with TLE via next-generation sequencing analysis in a Chinese population. High-throughput data gave us novel insights that circR-NAs were differentially expressed, and several key circRNAs showed potential as candidate makers for TLE-HS. RT-PCR analysis was performed to confirm the results of the RNA-seq and for further validation. Then, we performed an integrated analysis of circRNAs and mRNAs expression to discover novel key components and pathways involved in TLE-HS. A co-expression network was constructed to identify the circR-NAs and mRNA interaction. This study may provide molecular information based on pathogenesis in Type 1 TLE patients.

Patient Selection and Tissue Collection
All samples were obtained in the First Affiliated Hospital of Harbin Medical University (Heilongjiang, China) during anteromesial temporal resections. The study was approved by the Ethics Committee Board in this hospital. All samples were confirmed by immunohistochemistry according to the ILAE criteria (Blumcke et al. 2013). Full epileptological assessment including 3.0 T magnetic resonance imaging of the hippocampus, video-electroencephalographic monitoring, electroencephalogram, cognitive and neuropsychological testing, and interictal fluorodeoxyglucose-positron emission tomography combined with computed tomography. Selective amygdalohippocampectomies or standard anterior temporal lobectomies with amygdalohippocampectomies were performed on the TLE patients by the same surgeon. The excised hippocampal specimen was divided into two parts along the long axis in the operating room. One part was put in the liquid nitrogen immediately, the other part was fixed in 4% paraformaldehyde overnight at 4 °C and then dehydrated and embedded in paraffin. The body of the hippocampus was resected in each sample for next-generation sequencing and RT-PCR. A total of 7 ILAE-1 and 7 no-HS samples were collected. We randomly chose three ILAE-1 samples and three no-HS samples and prepared them for RNA-Seq. Some samples were collected in our previous study investing miRNAs . The clinical characteristics of these patients are listed in Table S1.

Immunohistochemistry
Tissue slices were taken at 4 μm from the paraffin-embedded hippocampal tissue. The steps are the same as those mentioned in our previous article . Briefly, the slices were dewaxed and hydrated, and treated with 3% hydrogen peroxide. Then antigen repair was performed in sections with high pressure. The slices were incubated in NeuN monoclonal antibody (Zhongshan Golden Bridge Biotechnology, Beijing, People's Republic of China) for 1 h. Finally, DAB was used as a chromogen ( Figure S1).

RNA Isolation and Quality Control
Total RNA was extracted from the frozen hippocampal tissue using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA) according to the instructions. The purity and quantity of RNA were determined using a Nano-Photometer® spectrophotometer (IMPLEN, CA, USA) and a Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

RNA Library Preparation and Sequencing
A total amount of 5 μg RNA per sample was used as input material for the RNA sample preparations. First, Total RNA was treated with the Epicentre Ribozero™ rRNA Removal Kit (Epicentre, USA), to remove all the rRNAs. The remaining RNAs were processed using the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer's recommendations. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia). After cluster generation, the libraries were sequenced on an Illumina Hiseq 4000 platform and 150 bp paired-end reads were generated. The raw sequencing data have been deposited into sequence read archive (SRA) database under accession number PRJNA699348.

Sequence Data Analysis
Clean reads were obtained by removing reads containing adapter, reads on containing ploy-N and low-quality reads from raw data. All the downstream analyses were based on clean reads with high quality. Next, the clean reads were mapped to the reference sequence using Bowtie (Langmead and Salzberg 2012). The circRNA were detected and identified using find_circ (Memczak et al. 2013) and CIRI2 (Gao et al. 2018). The raw counts were first normalized using TPM.
Differential expression analysis of two groups was performed using the DESeq R package (1.10.1). The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P-value found by DESeq were assigned as differentially expressed. An adjusted p value of < 0.05 was considered statistically significant.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Pathway Analysis
Gene Ontology (GO) enrichment analysis for genes of differentially expressed were implemented by the GOseq R package, in which gene length bias was corrected (Young et al. 2010). GO terms with corrected P value less than 0.05 were considered significantly enriched by differential expressed genes. KEGG is a database resource for understanding highlevel functions and utilities of the biological system (http:// www. genome. jp/ kegg/). We used KOBAS software to evaluate the statistical enrichment of differential expression genes in KEGG pathways (Mao et al. 2005).

Integration of Protein-Protein Interaction (PPI) Network and Module Analysis
The STRING database (http:// string-db. org) is an online network that critically assesses protein-protein interactions (PPI).

Normalized expression level = Mapped reads * 1000000∕Total reads
The present study inputs the DEGs on STRING in order to obtain the PPI network from it (confidence score > 0.400). Subsequently, the PPI network was constructed using Cytoscape software (http:// www. cytos cape. org). Cytoscape is a type of software for visualizing molecular interaction networks and biological pathways, with annotations, gene expression profiles and other datasets (Shannon et al. 2003). The present study used 4 plugins of Cytoscape: MCODE, CytoHubba, ClueGo, and CluePedia. MCODE is a Cytoscape application that clusters a given network based on topology in order to produce highly interconnected subgraphs of molecular complexes, and parts of pathways (Bader and Hogue 2003). Then, we used the CytoHubba plugin to calculate the degree of each node (Chin et al. 2014). Twelve topological algorithms(MCC, DMNC, MNC, Degree, EPC, BottleNeck, EcCentricity, Closeness, Radiality, Betweenness, Stress, ClusteringCoefficient) were carried out to identify the hub genes. ClueGO is a Cytoscape plug-in that can divide large clusters of genes into functional groups based on GO, KEGG, Wiki Pathways and Reactome (Bindea et al. 2009). CluePedia is a ClueGO plugin that can integrate information regarding genes, proteins into a network with ClueGO terms/pathways (Bindea et al. 2013).
In the present study, MCODE was used to create clusters. The main parameters were as follows: Node Density Cutoff: 0.1; node score cutoff: 0.2; Degree Cutoff: 2; K-Core: 2; and Maximum Depth: 100. Subsequently, the ClueGo and Cluepedia were used to analyze the clusters and visualize the functional groups based on GO and KEGG. The evidence codes used were as follows: #Genes in KEGG, GO_MF, GO_CC and GO_BP; P < 0.05.

Quantitative Real-Time Polymerase Chain Reaction (RT-PCR) Validation
Total isolated RNAs from ILAE-1 and no-HS groups were reversely transcribed using a EZNA™ Total RNA kit II (Omega Bio-Tek) and single-stranded cDNA was synthesized by Reverse Transcription kit (Toyobo Life Science, Japan). GAPDH was amplified as internal control and the relative amount of each circRNA and mRNA to GAPDH was calculated using the Eq. 2 − ΔCT. In detail, based on combinational consideration of the fold change, raw data, p value, 5 circRNAs and 3 mRNAs were selected for RT-PCR verification by a SYBR Green Real-Time PCR Master Mix (Toyobo Life Science, Japan) with 3 replication each. The sequences of primers are listed in Table S2.

Construction of the circRNA-mRNA Co-Expression Network
To identify the co-expressed mRNAs with DE-circRNAs, we calculated the correlation coefficient between DE-circRNAs 1 3 and DE-mRNAs at expression levels using the Pearson correlation test. The raw P-values () were corrected by multiple hypotheses using a permutation approach (Yao et al. 2015): keeping the circRNA expression constant, randomly disturbing the mRNA expression 1,000 times, and recalculating the Pearson correlation between the circRNA expression and disturbed mRNA expression. The permutation P-values () were used to calculate the empirical P-values () by the following formula (which introduces a pseudo-count of (1): where num P p < P r represents the number of P-values less than P r in the P p . The circRNA-mRNA pairs with P e < 0.01 were considered as co-expression relationships.
On the basis of the co-expressed circRNA-mRNA pairs, we constructed a circRNA-mRNA co-expression network. Each circRNA or mRNA represents a node, and two RNAs were connected by an edge, indicating a close correlation. Next, we identified the hub nodes by calculating the degree of each node (d), which reflected the importance of the node in the network. The normalized degree of each node was computed (Hu et al. 2019), as follows: where D i denotes the normalized degree of the ith node, d i represents the degree of the ith node, and max(d) indicates the maximum degree in the network. The higher the D value was, the more important role the node in the network played. We considered the nodes with D > 0.6 was the hub nodes in the network. The degree of each node was calculated using the "degree" function in the R package "igraph".

Differential Expression of CircRNAs and mRNAs Between ILAE-1 and no-HS Tissues
To investigate molecular mechanisms in TLE-HS, hippocampus tissues obtained from six patients were tested and total RNA was isolated for whole transcriptome analysis with RNA sequencing. The clean reads were mapped to the human reference genome using TOPHAT and were spliced into putative transcripts. Bioinformatics analysis was performed to identify differentially expressed mRNAs and circRNAs between the HS ILAE type 1 and no-HS groups. Gene expression variation was visualized by the volcano plots, as shown in Fig. 1. Every (1) P e = num P p < P r + 1 1001 , D ∈ (0, 1), point represents an mRNA or circRNA. 341 mRNAs (229 upregulated and 112 downregulated) and 131 circRNA (94 upregulated and 37 downregulated) were identified between two groups using fold-change filtering (|log2(fold change)|> 1) and p value < 0.05. The most marked 10 upregulated and downregulated circRNAs of ILAE-1 compared with no-HS from our NGS results were listed in Table 1. The hierarchical cluster showed the overview of gene expression. The result showed that the expression of circRNAs and mRNAs were different between these two groups. The depth of colour represents the expression level of mRNAs and circRNAs; green represents and low red represents high (Fig. 2).

Function Analysis of Differentially Expressed mRNAs Involving ILAE-1 and no-HS Tissues
We performed GO and KEGG pathway analysis of differentially expressed mRNAs between these two groups.
The top 30 GO terms were performed based on three categories: the cellular component (CC), molecular function (MF), and biological process (BP). For the BP group, the most meaningful GO terms were vesicle-mediated transport ( Fig. 2A). For the MF group, the main represented GO term was protein binding (Fig. 2B). For the CC group, the main represented category was cytosol (Fig. 2C). The KEGG pathways were listed in the Fig. 2D, and the pathway of MAPK signaling contained most of genes.
To further elucidate the functional relationship among different expression genes(DEGs), a protein-protein interaction network was generated with the STRING online database using a combined scoring method. A total of 441 known or predicted interactions were formed among 211 DEGs(PPI enrichment P-value = 7.83e-08), and the detailed results of PPI network are provided in supplement file. Subsequently, the PPI network was constructed using Cytoscape software (Fig. 3A). The hub genes in the PPI network were selected using the cytohubba plugin. We identified a total of four hub genes, including UBE3A, NCAM1, LMO7, NTRK2 and FYN, with UBE3A, NCAM1 and FYN being highly expressed and the rest being lowly expressed in HS ILAE type 1 group. Finally, top five overlapping genes were obtained. MCODE plugin was used to identify the modules in the network; modules including at least 10 nodes were selected (Fig. 3B-C). Further, GO and KEGG analyses of the modules were demonstrated by ClueGo/CluePedia plugin. There were mainly four functional modules enriched: ER to Golgi transport vesicle membrane, GABA-gated chloride ion channel activity, CENP-A containing nucleosome assembly, Basal transcription factors ( Figure S2).

Construction of a Co-Expression Network
To date, the functions of most circRNAs have not been determined. Therefore, we constructed a circRNA-mRNA co-expression network to identify the critical circRNA in HS ILAE type 1 (Fig. 4). The candidates to validate were chosen based on the following criteria: |log2(fold change)|> 4and P-value < 0.01. A total of 14 DE-circR-NAs (4 down-regulated and 10 up-regulated) and 42 DE-mRNAs (8 down-regulated and 34 up-regulated) were selected. There were seven hub nodes were identified in the co-expressed network, which included five circRNAs (hsa_circ_0025349, hsa_circ_0002405, hsa_circ_0004805, hsa_circ_0032254, and hsa_circ_0032875) and three mRNAs (FYN, SELENBP1, and GRIPAP1).

Validation of Dysregulated circRNAs and the Corresponding mRNAs Between ILAE-1 and no-HS Tissues
To validate the results of the RNA-seq, we chose a total of five circRNAs (hsa_circ_0025349, hsa_circ_0002405, hsa_circ_0004805, hsa_circ_0032254, and hsa_ circ_0032875) and three mRNAs (FYN, SELENBP1, and GRIPAP1) for further RT-PCR. The expression of   (Fig. 5). For all eight selected targets, the expression results calculated by the RT-PCR analysis were identical to those found by RNA-Seq, confirming the reliability of the RNA-Seq results.

Discussion
In adults, TLE is the most common form of drug-refractory epilepsy and surgery is the most effective therapy for this condition. Hippocampus in these patients often shows hippocampal sclerosis, a hallmark pathophysiological abnormality, which is characterized by segmental pyramidal cell loss and astrogliosis. In 2013, the ILAE developed a consensus classification system based on the qualitative histopathologic assessment of hippocampal subfields. The classification includes patients with ILAE HS type 1 (the most common type, neuronal loss in CA1 and CA4 predominantly), ILAE HS type 2 (neuronal loss in CA1) and ILAE HS type 3 (neuronal loss in CA4) (Blumcke et al. 2013). In our previous research (Na et al. 2015), we identified that ILAE HS type 1 has the best outcome, supporting the notion that TLE-HS is a heterogeneous condition. By focusing particularly on HS ILAE type 1, the dysregulated mRNA and circRNAs specific to this subtype could be evaluated with higher resolution than possible with a previous study of TLE-HS.
In this study, we conducted RNA sequencing of hippocampal tissue from HS ILAE-1 and no-HS groups to characterized the molecular mechanism related to HS. Primarily, the differential expression analysis identified 341 mRNA transcripts (229 upregulated and 112 downregulated) and 131 circRNA transcripts (94 upregulated and 37 downregulated) between these two groups. GO and KEGG pathway analyses were carried out to analyze the biological functions of the dysregulated mRNAs. Our GO results suggested that the differentially expressed genes were mostly enriched in vesicle-mediated transport protein binding and cytosol. Our KEGG pathway analysis showed that the mitogenactivated protein kinase (MAPK) signal pathway was the most enriched and the prion disease pathway was the most significantly altered. MAPKs are important signal transducers for many kinds of stimulation inside and outside the cell. MAPKs have been found to play a crucial role in synaptic plasticity along with neuronal damage and survival in neurodegenerative disease (Correa and Eales 2012). There is some evidence that MAPK plays a vital role in mediating the toxic effects of prion on synaptic structure and function in the neuronal system (Fang et al. 2018). We also mapped the differential protein interaction using the STRING database. Visualization of the network was generated by Cytoscape and MCODE was used to detect the gene modules in the PPI network. ClueGO analysis of network shows genes binned into ER to Golgi transport vesicle membrane, GABA-gated chloride ion channel activity, CENP-A containing nucleosome assembly, and Basal transcription factors pathway. Five hub genes were selected from the PPI network including UBE3A, NCAM1, LMO7, NTRK2, and FYN. UBE3A is an E3 ubiquitin ligase. Deletion or replication can lead to different neurodevelopmental diseases such as Angelman. UBE3A is highly expressed in immature hippocampal neurons in growth cones (Trezza et al. 2019), while growth cones have a major role in neuronal guidance. Therefore this may be related to aberrant neurogenesis in epilepsy (Lybrand et al. 2021). NCAM1 is a member of the cell adhesion molecules, mainly expressed in the nervous system, involved in the regulation of neuronal function and playing a key role Fig. 4 The correlation between the detected circRNAs and mRNAs. A Network showing the co-expressed relationships between circR-NAs and mRNAs. The yellow nodes represent the circRNAs. The blue nodes represent the mRNAs. The diamond nodes represent the hub nodes in the network. B The Pearson correlation matrix between the hub nodes. The size and color intensity of the circle denotes the correlation coefficient value. Positive correlations are visualized in red and negative correlations in a blue scale. P-value is marked in the matrix if P-value is below 0.05, which represents statistical significance 1 3 in synaptic plasticity. NCAM-1 in cerebrospinal fluid has been reported as a potential biomarker for drug-refractory epilepsy (Wang et al. 2012). As an emerin-binding protein, LMO7 has been little studied with epilepsy. However, LMO7 is consistently highly expressed in reactive oxygen species-resistant neuronal cells, exhibiting some neuronal protective effects (Maczurek et al. 2013). And the continuous release of reactive oxygen species is a key factor leading to neuronal deficiency in hippocampal sclerosis (Williams et al. 2015), which explains that LMO7 is lowly expressed in the hippocampal sclerosis group. NTRK2 encodes a member of the neurotrophic tyrosine receptor kinase family. Almoguera et.al identified NTRK2 as a susceptibility genes in drug-resistant epilepsy (Almoguera et al. 2019). NTRK2 is highly expressed in granule cells and is an important regulator of granule cell morphology. NTRK2 shapes hippocampal circuits by controlling early GABA signaling (Badurek et al. 2020), and low NTRK2 expression may exacerbate hippocampal circuit disorders following hippocampal sclerosis in epileptogenesis.
CircRNA can bind to mRNA and regulate gene expression at post-transcriptional levels. We assume that a combined analysis of circRNA and mRNA could be useful to evaluate the impact of circRNAs on epilepsy. We chose five circRNAs and three corresponding mRNAs in this network for further validation. RT-PCR results revealed that selected dysregulated circRNAs and mRNAs have distinct expressions between these two groups. The expression of hsa_circ_0032254 and hsa_circ_0032875 was correlated with FYN, SELENBP1 and GRIPAP1. The hsa_ circ_0032254 is related to gastric cancer (Bu et al. 2020). However, there is no report in epilepsy. The parental gene of hsa_circ_0032254, known as GPHN, caught our attention. The GPHN encodes a neuronal assembly protein, Gephyrin. In TLE specimens and rat models, Gephyrin was found low expressing. Previous studies have shown that mutation in GPHN is a risk factor for epilepsy (Tyagarajan and Fritschy 2014). A possible mechanism is that the properly folded and stably expressed mutant gephyrin possibly can interact with all its binding partners. Therefore, it has the capability to reduce the clustering and aggregation of a variety of postsynaptic proteins other than GABAARs and GlyRs, thereby serving dominant-negatively on downstream pathways. This would most likely impair the plasticity and modulation of gephyrin scaffolds, which are critical for the development of an inhibitory neuronal circuit (Dejanovic et al. 2017). EML5 splicing hsa_circ_0032875 and is a microtubule-associated protein that is expressed in the central nervous system. Previous researches showed that EML5 is highly expressed in pharmacoresistant epilepsy and may be related to the extension of neuronal dendrites and axons (Sun et al. 2015). It was assumed that circRNAs and their parent mRNAs are homologous and could have the same biological function.
As far as mRNA expression is concerned, three candidate mRNAs were found to be differentiated between ILAE-1 and non-HS groups. The Fyn, a non-receptor Src family of tyrosine kinase, acting as potential mediators of neuroinflammation (Tamura et al. 2001). It is reported that the Fyn was upregulated in the microglia and oligodendrocytes in the hippocampus during epileptogenesis. Silencing Fyn will suppress epileptogenesis through reduce neuroinflammatiom and apoptosis (Sharma et al. 2018;Luo et al. 2020). All these results indicate that SELENBP1 and GRIPAP1 in HS ILAE type 1 may affect the regeneration of neurons compared with those in no-HS type; however, this aspect of the study has not been reported yet.
The present study was the first to characterize the mRNA and circRNA expression profile in TLE patients between ILAE-1 and no-HS groups. However, several limitations need to be addressed. Firstly, the samples enrolled in this study were limited and we need a large number of subjects to verify our results and further decrease the scale of cir-cRNA profile as a diagnostic biomarker. Secondly, although several specific mRNAs and circRNAs had been discovered, their detailed mechanism in the progression of TLE and the relationship with prognosis was not investigated. Further experiments should be carried out. Lastly, the absence of patients with HS ILAE type 2 and 3, which should have been compared in this study. Our team is collecting these two rare types of specimens for further study.