Patients
Four distinct epilepsy pathologies were considered in this study, namely TLE-HS, FCD IIa, FCD IIb and TSC cortical tubers. In addition, age- and tissue-matched control tissue samples were collected (control cortex n = 14; control hippocampus n = 13). Brain tissues included in this study were obtained from the archives of the Departments of Neuropathology of the Amsterdam UMC (Amsterdam, The Netherlands) and the UMC Utrecht (Utrecht, The Netherlands) (Supplementary Table 6). Cortical and hippocampal brain samples were obtained from patients undergoing surgery for intractable epilepsy and diagnosed with FCD type II (n = 17 FCD IIa, n = 33 FCD IIb), TSC cortical tubers (n = 21) and TLE-HS (n = 64), respectively (Table 2; more details in Supplementary Table 6).
All cases were reviewed independently by two neuropathologists (A.E. and A.M.). Patients who underwent implantation of strip and/or grid electrodes for chronic subdural invasive monitoring before resection and patients who underwent previous resective epilepsy surgery were excluded from this study. The classification of hippocampal sclerosis (HS) was based on analysis of microscopic examination as described by the International League Against Epilepsy6. The diagnosis of FCD was confirmed according to the international consensus classification system proposed for grading FCD9. All patients with cortical tubers fulfilled the diagnostic criteria for TSC cortical tubers (including genetic analysis for the detection of germline mutations)39. All FCD type II samples underwent deep sequencing using DNA extracted from snap-frozen surgical brain tissue targeting 13 genes (FCD panel SoVarGen, South Korea); analysis for replicated data was performed in accordance with a previous study40 (Supplementary Table 7).
Control material was obtained at autopsy from age- and brain area-matched control samples that were obtained at autopsy from individuals without a history of seizures or other neurological disease (Table 2; more details in Supplementary Table 6). Brain tissue was frozen and kept at − 80°C (for molecular analysis) or fixed in 4% paraformaldehyde and embedded in paraffin (FFPE) for histological analysis. All procedures received prior approval by the local ethics committee of the contributing medical centers, and were conducted in accordance with the guidelines for good laboratory practice of the European Commission.
RNA isolation
For RNA isolation, human tissue was homogenized in 700 µl Qiazol Lysis Reagent (Qiagen Benelux, Venlo, The Netherlands). Total RNA including the microRNA (miRNA) fraction was isolated using the miRNeasy Mini Kit (Qiagen Benelux, Venlo, The Netherlands) according to the manufacturer’s instructions. The concentration and purity of RNA was determined at 260/280 nm using a Nanodrop spectrophotometer (Ocean Optics, Dunedin, FL, USA) and RNA integrity was assessed using a Bioanalyser 2100 (Agilent Technologies, Santa Clara, CA, USA).
RNA-Seq library preparation and sequencing
All library preparation and sequencing were performed by GenomeScan (Leiden, The Netherlands). The NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA) was used to process the samples. Sample preparation was performed according to the protocol ‘NEBNext Ultra II Directional RNA Library prep Kit for Illumina’ (NEB #E7760S/L). Briefly, mRNA was isolated from total RNA using oligo-dT magnetic beads. After fragmentation of mRNA, cDNA synthesis was performed. Next, sequencing adapters were ligated to the cDNA fragments followed by polymerase chain reaction amplification. Clustering and DNA-sequencing was performed using the NovaSeq6000 (Illumina, Foster City, CA, USA) in accordance with manufacturers’ guidelines. All samples underwent paired-end sequencing of 150 nucleotides in length; the mean read depth per sample was 47 million reads.
The Decontamination Using Kmers (BBDuk) tool from the BBTools suite was used for adapter removal, quality trimming and removal of contaminant sequences (ribosomal or bacterial)41. A phred33 score of 20 was used to assess the quality of the read, with any read shorter than 31 nucleotides in length excluded from the downstream analysis.
Reads were aligned directly to the human GRCh38 reference transcriptome (Gencode version 33)42 using Salmon v0.11.343. Transcript counts were summarized to the gene level and scaled using library size and average transcript length using the R package tximport44. Genes detected in less than 20% of the samples in any diagnosis cohort and with counts less than six across all samples were filtered out, resulting in 28,366 genes for downstream analysis. The gene counts were then normalized using the weighted trimmed mean of M-values method with the R package edgeR45. The normalized counts were then log2 transformed using the voom function from the R package limma46.
Unsupervised hierarchical clustering and discriminant analysis on principal components
Unsupervised hierarchical clustering based on principal components was used to identify underlying structure in the gene expression matrix using the stats R package16. Next, a discriminant analysis of principal components (DAPC) was performed using optim.a.score to identify the optimal number of principal components to retain as implemented by the adegenet R package16,17.
Identification of gene coexpression modules
The details of the module identification workflow are described by Srivastava et al.10. Briefly, coexpression networks were constructed per epilepsy cohort using hierarchical clustering of normalized gene expression. First, as healthy matching control samples were age-matched across the general sample set, the age distribution was assessed per cohort before applying the workflow. In addition, any outliers due to area of resection or library preparation were removed. Next, only genes showing high variability across samples were retained (median absolute deviation [MAD] ≥ 0.25). For all remaining genes, the 1-Spearman rank correlation was computed for all gene pairs47–49 and used to construct the adjacency matrix (soft-thresholding power = 6)50. Unsupervised hierarchical clustering using Ward’s method identified the clusters of genes51 (from K = 1–200). The optimal number (Kx) was defined based on the second derivative of percentage of the variance explained (R²) per K52. Next, a leave-one-out bootstrapping procedure was implemented to assess the effect of samples on the stability and robustness of gene coregulation modules. For each permutation, gene coexpression modules were identified using the above-mentioned workflow and records of gene module membership. Cluster membership was used to construct the similarity matrix to identify genes assigned to the junk module based on an arbitrary threshold (50% assigned to junk module). The remaining genes were clustered based on the similarity matrix to obtain the coexpression modules. Finally, the modules were divided using (anti-)correlation of genes within the module. Based on the relative over- or underexpression of the module’s genes compared with healthy control samples, each submodule was assigned an ‘o’ or ‘u’ suffix, respectively.
To ensure the robustness of the identified modules, coexpression modules were only assembled in epilepsy cohorts with greater than 20 samples. An additional joint analysis was performed across all mTORopathies (FCD IIa, FCD IIb and TSC cortical tubers). The presence of outliers related to technical covariates was assessed using principal component analysis regression and removed from further analyses.
Differential coexpression
For each module the correlation between gene expression was calculated in both healthy controls and epilepsy patients to obtain the difference in median R². The empirical P value was estimated for each module by comparing the difference in median R² to the null distribution generated by performing 10,000 permutations of samples across cohorts10,53.
Phenotype association to module eigenGene
The relationship between module expression and the different reported phenotypes was explored using a linear model between each module’s eigenGene and the covariate: HS subtype, log10 of self-reported seizure frequency, sex, age, duration, sequencing group and library preparation batch. As duration also depends on the age of the patients, age was made an additional covariate when assessing association to duration.
Functional annotation using enrichment analysis
The modules were functionally annotated using multiple pathway resources (MetaBase, Reactome and GO) as well as cell type enrichment based on marker gene signatures derived from PanglaoDB54. A hypergeometric test was used to assess the significance of enriched pathway terms or marker gene signatures using a false discovery rate (FDR) correction to rectify for multiple testing using all expressed genes as a background55.
CRAFT framework: in silico causal reasoning
Candidate upstream regulators for the identified gene coexpression modules were predicted using the CRAFT framework. Srivastava et al.10 defined a causal reasoning framework that utilizes the direction of effects between the three components of the system, namely CMPs, TFs and target genes. The interactions between these three components and the direction of these interactions were obtained from the Clarivate Analytics MetaBase® (version 6.15.62452, https://clarivate.com/products/metacore/), a meta-database of manually curated literature-based contextual biological interactions. Details of the module identification workflow have been described by Srivastava et al.10. Briefly, all expressed membrane receptors, TFs and target genes from MetaBase were identified. Next, for each TF the set of target genes was retrieved as well as its activity (activation, inhibition, unspecified) and upstream membrane receptors affecting a TF and their effect were obtained using MetaBase® defined canonical linear pathways. The overall effect of the membrane receptor on the underlying module was defined by combining the separate effects of CMP-TF and TF-gene. The significance of effect of a regulator (TF or CMP) on a module was subsequently assessed by testing the overlap between genes under the control of the regulator and the genes belonging to a module (hypergeometric test), taking all expressed genes as the universe. FDR was calculated using Benjamini-Hochberg correction of enrichment P values, taking into account the total number of enrichment tests performed in testing55.
Identification of shared epilepsy regulations based on gene coexpression modules
The subsequent paragraph details the identification of specific epilepsy regulations as captured by gene coexpression modules in the independent epilepsy cohorts. Although different structural epilepsies are studied, similar pathways or mechanisms may still be dysregulated. To identify shared epilepsy regulations, the amount of gene content overlap between the gene coexpression modules from each epilepsy cohorts was identified using the inclusion index:
\(inclusion index = \frac{length\left(intersect\right(x,y\left)\right)}{min\left(length\right(x),length(y\left)\right)}\)
with x and y as two gene coexpression modules. Next, unsupervised hierarchical clustering based on Ward’s method was used to identify modules that showed overlap in gene content51 using the silhouette method to identify the optimal number of clusters. The analyses were performed with the stats and factoextra R packages56. By design, within an epilepsy cohort, a gene can only belong to one coexpression module. Therefore, the intersect between gene coexpression modules across epilepsy cohorts was defined as those genes occurring in at least one module per epilepsy cohort. This gene intersection was subsequently submitted to a hypergeometric test to obtain functional annotation with pathway resources (MetaBase, Reactome, GO) as well as cell type enrichment based on marker gene signatures derived from PanglaoDB54. Finally, the conservation of gene coexpression in other epilepsy cohorts and healthy control tissue was assessed with the same permutation approach as for differential coexpression analysis.
Immunohistochemistry
Human brain tissue was fixed in 10% buffered formalin and embedded in paraffin. Paraffin embedded tissue was sectioned at 6 µm, mounted on pre-coated glass slides (Star Frost, Waldemar Knittel Glasbearbeitungs, Braunschweig, Germany) and processed for immunohistochemical staining. Immunohistochemistry was carried out as previously described20 on samples from patients as reported in Supplementary Table 6. The following antibodies and dilutions were applied: SP-1 (SP-1, rabbit monoclonal, Abcam, ab124804, 1:200) and lysine-specific demethylase 1 (LSD-1) (LSD-1, rabbit polyclonal, Cell Signaling Technology, Cat#2139S, 1:200) incubated overnight at 4°C. For double labeling of SP-1 and LSD-1, sections were incubated with NeuN (mouse monoclonal, clone MAB377; Chemicon, Temecula, CA, USA; 1:2,000), glial fibrillary acidic protein (GFAP; mouse monoclonal, clone GA5, Sigma-Aldrich, St. Louis, MO, USA; 1:4,000) and HLA-DP/DR/DQ (HLA-II, mouse monoclonal, clone CR3/43, Agilent Technologies, Santa Clara, CA, USA; 1:100) antibodies, after incubation with the primary antibodies overnight at 4°C. For detection, sections were first incubated with Brightvision poly-alkaline phosphatase-anti-rabbit (DVPR55AP, Immunologic, Duiven, The Netherlands) for 30 min at room temperature, and washed with phosphate-buffered saline and then with Tris–HCl buffer (0.1 M, pH 8.2) to adjust the pH. Alkaline phosphatase activity was visualized with the alkaline phosphatase substrate kit III Vector Blue (SK-5300, Vector Laboratories Inc., CA, USA). After washing in phosphate-buffered saline, sections were secondly incubated with Brightvision poly-horseradish peroxidase-anti-mouse (DPVM55HRP, Immunologic, Duiven, The Netherlands) for 30 min at room temperature. Signal was detected using the chromogen 3-amino-9-ethylcarbazole (AEC, Sigma- Aldrich, St. Louis, MO, USA) in 0.05 M acetate buffer filtered substrate solution. Sections incubated without the primary antibodies or with the primary antibodies followed by heating treatment were essentially blank.
References
39. Northrup, H., et al. Updated international tuberous sclerosis complex diagnostic criteria and surveillance and management recommendations. Pediatr Neurol 123, 50–66 (2021).
40. Sim, N.S., et al. Precise detection of low-level somatic mutation in resected epilepsy brain tissue. Acta Neuropathol 138, 901–912 (2019).
41. Bushnell, B., Rood, J. & Singer, E. BBMerge - Accurate paired shotgun read merging via overlap. PLoS One 12, e0185056 (2017).
42. Harrow, J., et al. GENCODE: the reference human genome annotation for The ENCODE Project. Genome Res 22, 1760–1774 (2012).
43. Patro, R., Duggal, G., Love, M.I., Irizarry, R.A. & Kingsford, C. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods 14, 417–419 (2017).
44. Soneson, C., Love, M.I. & Robinson, M.D. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Res 4, 1521 (2015).
45. Robinson, M.D., McCarthy, D.J. & Smyth, G.K. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140 (2010).
46. Ritchie, M.E., et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43, e47 (2015).
47. Peterson, L.E. CLUSFAVOR 5.0: hierarchical cluster and principal-component analysis of microarray-based transcriptional profiles. Genome Biol 3, SOFTWARE0002 (2002).
48. van Houte, B.P. & Heringa, J. Accurate confidence aware clustering of array CGH tumor profiles. Bioinformatics 26, 6–14 (2010).
49. Otto, B., et al. Transcription factors link mouse WAP-T mammary tumors with human breast cancer. Int J Cancer 132, 1311–1322 (2013).
50. Zhang, B. & Horvath, S. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol 4, Article17 (2005).
51. Ward, J.H. Hierarchical grouping to optimize an objective function. J Am Stat Assoc 58, 236–244 (1963).
52. Tibshirani, R., Walther, G. & Hastie, T. Estimating the number of clusters in a data set via the gap statistic. J R Stat Soc Ser B Stat Methodol 63, 411–423 (2001).
53. Choi, Y. & Kendziorski, C. Statistical methods for gene set co-expression analysis. Bioinformatics 25, 2780–2786 (2009).
54. Franzén, O., Gan, L.M. & Björkegren, J.L.M. PanglaoDB: a web server for exploration of mouse and human single-cell RNA sequencing data. Database (Oxford) 2019, baz046 (2019).
55. Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol 57, 289–300 (1995).
56. Kassambara, A. & Mundt, F. Factoextra: extract and visualize the results of multivariate data analyses. R Package Version 1.0.7. (2020).