Identification of epilepsy-associated neuronal subtypes and gene expression underlying epileptogenesis

Abstract

Epilepsy is one of the most common neurological disorders, yet its pathophysiology is poorly understood due to the high complexity of affected neuronal circuits. To identify dysfunctional neuronal subtypes underlying seizure activity in the human brain, we have performed single-nucleus transcriptomics analysis of >110,000 neuronal transcriptomes derived from temporal cortex samples of multiple temporal lobe epilepsy and non-epileptic subjects. We found that the largest transcriptomic changes occur in distinct neuronal subtypes from several families of principal neurons (L5-6_Fezf2 and L2-3_Cux2) and GABAergic interneurons (Sst and Pvalb), whereas other subtypes in the same families were less affected. Furthermore, the subtypes with the largest epilepsy-related transcriptomic changes may belong to the same circuit, since we observed coordinated transcriptomic shifts across these subtypes. Glutamate signaling exhibited one of the strongest dysregulations in epilepsy, highlighted by layer-wise transcriptional changes in multiple glutamate receptor genes and strong upregulation of genes coding for AMPA receptor auxiliary subunits. Overall, our data reveal a neuronal subtype-specific molecular phenotype of epilepsy.

Fig. 1. Multipatient single-nucleus transcriptomic dataset of the epileptic and nonepileptic temporal cortex.

a Schematic representation of the experimental outline for droplet-based single-nucleus RNA sequencing using 10× Chromium on FANS-isolated neuronal nuclei. Each sample was processed separately by FANS and 10× Chromium cDNA library preparation. b UMAP representation of neuronal nuclei isolated from multiple epileptic and nonepileptic cortices, and cell-type annotations for principal neurons and GABAergic interneurons. The colors represent subtypes with the labels showing subtype names. c General and family-specific marker expression for principal cells and GABAergic interneurons with the colors proportional to log-normalized expression values. d, e Family- and subtype-specific markers for principal cells (d) and GABAergic interneurons (e). Columns and rows represent subtypes and marker genes, respectively. The color shows the log2-fold change of this marker in a given subtype relative to the average expression in the other subtypes. f Confirmation of the layer-specific expression of cardinal markers for principal neurons in the healthy temporal cortex by in situ hybridization (taken from Allen Brain Atlas)^,. Scale bar: 400 μm.

Fig. 2. Integration of epileptic and nonepileptic datasets and identification of disease-related neuronal subtypes.

a UMAP embedding with the integration of epileptic and nonepileptic datasets using Conos, colored by condition. b The total number of nuclei identified per subtype and condition. c Percentage of nuclei per subtype, showing compositional change across conditions. A notable decrease in L2/3 subtypes as well as Pvalb_Sulf1 was observed for epilepsy. d Similarity score, based on gene-expression correlation between neuronal subtypes in the epileptic and nonepileptic cortex that reveals disease-related subtype-specific transcriptomic changes in the epileptic tissue. A lower similarity score indicates larger differences between conditions. The red line indicates the median similarity score across all subtypes. The green line represents the 0 level that corresponds to “no difference observed”. e, f Analysis showing overrepresentation of differentially expressed (DE) genes between epileptic and nonepileptic datasets in genes identified in genetic studies in human patients and mouse models (e), and epilepsy genes identified in the largest epilepsy GWAS study to date (f). The odds ratio of the Fisher’s test is shown on the y scale with the bar height corresponding to the conditional maximal likelihood estimate and whiskers showing 95% confidence intervals. The red horizontal line shows an odds ratio equal to 1, which corresponds to “no difference observed”.

Fig. 3. Identification of epilepsy-associated pathways and transcriptomic shifts across neuronal subtypes.

a GO-term enrichment analysis ordered by neuronal subtype reveals both subtypes with large transcriptomic changes (>100 GO terms) and subtypes with only few or no enriched GO terms in the epileptic dataset. The total number of GO terms that passed the 0.05 threshold for the adjusted P value of the overrepresentation test is shown on the y axis. Colors of the stacked barplot represent the top-level GO term: biological process (BP), cellular component (CC), or molecular function (MF). b The major groups of GO terms clustered by their level of enrichment per subtype reveal common transcriptomic shifts across neuronal subtypes in the epileptic brain. Rows correspond to GO terms, ordered according to hierarchical clustering. Columns correspond to cell types. Colors represent –log10 of adjusted P values of the overrepresentation test, trimmed with the upper boundary of 10. c For the blue cluster in (b), the plot shows a UMAP embedding of the GO terms per each subtype. Each point corresponds to a single square on the heatmap in (b). The distances between points are proportional to the Jaccard distance of the enriched genes between two given GO terms in certain subtypes. Thus, points, which are close to each other on the plot, are represented by similar sets of the enriched genes. Left—colored by subtype, right—colored by GO term. The numbers on the right panel indicate GO terms in a subcluster outlined by the dashed line. d Heatmap showing neuronal subtypes grouped based on Jaccard similarity of the enriched “Biological Pathway” GO terms. Rows and columns correspond to cell types, and the intersection represents the weighted Jaccard similarity between the two subtypes. Bold lines separate high-order clusters; neuronal subtypes labeled by the orange and green colors correspond to GABAergic interneurons and principal neurons, respectively. Such a clustering allows to identify groups of subtypes, where each group might correspond to a local circuit/network.

