Identification of Specific Circular RNA Expression Patterns and MicroRNA Interaction Networks in Mesial Temporal Lobe Epilepsy

Circular RNAs (circRNAs) regulate mRNA translation by binding to microRNAs (miRNAs), and their expression is altered in diverse disorders, including cancer, cardiovascular disease, and Parkinson’s disease. Here, we compare circRNA expression patterns in the temporal cortex and hippocampus of patients with pharmacoresistant mesial temporal lobe epilepsy (MTLE) and healthy controls. Nine circRNAs showed significant differential expression, including circRNA-HOMER1, which is expressed in synapses. Further, we identified miRNA binding sites within the sequences of differentially expressed (DE) circRNAs; expression levels of mRNAs correlated with changes in complementary miRNAs. Gene set enrichment analysis of mRNA targets revealed functions in heterocyclic compound binding, regulation of transcription, and signal transduction, which maintain the structure and function of hippocampal neurons. The circRNA–miRNA–mRNA interaction networks illuminate the molecular changes in MTLE, which may be pathogenic or an effect of the disease or treatments and suggests that DE circRNAs and associated miRNAs may be novel therapeutic targets.


INTRODUCTION
Circular RNAs (circRNAs) are a unique class of non-coding RNA found in human tissues and are linked to diverse diseases (Chen et al., 2016(Chen et al., , 2018. CircRNAs originate from backsplicing of the precursor-mRNA (pre-mRNA) transcript. This alternative (back) splicing of pre-mRNA transcripts occurs when the spliceosome covalently joins the 3 end of a downstream exon to the 5 end of an upstream exon. Thus, circRNAs are covalently closed and lack polyadenylated 3 -tails or 5 -caps, making them exonuclease resistant and longer lived (∼48 h) than linear RNA (∼10 h; Curry-Hyde et al., 2020). The stability of circRNAs allows for packaging in extracellular vesicles, suggesting a role in cell-to-cell communication (Lasda and Parker, 2016). circRNAs may also act as microRNA (miRNA) sponges; that is, some circRNA sequences are complementary to miRNAs and may sequester these miRNAs, preventing them from binding to target mRNA (Hansen et al., 2013;Jeck et al., 2013).
Mesial temporal lobe epilepsy (MTLE) is a common human epilepsy and is often accompanied by hippocampal sclerosis (HS; Engel, 2001): a pathology associated with pharmacoresistance (Blümcke et al., 2013). HS is characterized by neuronal loss and gliosis, granule cell dispersion, and aberrant mossy fiber sprouting (axonal projection) maximal in the CA1 and dentate gyrus. The hippocampus is critical for short-term memory (Bauman et al., 2019). MTLE patients with HS (MTLE-HS), especially in the dominant (or both) hemisphere(s), exhibit memory deficits. Pharmacoresistant epilepsy and memory impairment may result from cell loss, synaptic reorganization, altered excitatory-inhibitory balance, and aberrant hippocampal discharges (Scharfman et al., 2003;Blümcke et al., 2009;Gelinas et al., 2016). Many MTLE patients undergo surgical resection if seizures are pharmacoresistant, with seizure-free rates of 60-80%, and reduced seizures in 95%. However, many surgical candidates do not undergo epilepsy surgery due to limited health care resources or fears of complications (Kang et al., 2016).
We recently reported linear RNA and small RNA transcriptome analysis in the temporal cortex and hippocampus of MTLE individuals. We found large transcriptomic changes for linear RNA in the hippocampus but not in the cortex. By contrast, differentially expressed (DE) small RNAs (e.g., miRNAs) were identified across both brain regions (Mills et al., 2020). Two other studies compare circRNAs in MTLE patients and controls, although brain regions from controls did not precisely match MTLE patients. The first, a microarray analysis, sampled a 5639 probe circRNA chip for human circRNA splicing sites. A > 3-fold DE was identified for 586 circRNAs; the top 10 dysregulated circRNAs were validated with real-time PCR (Gong et al., 2018). The second study found 254 DE down-regulated circRNAs and 188 DE up-regulated circRNAs (Li et al., 2018). These results may reflect tissue-specific differences in circRNA expression owing to the different source of tissues in MTLE patients versus controls.
Here, we investigated circRNA expression in the temporal cortex and hippocampus of MTLE and healthy individuals. Our results identify aberrant circRNA expression in MTLE individuals. We hypothesize that altered circRNA expression perturbs gene networks in MTLE by sequestering miRNAs from their target genes. Using bioinformatic techniques, we analyzed the miRNA sponging potential of DE circRNAs and assessed altered gene networks. To our knowledge, this is the first study to combine linear, circular, and miRNA transcriptome profiles in assessing MTLE molecular signals.

Data Acquisition
Our patient and control cohorts, tissue sampling, and RNA sequencing and analysis methodologies were reported in a previous study (Mills et al., 2020). This data set contains 14 postmortem control samples (6 temporal cortex and 8 hippocampal) and 24 surgical MTLE samples (7 cortical and 17 hippocampal). Demographic and clinical characteristics for each patient and control are summarized in Supplementary Table 1. Control tissues were obtained during the autopsy of age-matched individuals without neurological disease. The MTLE tissues were obtained at surgical resection. For total and small RNA-seq library preparations, ribosomal RNA (rRNA) was depleted, and paired-end sequencing was performed using HiSeq 4000 with read lengths of 151 nucleotides (nts) to a depth of 50 million and 20 million reads, respectively (Mills et al., 2020).

Read Quality Analysis
To assess the quality of sequencing, each forward and reverse read file was analyzed with FastQC (Andrews, 2018;version 0.11.8). Trimmomatic (Bolger et al., 2014;version 0.36) was used to remove low-quality reads and sequencing adaptors from each read. Trimmomatic was run in paired-end mode with a phred33 score. Low-quality leading and trailing reads of 3 nt were removed with a sliding window of 4:15 and a minimum length of 30 nt. After trimming, FastQC determined the quality of trimmed reads. The small RNA data sets were assessed with FastQC and sequencing reads trimmed with Trimmomatic, which was run in paired-end mode with a phred33 score. Low-quality leading and trailing reads of 3 nt were removed with a sliding window of 4:15 and minimum length of 17 nt.

Read Alignment
Sequencing reads were aligned to the human reference genome GRCh38 accessed from UCSC (Kent et al., 2002). Reads were mapped with STAR (Dobin et al., 2013;version 2.6.1b) to rapidly generate linear and chimeric reads and information pertaining to the uniquely mapped read quantity.

Linear Detection
To analyze linear and small RNA transcripts, the indexed.bam file generated by STAR was parsed to StringTie (Pertea et al., 2015;version 1.3.4d) with the GRCh38 Ensembl or miRbase annotation.gtf for the small RNA data sets. StringTie was run with the -e flag to identify an abundance of known linear transcripts and ignore novel transcripts. The StringTie Python script prepDE.py created a matrix file of each gene and its respective expression across all samples.

CircRNA Detection
We used the CIRCexplorer2 (Zhang et al., 2016;version 2.3.6) and DCC (Cheng et al., 2016;version 0.4.8) programs to detect circRNAs. Each program provides different degrees of functionality. For example, the output of CIRCexplorer2 provides information regarding the coding strand and how many exons compose a circRNA. Using DCC, we removed potential circRNAs with <2 reads and that were not expressed in at least 10 samples to reduce false positives.
With a custom Python script, the output files of these two programs were merged and filtered for counts per million mapped reads (CPM) values of ≥0.1. CPM normalizes circRNA abundance to account for the differences in the number of uniquely mapped reads across the samples. This filtration guarantees that the circRNAs detected by both tools are highly expressed and deemed as bona fide.

Analysis of Expression Distribution
To determine whether detected circRNAs and linear transcripts cluster with respect to tissue or disease type, we performed a t-distributed stochastic neighbor embedding (t-SNE) machine learning algorithm with the R (R Development Core Team, 2013) package Rtsne (Van Der Maaten et al., 2014).

Differential Expression Analysis
For the DE analysis, we used the R program edgeR (Robinson et al., 2009;Law et al., 2018). Two DE experiments were performed in which the control-cortex was compared to the MTLE-cortex and the control-hippocampus was compared to the MTLE-hippocampus. A circRNA or mRNA was removed if its read count was zero across ≥10 samples. DE analysis was carried out using the edgeR decideTests function; genes with a Benjamini-Hochberg adjusted p-value < 0.05 were considered statistically significant.

miRNA Response Elements Analysis
To identify miRNA response elements (MREs) within the DE circRNAs, the ENCORI (Li et al., 2014;Dori and Bicciato, 2019) Argonaut (AGO)-CLIP-seq database was queried. A custom Python script searched for miRNA binding sites located within the specific circRNA transcript coordinates. Using a curl command to access the ENCORI API, we collected the miRNA-mRNA interaction data for all miRNAs binding to each circRNA. Within this command, we filtered the output for ≥1 supporting CLIP-seq experiments and ≥2 miRNA targetpredicting programs. A custom Python script collected the read coverage determined by StringTie for each target gene's transcript for each sample.
We employed edgeR to identify DE transcript sets for circRNAs of interest. For example, if the circRNA was DE in the hippocampus, we performed a linear DE of the target genes of the miRNAs predicted to bind that circRNA between the controlhippocampus and the MTLE-hippocampus. The DE matrix was input to ggplot (Wickham, 2006), ggrepel (Slowikowski, 2017), and dplyr (Wickham et al., 2017) to generate volcano plots to visualize DE patterns. To convert ensemble transcript identifiers to more readable gene symbols, we used the biomaRt (Durinck et al., 2009) R package.

Gene Set Enrichment Analysis
To investigate the significance of DE linear transcripts, we utilized the SetRank (Simillion et al., 2017) GSEA algorithm. This algorithm increases the confidence of a gene ontology (GO) analysis by removing gene sets that are flagged as significant if their significance is due to an overlap with other gene sets. SetRank analyzed the tables of top-ranked genes generated by edgeR against a background gene set filtered for all GO databases, including biological processes, molecular function, cellular compartment, and REACTOME.

Read Mapping and CircRNA Detection
The average percentage of uniquely mapped reads for each sample was 90%, illustrating the high quality of RNA template preparation and sequencing. CIRCexplorer2 detected significantly more circRNAs than DCC owing to DCC's high stringency and low false-positive rate. After merging the circRNAs to find those commonly detected by both tools, 1515 circRNAs were detected across all samples; the quantity of circRNA was expected, as circRNA is highly abundant in the brain (Gokool et al., 2020). All uniquely mapped read values and circRNAs detected by CIRCexplorer2 and DCC and those detected by both tools are provided in Supplementary Table 1.

t-SNE Analysis
To assess sample clustering, we employed the R package Rtsne to perform a t-SNE machine learning analysis on the detected circRNA expression levels for each sample ( Figure 1A). CircRNA expression was clustered by brain region rather than disease versus control, suggesting region-specific functionality of circRNAs. A t-SNE analysis of detected linear transcripts revealed less organized clustering of samples ( Figure 1B). Cortex samples were clustered closely together, and hippocampus samples were less clustered.

CircRNA Differential Expression in MTLE
Two separate differential gene expression comparisons were performed: control-cortex versus MTLE-cortex and controlhippocampus versus MTLE-hippocampus. Within the cortex samples, there was a single DE circRNA that is derived from the calmodulin regulated spectrin associated protein 1 (CAMSAP1) gene. Comparisons between MTLE and control hippocampal tissue revealed eight DE circRNAs. Table 1 summarizes the metrics of the nine DE circRNAs in both tissues.
Next, we investigated expression patterns of linear transcripts expressed by genes co-expressing DE circRNAs. Figure 2 shows the distribution of circRNA and linear transcript expression levels within individual biological replicates. Only two genes (Figure 2), Homer protein homolog 1 (HOMER1), and NMD3 ribosome export adaptor (NMD3), revealed significant DE of linear and circular products between MTLE and control hippocampus.

Estimation of miRNA Binding Capacity by DE CircRNAs
The ENCORI database provides AGO-CLIPseq data for circRNAs with miRNA-AGO complex binding sites. Because ENCORI uses the GRCh37 coordinates, each DE circRNA backspliced junction (BSJ) coordinate was converted to GRCh37 using the NCBI remap tool 1 . Using a custom Python script, we identified miRNA binding within the start and end

CircRNA-miRNA-mRNA Interaction Analysis
We explored potential associations between DE circRNAs, the miRNAs predicted to bind these circRNAs, and linear RNA targets that interact with these miRNAs. For example, circRNA-HOMER1 was DE in the hippocampus, so DE analysis of mRNAs, targeted by the 12 miRNAs predicted to bind to circRNA-HOMER1 (Table 2), was conducted on hippocampus samples. The edgeR-generated adjusted p-values were plotted to visualize transcript distribution as a volcano plot ( Figure 3A).

Gene Set Enrichment Analysis
For the transcripts targeted by bound miRNAs, a GO analysis was performed with SetRank, a GSEA algorithm that removes falsepositive hits. By utilizing detailed significance tests and graph theory, SetRank addresses the gene overlap and multiple testing problems associated with most GSEA algorithms and allows for querying of multiple different databases simultaneously (Simillion et al., 2017). For the circRNA-HOMER1 interaction network, the proteins were revealed to be located in the intracellular organelle, cytoplasm, and nucleus cellular components and had the molecular functions of heterocyclic compound binding and biological process response to stimulus. For circRNA-NMD3, the proteins were attributed to the nuclear lumen cellular compartment and had the molecular functions of heterocyclic compound binding and catalytic activity. Gene set interaction networks for HOMER1 and NMD3 generated by SetRank and visualized within Cytoscape (Su et al., 2014) are shown in Figure 4 and represent how each of the GSEA terms relate to each other. Nodes are color-coded by their corrected p-values, and the thickness and direction of the edges represent the size of the intersection and points from the least to the most significant. A full list of enriched terms ranked in order of importance for each DE circRNA-miRNA-mRNA network is included in Supplementary Table 3. Additionally, SetRankgenerated GSEA membership files reveal to which GSEA terms each protein is attributed, and these are included as Supplementary Files 1A-F.

DISCUSSION
In this study, we present, for the first time, changes in the circular transcriptome landscape within the temporal cortex and  Table 2), which were predicted to interact with circRNA-HOMER1.
(B) Volcano plot of the linear transcripts that are targeted by the miRNAs, which were predicted to interact with circRNA-NMD3. The X-and Y-axes of this plot show the log fold change and -log 10 adjusted p-value, respectively, indicating how differentially expressed a transcript is when compared to controls. Red indicates the transcripts with adjusted p-values of <0.05, which suggests that these linear transcripts were affected by this set of miRNAs interacting with circRNA. Labeled transcripts are the top 10 most DE.

FIGURE 4 | (A)
Gene set interaction network for the miRNA targets composing the HOMER1 circRNA-miRNA-mRNA interaction network. (B) Gene set interaction network for the miRNA targets composing the NMD3 circRNA-miRNA-mRNA interaction network. Node color indicates the corrected p-value, edge thickness relates to the size of intersection between two gene sets, and edge direction points from the least significant to the most significant.
hippocampus of MTLE patients. An integrated bioinformatics analysis revealed circRNA-miRNA-mRNA regulatory networks that are specifically formed in MTLE brains, thus showing that multiple compartments of the transcriptome are perturbed in this type of epilepsy. Among the thousands of neuronally expressed circRNAs, we identified one DE circRNA in the cortex and eight DE in the hippocampus of MTLE patients. Both linear and circular transcripts for two genes, NMD3 and HOMER1, were DE in the hippocampus. The concurrent expression pattern suggests that these circRNAs may be co-expressed with linear transcripts. Alternatively, these circRNAs may be regulating their own transcription through interactions with U1 small nuclear ribonucleoprotein (snRNP) or by regulating RNA polymerase II-mediated transcription (Kristensen et al., 2019). We compared our detected circRNAs to the list of circRNAs detected by Li et al. (2018) and found overlapping circRNAs, which showed a similar direction of expression. For circRNAs detected in the cortex, there was an overlap of 23 circRNAs and, in hippocampus, an overlap of 26 circRNAs. None of our DE circRNAs were detected. The lack of complete overlap may be explained by differences in bioinformatic analysis and the lack of tissue specificity in the aforementioned study.
Our bioinformatics pipeline established that seven of the nine DE circRNAs could be classified as miRNA sponges, i.e., they contained sequences complementary to miRNAs. Our interaction network analysis determined the effects of circRNA expression relative to miRNA binding to target mRNAs. These DE target mRNAs revealed functions in signal transduction and transcription for neurons of the hippocampus, implying a potential role in the etiology of MTLE-HS.
For circRNA-NMD3 and the set of miRNAs binding to circRNA-FBXL4, there was an overlap of miRNAs from the same family (Supplementary Table 2). This may reflect sequence conservation across miRNAs of the same family because miRNA families often share sequence, structure, and function (Kaczkowski et al., 2009).
CircRNA-NMD3 is expressed from the NMD3 ribosome export adaptor gene locus. NMD3 encodes a protein that transports the ribosomal 60S subunit to the cytoplasm (Trotta et al., 2003). In MTLE, the dentate gyrus is where excitatory mossy fiber axons lose contact with mossy cells and migrate back into the molecular layer as the densely packed granule cells migrate and disperse (Schmeiser et al., 2017). CircRNA-NMD3 is shown to be DE in the hippocampus of MTLE patients. CircRNA-NMD3 interacts with 15 miRNAs, and further investigation of those miRNAs identified DE of their target linear transcripts. Of the 27 DE transcripts, 11 were up-regulated and 16 were down-regulated in the hippocampus. GSEA of those proteins revealed involvement in molecular function, cellular component, biological processes, and REACTOME pathways. The first ranked term, heterocyclic compound binding, might be of significance as these compounds have been used in multiple antiepileptic drugs (Wei et al., 2015). These DE proteins in the hippocampus of MTLE individuals may inhibit the binding of these compounds and may explain the high rates of pharmacoresistance associated with the disease. REACTOME pathways for these proteins include signal transduction and metabolism of proteins, which may be affecting synaptic stimulation. As similarly revealed by Mills et al. (2020), the DE transcripts are shown to have a function in the immune system. CircRNA-FBXL4 is expressed from the gene locus of the F-Box and leucine-rich repeat protein 4. The FBXL4 protein colocalizes with mitochondria to maintain mtDNA (Bonnen et al., 2013). CircRNA-FBXL4 exhibits binding sites for 57 miRNAs of which seven comprise the miRNA family hsa-miR-378. There were 83 DE transcripts, and 33 were up-regulated and 50 down-regulated. GSEA revealed that the proteins involved in this network have biological processes of protein localization and regulation of the protein modification process. REACTOME terms included gene expression and metabolism of proteins. The cellular compartment location of these proteins is shown to be in the endoplasmic reticulum, which has a clear relation to these protein-specific functions. The unspecific functions of these proteins may be due to the enrichment for miRNA binding sites on circRNA-FBXL4 and the subsequently large pool of mRNA targets.
We observed significant differences in the expression levels of circRNA-HOMER1 and linear transcripts between MTLE hippocampus and controls. The HOMER1 gene encodes proteins involved in neuronal activity that bridge Group I metabotopic glutamate receptors (mGluR1/5) with inositol 1,4,5-triphosphate receptors (IP3Rs) on the endoplasmic reticulum, synaptic calcium ion channels, and NMDA receptor signaling complexes (Aloisi et al., 2017). After synaptic activity, a short splice variant HOMER1a is expressed as an immediate early gene and may antagonize group I mGluR activation (Tappe and Kuner, 2006). HOMER1a is implicated in epileptogenesis (Potschka et al., 2002;Celikel, 2007;Wagner et al., 2013). CircRNA-HOMER1 has binding sites for 12 miRNAs. There were 37 transcripts shown to be DE with 12 up-regulated and 25 down-regulated. GSEA of these proteins similarly revealed a molecular function in heterocyclic compound binding and ion binding. This may have the effect of modifying the hippocampal neurons' ability to bind heterocyclic neurotransmitters and antiepileptic medications and maintain the ion exchange required for generating action potentials.
Our study identified nine DE circRNAs between the cortex and hippocampus of individuals with MTLE. Of these, seven exhibited miRNA sponging characteristics, which likely have phenotypic roles in MTLE by contributing to, or resulting from, HS. Dissecting the effects of each miRNA to these interaction networks can advance our understanding of the individual roles played by miRNA sponging in MTLE pathogenesis and identify targets for further research. Future studies should investigate the effects of knocking down these target circRNAs in neuronal cell lines using in vivo circRNA-specific knockdown (Zimmerman et al., 2020) and assessing the subsequent changes in miRNA activity. If unbound miRNAs regulate mRNA translation, quantitative changes in protein expression should result. Additionally, more experimental validations are required to define the mechanisms of circRNA-miRNA binding and how the secondary RNA structure of a circRNA can facilitate these interactions. For circRNAs that do not interact with miRNAs, other mechanisms of action should be explored, such as RBP binding and protein scaffolding.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: European Genome-phenome Archive (EGA) which is hosted by the EBI and the CRG, under the accession number: EGAS00001003922.

AUTHOR CONTRIBUTIONS
MJ, EA, and OD conceived of the study. JM, SD, DF, MT, CS, and RT contributed to primary clinical and analytical design of the study. LG performed the data analysis. AC-H assisted with the data analysis. LG, JM, EA, OD, and MJ wrote the manuscript. AC-H, SD, and DF provided feedback to the manuscript. All authors read and approved the final manuscript.