Fig. 4. Identification of signaling pathways and genes in cortical neuronal subtypes that might underlie seizure activity.

a The total number of enriched DE genes for GO terms “action potential”, “glutamate receptor signaling pathway”, and “AMPA glutamate receptor complex” in each neuronal subtype. b–d Expression level of DE genes found in the GO terms “action potential”, “glutamate receptor signaling pathway”, and “AMPA glutamate receptor complex”, respectively, ordered by neuronal subtype. The color of the points represents Z scores of differential expression between conditions, with the blue color showing downregulation and red colors indicating upregulation of gene expression in epilepsy. The size of the points corresponds to the average expression level of a gene in a given cluster. Points with low Z scores have higher transparency.

Fig. 5. Complex layer-wise dysregulation of genes involved in glutamate-mediated excitation in the cortex of epileptic patients.

a Notable layer-wise upregulation or downregulation of genes involved in glutamate-mediated excitation. b Overview of smFISH that shows CKAMP44 (SHISA9) and RORB mRNAs in the temporal cortex of control (nonepileptic) and epileptic brains. The RORB probe was used to label layer L4, thereby allowing the identification of the cortical layer structure (together with DAPI staining). Note the higher number of CKAMP44 mRNA molecules across all layers in the epileptic brain. c, d Representative image and quantification of CKAMP44 mRNA molecules in control (nonepileptic) and epileptic temporal cortices demonstrate upregulation of CKAMP44 in RORB-positive neurons belonging to L4_Rorb and L5–6_Fezf2_Lrrk1 principal neurons (one-tailed Mann–Whitney test, 98 cells in two sections of two control brains and 215 cells in three sections of three epileptic brains).

Fig. 6. Layer-wise dysregulation of gene expression of glutamate receptor subunits in the cortex of epileptic patients.

a Overview of a brain section showing DAPI labeling that was used for identification of the cortical layers. b–d Representative images of the temporal cortices of control (nonepileptic) and epileptic brains showing an increased number of mRNA molecules for GRIA1 in L2–3 (b), GRIA1 in L5–6 (c), and GRIN3A in L5–6 (d). e–g Quantifications showing the increase of mRNA molecules in the epileptic temporal cortex for GRIA1 in L2–3 (e: one-tailed Mann–Whitney test, 136 cells in three sections from three control brains and 113 cells in three sections from three epileptic brains), GRIA1 in L5–6 (f: one-tailed Mann–Whitney test, 145 cells in three sections from three control brains, and 117 cells in three sections from three epileptic brains) and GRIN3A in L5–6 (g: one-tailed Mann–Whitney test, 130 cells in two sections from two control brains and 143 cells in three sections from three epileptic brains).

Fig. 7. Dysregulation of genes underlying hypoinhibition in the cortex of epileptic patients.

a, b Visualization of GAD1 and GAD2 expression to assess hypoinhibition of GABAergic interneurons revealed a general decrease in expression of these genes involved in GABA synthesis. Total-count normalized expression level is shown on the y axis. Points represent median expression values, averaged over all samples for each subtype. Whiskers show 25 and 75% levels of expression, also averaged over all samples. c Decreased expression of cannabinoid receptor 1 (CNR1) in Vip and Id2 non-Lamp5 subtypes. d Representation image showing smFISH analysis of CNR1 transcript abundance in VIP-positive GABAergic interneurons in a nonepileptic cortex. e Quantitative analysis of CNR1 transcript abundance by smFISH demonstrating a decrease in CNR1 expression in VIP-positive GABAergic interneurons in epileptic cortex (one-tailed unpaired t test, 22 cells in three sections from three control brains and 43 cells in five sections from five epileptic brains).

Fig. 8. Neuronal subtypes most affected by epilepsy based on integrative analysis.

a Six metrics—expression similarity, cell-type composition, number of changed GO terms, enrichment in GWAS genes, enrichment in epilepsy genes, and number of DE genes—are aggregated into a single score. The six metrics are shown on the y axis, ordered by the metric weight, where the weight is represented by the blue vertical colorbar on the right and by the color transparency of the rows. Cell types grouped into interneurons/principal neurons are shown on the x axis, ordered by the total score (green horizontal colorbar on the top). The colors on the heatmap represent the strength of the effect: 0—not affected, 1—affected, 2—highly affected, 3—most affected. For a rational of how we assigned the level of effect and weight for each metric, see “Methods”. b Clustering of cell subtypes based on enrichment of certain GO terms (similar to Fig. 3d). Enlarged inlet shows co-clustering of cell subtypes with the most affected transcriptomes (based on the metric in Fig. 8a).