Article Sublayer- and cell-type-specific neurodegenerative transcriptional trajectories in hippocampal sclerosis Graphical abstract Authors Elena Cid, Angel Marquez-Galera, Manuel Valero, ..., Angel Barco, Jose P. Lopez-Atalaya, Liset M. de la Prida Correspondence jose.lopezatalaya@csic.es (J.P.L.-A.), lmprida@cajal.csic.es (L.M.d.l.P.) In brief Cid et al. use a combination of techniques from single-cell recordings in vivo to bulk tissue and single-nucleus RNA-seq to reveal cell-type-specific segregation of neurodegenerative and epileptogenic processes in experimental epilepsy. They find that segregated disease-specific gene expression changes run differentially across cell types and sublayers of the hippocampal CA1 region. Highlights d Superficial CA1 pyramidal neurons are overactive in epileptic rodents d Neurodegeneration-related gene signatures progress radially in the epileptic CA1 d Larger vulnerability of superficial CA1 pyramidal cells in experimental epilepsy d Cell-type-specific disease trajectories characterize hippocampal sclerosis Cid et al., 2021, Cell Reports 35, 109229 June 8, 2021 ª 2021 The Author(s). https://doi.org/10.1016/j.celrep.2021.109229 ll ll OPEN ACCESS Article Sublayer- and cell-type-specific neurodegenerative transcriptional trajectories in hippocampal sclerosis Elena Cid,1,8 Angel Marquez-Galera,2,8 Manuel Valero,1 Beatriz Gal,1,3 Daniel C. Medeiros,1 Carmen M. Navarron,2 Luis Ballesteros-Esteban,1 Rita Reig-Viader,4,5 Aixa V. Morales,1 Ivan Fernandez-Lamo,1 Daniel Gomez-Dominguez,1 Masaaki Sato,6 Yasunori Hayashi,6,7 Àlex Bayés,4,5 Angel Barco,2 Jose P. Lopez-Atalaya,2,* and Liset M. de la Prida1,9,* 1Instituto Cajal, CSIC, 28002 Madrid, Spain 2Instituto de Neurociencias, Universidad Miguel Hernández-Consejo Superior de Investigaciones Cientı́ficas (UMH-CSIC), 03550 Sant Joan d’Alacant, Alicante, Spain 3Universidad Europea de Madrid, 28670 Villaviciosa de Odón, Madrid, Spain 4Institut d’Investigació Biomèdica San Pau, 08041 Barcelona, Spain 5Universitat Autònoma de Barcelona, 08193 Bellaterra, Cerdanyola del Vallès, Spain 6RIKEN Brain Science Institute, Wako, 351-0198 Saitama, Japan 7Department of Pharmacology, Kyoto University Graduate School of Medicine, 606-8501 Kyoto, Japan 8These authors contributed equally 9Lead contact *Correspondence: jose.lopezatalaya@csic.es (J.P.L.-A.), lmprida@cajal.csic.es (L.M.d.l.P.) https://doi.org/10.1016/j.celrep.2021.109229 SUMMARY Hippocampal sclerosis, the major neuropathological hallmark of temporal lobe epilepsy, is characterized by different patterns of neuronal loss. The mechanisms of cell-type-specific vulnerability and their progression and histopathological classification remain controversial. Using single-cell electrophysiology in vivo and im- mediate-early gene expression, we reveal that superficial CA1 pyramidal neurons are overactive in epileptic rodents. Bulk tissue and single-nucleus expression profiling disclose sublayer-specific transcriptomic signa- tures and robust microglial pro-inflammatory responses. Transcripts regulating neuronal processes such as voltage channels, synaptic signaling, and cell adhesion are deregulated differently by epilepsy across sub- layers, whereas neurodegenerative signatures primarily involve superficial cells. Pseudotime analysis of gene expression in single nuclei and in situ validation reveal separated trajectories from health to epilepsy across cell types and identify a subset of superficial cells undergoing a later stage in neurodegeneration. Our findings indicate that sublayer- and cell-type-specific changes associated with selective CA1 neuronal damage contribute to progression of hippocampal sclerosis. INTRODUCTION The most common form of hippocampal sclerosis (type 1; 60%–80% of TLE cases) shows severe neuronal loss of CA1, Epilepsies are brain disorders characterized by an enduring pre- CA3, and CA4 pyramidal neurons and milder loss in CA2, disposition to generate seizures with emotional and cognitive with variability along the anteroposterior axis (Blu €mcke et al., associated comorbidities. Despite significant therapeutic ad- 2013; Wyler et al., 1992). Other cell types, including microglia vances, one-third of affected individuals remain resistant to and astrocytes, are also affected (Blu € mcke et al., 2013; pharmacotherapy (Chen et al., 2018). Temporal lobe epilepsy Morin-Brureau et al., 2018). In contrast, type 2 hippocampal (TLE), the most prevalent form of pharmacoresistant epilepsy, sclerosis (10%–20% of cases) is associated with predominant is frequently associated with hippocampal sclerosis (Blu €mcke CA1 neurodegeneration and minimal loss in other regions et al., 2013). Hippocampal sclerosis is characterized by specific (Blu€mcke et al., 2013). Given the disparities between clinical se- patterns of neuronal loss affecting different hippocampal sub- ries, there is no consensus on whether neuronal loss pro- fields from the CA1 to CA3/4 areas, the hilus of the dentate gyrus, gresses over time (Blu € mcke et al., 2009; Davies et al., 1996). and superficial layers of the entorhinal cortex (Blu €mcke et al., In addition, individual variabilities and anatomical inhomogene- 2009; Du et al., 1995; de Lanerolle et al., 2003). Factors such ities complicate classification (Coras et al., 2014; Rodrigues as epilepsy history, age of onset, and the relationship with early et al., 2015; Saghafi et al., 2018). For instance, patchy neuronal precipitating events may all influence the degree and severity of loss has been described in the CA1 region in some cases, hippocampal sclerosis (Davies et al., 1996). whereas in others, it seems to adopt a more laminar profile Cell Reports 35, 109229, June 8, 2021 ª 2021 The Author(s). 1 This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). ll OPEN ACCESS Article Figure 1. Differential responses of CA1 pyramidal cells during epileptiform activities (A) Intracellular and multi-site local field potential (LFP) recordings from epileptic rats. Cells were identified with streptavidin and tested against Calbindin. Scale bar, 100 mm. SO, stratum oriens; SP, stratum pyramidale; SR; stratum radiatum. (B) Intracellular activity during IID at different membrane potentials (gray traces). Traces are aligned by the peak of IID recorded at the SR. HFOs were recorded at the SP. (C) Responses of the cell shown in (A) during SPW fast ripples. (D) Deviation from the resting membrane potential (RMP) recorded in individual cells during IIDs (n = 12 cells). Red colors reflect depolarization and blue hy- perpolarization. Cells are ranked by their distance to the SR and classified as deep and superficial (Sup) (subplot at right). The discontinuous line marks sublayer limits. (E) Same as in D for SPW-fast ripples (n = 19 cells). (F) Mean membrane potential responses around IID events showed differences between deep and Sup cells (Friedman c2 [1,333] = 28.7, p < 0.001). Data are from 6 deep and 6 Sup CA1 pyramidal cells. (G) Same as in (G) for SPW fast ripple events. Note the larger after-event depolarization in Sup cells (Friedman c2 [1,333] = 14.67, p < 0.0001). Data are from 8 deep and 11 Sup CA1 pyramidal cells. (H) cFos immunoreactivity in representative CA1 sections from one rat exhibiting bilateral forelimb clonus after sound stimulation (seizure) versus a non-stim- ulated epileptic rat with no observed seizure (basal). Scale bar, 80 mm. (I) Intensity of c-Fos from all pyramidal cells in one confocal section per rat (small dots) and mean data per animal (larger dots). Significant interaction was confirmed for animals between groups and sublayers (p = 0.0044; n = 4 epileptic basal, n = 3 epileptic seizures). Post hoc differences: *p < 0.05. (J) Schematic of electrophysiological and histopathological findings. (Prada Jardim et al., 2018). In some affected individuals, cell tum) and deep (closer to oriens) CA1 pyramidal neurons loss concentrates in the CA4 region and the dentate gyrus project differentially and diverge in their participation of and is frequently integrated in dual pathologies; e.g., TLE and sharp-wave ripple activity, theta-gamma oscillations, and malformations of cortical development, classified as type 3 hip- behavioral-cognitive correlates (Soltesz and Losonczy, 2018; pocampal sclerosis (3%–7%) (Mathern et al., 1997). The mech- Valero and de la Prida, 2018). In spite of data suggesting critical anisms underlying specific vulnerability of diverse cells and regionalization of CA1 neuronal responses to ischemia, anoxia, their role in the histopathological landscape and clinical signif- and epilepsy (Morris et al., 1995; Valero et al., 2017; Wittner icance remain unknown. et al., 2002), little is known about their clinical relevance and Recent techniques and methods operating at single-cell res- potential relationship with neuronal vulnerability. Understanding olution point to an exquisite cell-type-specific organization that the effect of cellular diversity and transcriptional changes in ep- is instrumental for brain function (Habib et al., 2016; Zeisel ilepsy progression may lend insights into more specific mech- et al., 2015). In the hippocampus, the CA1 region is organized anisms toward new diagnostic and therapeutic opportunities radially in two distinct sublayers with characteristic gene (Pfisterer et al., 2020). expression gradients along the anteroposterior and proximo- Here we combine gene expression profiling at the single-nu- distal axes (Cembrowski et al., 2016a; Dong et al., 2009; Slo- cleus and microdissected tissue levels with single-cell electro- mianka et al., 2011). Functionally, superficial (closer to radia- physiology to disclose epileptogenic and neurodegenerative 2 Cell Reports 35, 109229, June 8, 2021 ll Article OPEN ACCESS changes running differentially across CA1 sublayers in an exper- Sublayer regionalization of epilepsy-associated imental model of hippocampal sclerosis. Our study highlights the transcriptional responses importance of leveraging cell type specificity to better under- To investigate how the differential responsiveness of deep and stand the phenotypic complexities accompanying hippocampal superficial cells may relate with distinct transcriptional re- sclerosis in epilepsy. sponses across sublayers, we performed RNA sequencing (RNA-seq) analysis of laser capture microdissection (LCM) sam- RESULTS ples from control (saline-injected) and epileptic rats (n = 3 repli- cates each; Figures 2A, S2A, and S2B; Table S2). Large activity burden in superficial CA1 pyramidal cells Analysis of LCM RNA-seq data revealed sublayer-specific during epileptiform activities genes common to control and epileptic animals (absolute log2 Single CA1 pyramidal cells were recorded intracellularly from fold change [FC] > 0.5; adj. p < 0.01; Figure 2B; superficial versus anesthetized epileptic rats to evaluate their intrinsic excitability deep samples). Among these, we retrieved bona fide marker and activity during interictal discharges (IIDs), sharp-wave genes of superficial (Calb1, Grm1, and Syt17) and deep CA1 py- (SPW) fast ripples, and ictal discharges (Figures 1A–1C; Figures ramidal cells (Ndst4), consistent with previous data for mice S1A and S1B). Recorded cells were identified post hoc (with (Cembrowski et al., 2016b, 2016a) (Figure S2C; see Figures streptavidin) and immunostained against Calbindin (Calb) to S2D and S2E for validation by in situ hybridization). Sublayer classify them as deep (negative) or superficial pyramidal cells gene expression analysis confirmed preserved regionalization (positive) (Figure 1A). SPW-associated high-frequency oscilla- in epileptic rats with only a subset of uncorrelated transcripts tion (HFO) events were detected automatically and classified (FC > 0.5; adj. p < 0.01; Figure 2C, red dots; see Figure S2F for as ripples (100–150 Hz), fast ripples (>150 Hz) and IID using less stringent criteria: adj. p < 0.1). Strikingly, we found that amplitude and spectral information (Figures S1C and S1D), as some gene markers of other cell types in the hippocampus, before (Valero et al., 2017). including interneurons (Sst, Vip, and Kit), oligodendrocytes Intracellular activities recorded during IID events (n = 12 cells; (Mbp, Plp1, and Mobp), astrocytes (Aqp4, Gja1, and Gfap), Figure 1B) and SPW fast ripples (n = 19 cells; Figure 1C) were and homeostatic microglia (Tgfbr1 and Csf1r), also exhibited typically associated with consistent depolarization from the sublayer differences, consistent with heterogeneous distribution resting membrane potential (RMP), and firing of all pyramidal of cell types across the CA1 radial axis (Figures 2B, 2C, and S2F) cells was examined (Figures 1D and 1E). A temporal analysis (see below). of membrane potential changes (30-ms bins) showed differ- Next, we investigated the transcriptional changes associated ences between cell types during both IID events (Figure 1F; with epilepsy across CA1 sublayers (epileptic versus control Friedman c2 [1,333] = 28.7, p < 0.001) and SPW fast ripples (Fig- samples per sublayer). We identified 103 significantly differen- ure 1G; Friedman c2 [1,333] = 14.67, p < 0.0001), with superficial tially expressed genes (DEGs) in the deep sublayer of epileptic cells consistently showing larger depolarization at the event versus control rats versus 376 DEGs in the superficial sublayer peak. Post-event membrane potential responses recorded (FC > 0.5; adj. p < 0.01; Figure 2D). Although epilepsy-associated showed effect for the type of events (F(2,57) = 7.17, p = transcriptional changes run in similar directions in both sublayers 0.0018), sublayer (F(1,57) = 5.44, p = 0.0236) and interaction (i.e., no counter-regulated genes), they were more exacerbated (F(2,57) = 4.89, p = 0.0113). Post hoc Tukey tests confirmed sig- in the superficial CA1 (Figures 2D and S2H). Overall, we retrieved nificant smaller hyperpolarization (p = 0.0020) and a higher firing more upregulated than downregulated genes, particularly in the rate (p = 0.045) following IID events in superficial compared with superficial (73% DEGs) versus the deep sublayer (57% DEGs) deep cells. We found some minor differences of intrinsic excit- (Figure 2D, Venn diagrams on the right). Because of the small ability between cell types and groups (Table S1). size of the gene set differentially regulated by epilepsy in the We reasoned that superficial pyramidal cells should be more deep sublayer, a less stringent threshold was used to perform steadily activated during seizures, given their poor post-event functional enrichment analysis with gene sets of equal sizes. hyperpolarization. We exposed 3 epileptic rats to high-pitched We first required the genes to show a modest level of statistical sounds to promote convulsive seizures (pulses of 95–100 dB significance (adj. p < 0.1) and ranked them by FC to retain the top and 1- to 20-s duration at random intervals of 1–20 s during 250 upregulated genes and top 250 downregulated genes. 10 min) (seizure group). This group was compared with 4 unsti- Functional enrichment analyses of these gene sets of equal mulated epileptic rats without any observed clinical seizure size also revealed marked differences (Figure 2E). (basal group). Animals were sacrificed after 1 h to evaluate Downregulated DEGs in epilepsy in the superficial CA1 sub- expression of the immediate-early gene cFos (Figure 1H). Statis- layer show robust functional association with many Gene tics per cell confirmed effects for interaction between groups Ontology (GO) terms (adj. p < 0.05) linked to neuronal pro- and sublayer (p < 0.0001; 2-way ANOVA). Significant interaction cesses, such as synaptic transmission and neuron projection, was confirmed for animals (p = 0.0044) with post hoc differences including many well-known modulators of epileptogenic pro- between groups for superficial cells only (p = 0.0395) and be- cess, such as Cacnb1, Cacnb4, Cacng2, Kcnq3, Kcnab1, tween deep and superficial cells within the seizure group (p = Scn4b, Grin2a, Nptn, and Gria4 (Winden et al., 2015) (Figure 2E). 0.0117) (Figure 1I). Most of these Gene Ontology (GO) terms were significantly en- These results suggest higher responsiveness of superficial riched only in the superficial CA1 sublayer (Figure 2E). CA1 pyramidal cells during epileptiform activities in chronically Regarding upregulated genes in epilepsy, there were no GO epileptic rats (Figure 1J). terms significantly enriched using the gene set from the deep Cell Reports 35, 109229, June 8, 2021 3 ll OPEN ACCESS Article Figure 2. Bulk tissue gene expression profiling of the epileptic hippocampal area CA1 reveals regionalized transcriptional response (A) Representative laser capture microdissection (LCM) sampling of the Sup CA1 sublayer. Scale bar, 100 mm. (B) Heatmaps of DEGs in deep (left) and Sup CA1 sublayers (right) in three replicates from control (C) and epileptic (E) rats (FC > 0.5, adj. p < 0.01). Arrowheads point to bona fide gene markers of Sup (Calb1, Grm1, and Syt17) and deep (Ndst4 and Nr4a2) pyramidal neurons. Note the presence of gene markers for in- terneurons (e.g., Vip, Sst, Sema3c, and Kit) and oligodendrocytes (Plp1, Mbp, Mobp, Mal, Enpp2, Cldn11, and Ermn). (C) Scatterplot of DEGs between Sup and deep CA1 sublayers in control or epileptic rats (FC > 0.5, adj. p < 0.01). Bona fide gene markers of Sup (Calb1, Grm1, and Syt17) and deep (Ndst4, Nr4a2, and Col11a1) pyramidal neurons are highlighted (black text and bold font). Also shown are marker genes for interneurons (Sst, Vip, and Kit), oligodendrocytes (Mbp, Mobp, and Plp1), astrocytes (Aqp4 and Gja1), and microglia (Tgfbr1). Note the presence of a subset of uncorrelated transcripts, including canonical markers of microglia cells (Trem2, Irf8, Fcgr2b, and Cd68), at the Sup sublayer in epileptic rats (red). MLE, maximum-likelihood estimate. (D) MA plots showing epilepsy-associated significantly upregulated (red) and downregulated (blue) genes in deep and Sup CA1 sublayers (FC > 0.5, adj. p < 0.01) (left). Venn diagrams show the overlap of upregulated (top) and downregulated (bottom) genes between epileptic and control in Sup and deep sublayers of CA1. (E) GO analysis of the top 250 upregulated and downregulated genes in deep and Sup sublayers in epilepsy (adj. p < 0.1 and then ranked by FC). sublayer (adj. p < 0.05). Strikingly, significantly upregulated tions within the CA1 sample and/or changes in cell type compo- transcripts in the superficial sublayer of epileptic CA1 show a sition in response to epilepsy. strong functional association with immune and inflammatory re- To address this point, we reanalyzed a published dataset of sin- sponses, including cytokine production, cell migration, cell pro- gle-cell RNA-seq (scRNA-seq) from 3,005 barcoded individual liferation, and programmed cell death (e.g., Tlr3, Tlr7, C1qa, cells from the mouse somatosensory S1 cortex and CA1 regions Csf1r, Tgfb1, Cd74, Cd86, Itgax, Itgam, Itgb2, Casp1, Casp4, (Zeisel et al., 2015; Figure S2I). We first performed unsupervised Nlrc4, Pycard, and Tnfrsf1b) (Figure 2E). Upregulated genes in consensus clustering to identify major cell types in CA1 and ca- the superficial CA1 sublayer of epileptic animals include core nonical marker genes (Figures S2J and S2K). Next, to identify components of the inflammasome, such as Nlrc4, Casp1, and cell-type-specific gene signatures in bulk tissue RNA-seq data, Casp4, which, upon activation, can promote pyroptosis, a we used scRNA-seq expression levels for the top 250 DEGs be- pro-inflammatory form of lytic cell death (McKenzie et al., tween CA1 sublayers of control and epileptic rats to summarize 2020). The transcript levels of the cysteine protease Casp3, a the cells by linear dimensionality reduction with principal-compo- central effector of apoptosis, was also found to be upregulated nent analysis (Figure 3A). In addition, we also performed pairwise specifically in the epileptic CA1 (global effect, adj. p < 0.05; su- correlations and hierarchical clustering for these gene sets (Fig- perficial, adj. p = 0.0612; deep, adj. p = 0.6443). Thus, differen- ure 3A). In control animals, these analyses revealed that gene tial gene expression analysis suggests that the nature and signatures of distinct cell types, including pyramidal neurons, as- severity of transcriptional responses to TLE may segregate trocytes, and oligodendrocytes, were prominent in the CA1 deep radially across CA1. layer, whereas genes enriched in the superficial sublayer were mostly associated with pyramidal neurons and interneurons. In Cell type deconvolution of LCM RNA-seq data suggests sharp contrast, we noticed strong upregulation of a core gene sublayer-specific neurodegeneration associated with signature that included many known microglial markers (Micro; microglia e.g., Tyrobp, Trem2, Irf8, Irf5, Fcgr2b, Cd74, Cd68, Itgam, and Given the sublayer heterogeneity of biological processes re- Arpc1b) and was specific to the superficial layer of the epileptic vealed by bulk-tissue LCM RNA-seq, we reasoned that they CA1 (Figures 3A and 3C). We further confirmed these findings us- might reflect different contributions by discrete cellular popula- ing a different scRNA-seq dataset from the Allen Brain Map portal 4 Cell Reports 35, 109229, June 8, 2021 ll Article OPEN ACCESS Figure 3. Microglia subpopulations underlie the transcriptional signature of Sup CA1 in epilepsy (A and B) Subpopulation signatures inferred by deconvolution of bulk tissue transcriptome profiles. Gene sets were the top 250 DEGs between Sup and deep CA1 sublayers in epileptic and control rats (as indicated) identified in bulk tissue RNA-seq. For the selected genes, normalized expression was retrieved from publicly available scRNA-seq data from the mouse CA1 hippocampal region (Zeisel et al., 2015) (A) or the Allen Brain Map portal (Mouse Whole Cortex and Hippocampus SMART-seq [2019] with 10x-SMART-seq taxonomy [2020]) (B), and single cells were summarized by linear dimensionality reduction using principal-component (legend continued on next page) Cell Reports 35, 109229, June 8, 2021 5 ll OPEN ACCESS Article (Mouse Whole Cortex and Hippocampus SMART-seq [2019]) with trated at the CA1 region in the systemic kainate model of TLE in 10x-SMART-seq taxonomy (2020) (Figure 3B). To evaluate this in rats, more consistent with sclerosis type 2 (Inostroza et al., situ, we performed immunofluorescence staining against protein 2011). To precisely evaluate CA1 sublayers, we combined immu- markers for microglia (Iba1), astrocytes (glial fibrillary acidic pro- nostaining against the CA1 marker Wfs1 with the superficial tein [GFAP]), and oligodendrocytes (Olig2) and quantified their marker Calbindin (Figure 4A). Proximal, intermediate, and distal density across sublayers (Figure 3D, left). We confirmed the pres- CA1 segments were evaluated separately in tissue from 10 control ence of physiological segregation of glial cell types across the and 10 epileptic rats (1 section per animal from 3.2 to 4.8 mm CA1 sublayers (control) and changes in TLE, with between-groups from the bregma). To avoid effects of hippocampal atrophy over differences reaching significance for microglia (Figure 3D, right). cell density per area, we chose to evaluate the linear density of cells Prompted by these results, we dissected the contribution of mi- along the stratum pyramidale (SP) (cells per millimeter). croglia-associated transcripts in our sample. We noticed that most Cell loss in all three segments was mostly restricted to Calb+ su- of the uncorrelated transcripts in the comparison of sublayer-en- perficial CA1 pyramidal cells in chronically epileptic rats (Figure 4B; riched genes across conditions were highly expressed in microglia 3-way ANOVA for group, sublayer, and proximodistal effects, all (red dots in Figures 2C and S2F; see microglia-specific genes in significant at p < 0.0001; post hoc Tukey test, p < 0.001). Because Figure S2G). Thus, we evaluated their potential functional effect Calbindin immunoreactivity may be affected during epileptogene- by building the microglial sensome (Hickman et al., 2013), which sis (You et al., 2017), we exploited the transgenic mouse line was upregulated profoundly in the superficial sublayer of epileptic Thy1.2-G-CaMP7-DsRed2 with restricted expression on deep py- rats (Figure S3A). Moreover, evaluation of microglia-neuron inter- ramidal cells (Figure 4C). Using this line, we estimated superficial actions via ligand-receptor pairing confirmed sublayer-specific ef- cells as Wfs1+/DsRed2 and deep cells as Wfs1+/DsRed2+. Esti- fects (Figure S3B). For example, in the epileptic superficial CA1 mated global error by off-target Calb+/DsRed2+ expression was subfield, we found dysregulated expression of the transcript en- at 12% (10% proximal, 10% intermediate, 16% distal). We found coding the CD200 receptor, whose expression in neurons pro- similar regionalization of neuronal loss at the intermediate and vides a ‘‘don’t eat me’’ signal that reduces microglial activation distal segments of chronic epileptic mice using co-localization be- (Hoek et al., 2000). Similarly, Cx3cl1, a neuronal chemokine that tween DsRed2 and Wfs1 (Figure 4D, left; n = 3 control, n = 5 dampens the microglia inflammatory response and neurotoxicity, epileptic; 6–8 weeks post-status; significant interaction between was also strongly downregulated, whereas C3 and Csf1 tran- group and sublayer, p = 0.0006; post hoc t test, p < 0.01). scripts were upregulated in the superficial CA1 sublayer (Biber Given a potential interaction between expression of Calbindin et al., 2007; Luo et al., 2013). In contrast, the epileptic deep CA1 and aging (Iacopino and Christakos, 1990), we evaluated tempo- sublayer showed a much weaker microglia signature, and only ral trends of CA1 cell loss in epileptic rats, where a time-depen- Cx3cl1 was found to be downregulated (Figure S3B). dent decrease of Calb+ cells was appreciated (r = 0.52, p = Our analyses reveal heterogeneous distribution of resident cell 0.097, Pearson correlation), but not for control animals (p = types across the CA1 radial axis of the normal hippocampus. 0.66; Figure 4E). We also examined the sublayer-specific vulner- Strikingly, we identified critical changes in the epileptic hippo- ability of CA1 neurons in the lithium-pilocarpine model in rats campus, including local accumulation of reactive microglia and (Curia et al., 2008), which is more consistent with type 1 hippo- changes of key modulators of the neuronal immune response campal sclerosis (Inostroza et al., 2011). Notably, pilocarpine- specifically at the superficial sublayer. treated rats showed a similar reduction of Calb+ CA1 pyramidal cells early along epileptogenesis (p = 0.0003, n = 5 rats, 2– CA1 hippocampal sclerosis is degenerative and 13 weeks post-status; Figure 4E, note full cell loss after 10 weeks sublayer specific post-status, arrowhead). Finally, to fully exclude an interaction Given the results described above, we sought to evaluate the with Calbindin immunoreactivity, we evaluated Wfs1+ cell den- anatomical landscape of CA1 hippocampal sclerosis. Previously, sity at deep and superficial sublayers by relying only on anatom- we reported a significant reduction of linear cell density concen- ical criteria (location within the sublayer) and found consistent analysis (PCA) (top panels). Cells are colored by population membership (STAR Methods). Also shown are heatmaps of pairwise correlation for selected gene sets using scRNA-seq data, as indicated above (bottom panels). Bona fide markers of distinct cell types were detected in clusters of highly correlated genes rep- resenting cell type gene signatures convoluted in the bulk tissue RNA-seq (Pyr, pyramidal cell; Inter, interneuron; ODC, oligodendrocytes; Astro, astrocytes; Endo, endothelial cells; Micro, microglia; Mural, mural cells). Gene signatures of different types of cells (e.g., Pyrs, Astros, and ODCs) in the deep CA1 sublayer of control and epileptic rats lead to segregation of the cells in the corresponding PCAs and heatmap of the pairwise correlation matrix. Note the radial distribution of distinct cell types in the CA1 region in control and epileptic rats. Also note the presence of a strong gene signature of Micros in the Sup CA1 sublayer of epileptic rats (arrowhead). Names of highly correlated genes enriched in Micros are shown (Zeisel et al., 2018). (C) Heatmaps of significantly differentially expressed cell marker genes (Figures S2I–S2K) in epilepsy in the deep or Sup sublayer (left panels), heatmap showing the sublayer effect in control and epileptic groups (paired design) for cell-marker genes modulated in epilepsy (center panels), and cell types associated with each marker gene (right panel). Pyr, pyramidal cells. *adj. p < 0.1. (D) Representative sections of control and epileptic rats immunostained against GFAP (Astros), Olig2 (ODCs), and Iba1 (Micros) (left panels). The leftmost section for each condition shows co-localization between the cell type marker (cyan), Calb (purple), and bisbenzimide (Bis; white). Bar plots show quantification of astrocyte (GFAP+), ODC (Olig2+), and Micro (Iba1+) linear cell density across sublayers and groups. Shown are significant differences in interaction between sublayer and groups for Astros (F(1,5) = 10.1, p = 0.022) and Micros (F(1,5) = 7.4, p = 0.042) with post hoc differences at deep and Sup layers in the control group, the effect of the sublayer for ODCs (F(1,5) = 23.8, p = 0.0045), and differences between groups only for Micro cells. Data are from 3 control and 4 epileptic rats. *p < 0.05, **p < 0.01, post hoc tests. Scale bar, 60 mm. 6 Cell Reports 35, 109229, June 8, 2021 ll Article OPEN ACCESS Figure 4. Hippocampal sclerosis is sublayer and cell type specific (A) Immunostaining against the CA1 marker Wfs1 co-localized with Calb at proximal, intermediate, and distal segments in control and epileptic rats. Scale bar, 50 mm. (B) CA1 linear cell density confirmed neuronal loss affecting mainly Calb+ cells. Data are from 10 control and 10 kainate-treated epileptic rats (1 section per rat, 3.2 to 4.8 mm from the bregma, 203). Significant 3-way ANOVA for group (F1,120) = 32.6, p < 0.0001; proximodistal (F(2,120) = 11.7, p < 0.0001) and sublayer (F(1,120) = 424,1; p < 0.0001). Post hoc Tukey test: **p < 0.01, ***p < 0.005. Scale bar, 50 mm. (C) Hippocampal sections from representative control and epileptic transgenic (TG) mice expressing G-GaMP7-DsRed2 in Calb deep Pyrs. (D) Linear density data from 3 control and 5 epileptic TG mice (1 section per mouse at 203). Significant 3-way ANOVA for groups (F(1,36) = 7.8, p = 0.008) and sublayers (F(1,36) = 8.9, p = 0.0051), no proximodistal effect. Post hoc t test: **p < 0.01, ***p < 0.005. (E) Left: experimental timeline. Right: mean density of Calb+ and Calb CA1 Pyrs as a function of time fitted by linear trends. Data are from 10 control, 10 epileptic (kainate), and 5 epileptic (lithium-pilocarpine) rats. Arrowheads point to a lithium-pilocarpine rat with full cell loss in CA1 10 weeks after status. (F) CA1 mean cell loss in rats quantified as deep versus Sup using Wfs1. Sup CA1 cells are more affected than deep cells. Significant effect for group in a 2-way ANOVA, F(2,48) = 5.9, p = 0.005. Post hoc unpaired t test: *p < 0.01, **p < 0.01. (G) Fluoro-Jade signal co-localized with Wfs1 and quantification of Fluoro-Jade+ cells across sublayers; n = 8 epileptic rats. Paired t test, *p < 0.01. Scale bar, 40 mm. (H) Density of deep (Calb, left) and Sup (Calb+, right) Pyrs against the density of Micros (Iba1+) counted at deep and Sup sublayers (data are from 6 epileptic rats). Significant Spearman decreasing monotonic trend only for Sup cells (r = 0.87, p = 0.0236). (I) Schematic of epilepsy-associated histopathological findings indicating larger neuronal loss (blue ramp) and increased Micro cells (red) in the Sup sublayer. results in both TLE models (Figure 4F; 2-way ANOVA for group, superficial pyramidal neurons is strongly associated with local p = 0.005). accumulation of microglia. To further confirm sublayer-specific neurodegeneration, we Our results support the idea that hippocampal sclerosis re- combined Fluoro-Jade staining, which characteristically labels sults from specific interactions between different cell types in a degenerating cells, with Wfs1 immunostaining in sections from particular niche (i.e., superficial CA1 pyramidal cells and acti- epileptic rats (n = 8 epileptic rats) and again observed stronger vated microglia), leading to regionalized neurodegenerative sig- neurodegeneration in superficial sublayers (Figure 4G; paired t natures in the sclerotic CA1 (Figure 4I). test, p = 0.0017). Notably, double immunostaining against Calb and the microglial marker Iba1 revealed an inverse correlation Single-nucleus RNA-seq confirms cell-type-specific of the amount of microglia and the density of superficial pyrami- neurodegeneration of CA1 pyramidal cells dal neurons, consistent with a decreasing monotonic trend To gain more insights into the transcriptional activity patterns un- (Calb+; r = 0.87, p = 0.0236; Spearman correlation but not for derlying hippocampal sclerosis, we performed unbiased high- deep cells [Calb]; Figure 4H). This result indicates that loss of throughput RNA-seq of isolated single-nucleus RNA-seq Cell Reports 35, 109229, June 8, 2021 7 ll OPEN ACCESS Article (snRNA-seq) (Figure 5A). We analyzed the transcriptomes of related to neuronal survival. For instance, we found significant 6,739 single high-quality CA1 barcoded nuclei derived from 2 upregulation of death-associated protein kinase 1 (Dapk1) control and 2 epileptic mice (control, 3,661 nuclei; epilepsy, specifically in superficial CA1 neurons (Bialik and Kimchi, 3,078 nuclei). Consistent with recent reports (Habib et al., 2006; Figure 5G). Conversely, genes associated with the 2016; Zeisel et al., 2015), barcoded nuclei were automatically pro-survival receptors TrkB (Ntrk2) and TrkC (Ntrk3) and sorti- classified into 13 clusters, which were then aggregated into six lin-related VPS10 domain-containing receptor 3 (Sorcs3) were major cell classes: excitatory neurons, interneurons, oligoden- upregulated exclusively in deep cells (Huang and Reichardt, drocytes, microglia, oligodendrocyte precursor cells, and astro- 2001; Pfisterer and Khodosevich, 2017; Figure S6B). Other cytes (Figures S4A and S4B; see Figure S4C for control versus genes regulating neuronal vulnerability to insults were found epileptic and Figure S4D for unique molecular identifier (UMI) to be deregulated in both populations of CA1 pyramidal neu- levels and number of genes per cell type; Table S3). rons. These include the anti-apoptotic genes Akt3, Chl1, and We focused on pyramidal CA1 neurons, which represented Spag5 and the pro-apoptotic gene Nkain2 (He et al., 2020; 80% of the cell population (5,347 nuclei: 2,934 control and Katic et al., 2017; Tawarayama et al., 2018; Figures 5G and 2,413 epileptic). After two clustering rounds (Figure S4E), we S6B). identified deep and superficial pyramidal cells from control For some transcripts, we noted significant differential expres- (1,123 deep, 1,810 superficial) and epileptic (878 deep and sion in epilepsy-specific Pyr_ES cells and chose to focus on this 1,469 superficial) samples and a small population of pyramidal small cluster (Figure 5B, green). Compared with control pyrami- nuclei that were mostly epilepsy-specific (Pyr_ES; 67 nuclei: 66 dal cells, Pyr_ES cells showed transcriptional dysregulation of from epileptic CA1, 1 from control CA1) (Figure 5B). Notably, the GO families that were altered in superficial and deep CA1 clustering exhibited consistent distribution of markers for neurons in epilepsy (Figure S6C). To avoid confounding effects deep (e.g., Ndst4 and Col11a1) and superficial pyramidal neu- caused by the small sample size, we focused on transcripts iden- rons (e.g., Calb1 and Epha3) (Figures S4F and S4G). The spe- tified as differentially expressed between Pyr_ES and all other cific population of Pyr_ES cells showed total UMI values well epileptic cells that were well represented in both populations above the quality threshold (>500 UMIs) and a number of anno- (i.e., detected in at least 50% of cells per group; pct.1 > 0.5 tated genes lying within the control and epileptic ranges (Fig- and pct.2 > 0.5 over test result). Notably, many of these tran- ure S4H), excluding potential artifacts. Differential expression scripts were upregulated more in Pyr_ES neurons compared analysis of Pyr_ES cells against all detected cells, excluding with other epileptic cells, including the previously mentioned deep and superficial pyramidal neurons, revealed enrichment apoptotic-related genes Spag5, Chl1, and NKain2 (Figures 5G, of pyramidal CA1 marker genes such as Gria2, Rasgrf1, Camkv, 5H, and S6B). and Brd9, along with other transcripts enriched in excitatory To validate these findings at the single-cell level in situ, we neurons, including Epha6, Hs6st3, Cntnap2, Kcnh7, Kcnip4, used specific RNA-based probes (RNAscope) for sublayer-spe- and Meg3, which confirmed their pyramidal nature (Figures cific (Dcc) and group-specific (Dapk1 and Spag5) transcripts in S4I and S4J). Similar to our LCM RNA-seq observations at combination with Wfs1 immunostaining to label CA1 pyramidal the tissue level, we found consistent correlation of sublayer- cells in coronal brain sections of control and epileptic mice (n = specific genes in control and epileptic samples (Figures 5C, 5 control, n = 5 epileptic mice). As expected, we observed fluo- 5D, S5). rescence-labeled spots restricted to the somata of CA1 pyrami- Next, we sought to identify gene programs underlying the dal neurons (Figures 5I, 5J, S6D, and S6E). Quantification of differential vulnerability of deep and superficial CA1 pyramidal RNAscope revealed sublayer effects of Dcc (2-way ANOVA cell in response to epilepsy (adj. p < 0.05). Although statistical F(1,8) = 21.59, p = 0.0017) with significant post hoc differences evidence of this analysis is limited by the small sample size, it between layers in both groups (p < 0.05; see also distribution provided useful prospective hypotheses and identified many per cells in Figure 5K, small dots). In contrast, for Dapk1 epilepsy-related transcriptional changes. Many transcripts (F(1,8) = 12.99, p = 0.0069) and Spag5 (F(1,8) = 10.53, p = common to both cell types were associated with modulation 0.0118), we found that significant differences between groups of synaptic transmission, synapse organization, and regula- involved mainly the superficial sublayer, in line with our snRNA- tion of membrane potential (Figures 5E and 5F, common seq results (Figure 5K). genes in black). Single-nucleus differential expression anal- Our analysis confirmed profound transcriptional changes ysis also revealed striking differences between superficial related to neuronal excitability and neurodegeneration in and deep cells of genes involved in synaptic signaling and epileptic CA1 pyramidal neurons. These functions appear to ion membrane transport (Figures 5E and 5F). To exclude ef- be segregated across deep and superficial cell types, consis- fects of different sample sizes between groups and cell types, tent with electrophysiological, LCM RNA-seq, and histological we evaluated the top-ranked 300 DEGs and found similar GO data. Notably, snRNA-seq analysis revealed disease-associ- families to be involved differentially (Figure S6A). Deregulated ated genes that were highly cell type specific, including some genes in superficial cells in epilepsy showed a strong func- with key functions in neuronal survival. Hence, epilepsy- tional association that was related to the gamma-aminobuty- induced cell-type-specific regulation of genes with key roles ric acid (GABA) signaling pathway, calcium ion transport, po- in synaptic plasticity and survival is a potential molecular link tassium ion transport, regulation of glutamatergic synaptic between differential activity burden and vulnerability across transmission, axon growth and guidance, and synapse as- the CA1 radial axis, leading to regionalized hippocampal sembly. We also noted cell-type-specific changes in genes sclerosis. 8 Cell Reports 35, 109229, June 8, 2021 ll Article OPEN ACCESS Figure 5. snRNA-seq profiling of the normal and epileptic CA1 area (A) Nuclei were isolated from the CA1 region of adult mice and purified by flow cytometry for single-nucleus RNA-seq (snRNA-seq). (B) Uniform manifold approximation and projection (UMAP) plots of CA1 Pyrs subtypes segregated by condition (control and epileptic). Pyr_ES, epilepsy specific. (C) Heatmap showing normalized expression for principal gene markers for deep and Sup neurons in control and epileptic mice (96 enriched genes: 62 deep, 34 Sup; absolute log FC > 0.25; min.pct = 0.5; adj. p < 1030). Note the presence of a subset of cells with differential gene expression corresponding to the Pyr_ES cell population (green). (D) Significantly enriched genes in Sup and deep CA1 neurons in control mice and their relative level of enrichment in epileptic animals (adj. p < 0.05, 493 genes). Gene expression levels of CA1 neuronal subtype-specific enriched genes are preserved in epilepsy. (E) Venn diagram and heatmap of DEGs upregulated in epilepsy (adj. p < 0.05, Wilcoxon rank sum test). The heatmap shows normalized expression levels for upregulated DEGs that are common (34 genes), specific to deep (29 genes), or Sup (59 genes) CA1 neurons. Representative genes and associated significant GO terms are shown. FDR *adj. p < 0.05 (Fisher’s exact test). (F) Same as in (E) for downregulated genes. (G) Violin plot of normalized expression value for selected genes in deep and Sup neurons in epilepsy (blue) and control (black). Wilcoxon rank-sum test; *p < 0.05, **p < 1010, ***p < 1050, #p < 10100. Expression levels in Pyr_ES cells are also shown (green). (H) FC (left) and significance (right) for most upregulated genes in Pyr_ES compared with epileptic CA1 Pyrs (absolute log FC > 0.5 and min.pct = 0.5). (I) Immunostaining against the CA1 marker Wfs1 and multiplexed RNAscope for Dcc and Spag5 transcripts. Cells with their somata in the confocal plane are outlined. Note significant accumulation of Spag5 in an epileptic cell (green arrowhead). Scale bar, 10 mm. (J) Same as in (I) for Dapk1 and Spag5 transcripts. (K) Quantification of Dcc, Dapk1, and Spag5 per cell (small dots) and mouse (5 control and 5 epileptic); post hoc paired t tests, *p < 0.05. Cell Reports 35, 109229, June 8, 2021 9 ll OPEN ACCESS Article Pseudotemporal ordering reveals sublayer-specific DISCUSSION neurodegenerative progression Based on the data above, we speculated that epilepsy-related Our work identifies sublayer-specific transcriptional changes in responses and degenerative signals might evolve distinctly experimental hippocampal sclerosis. Using a combination of across cell types and sublayers. This may be especially critical techniques, we found that acquired TLE involves heterogeneous for cell death programs accompanying hippocampal sclerosis. biological processes running across deep and superficial CA1 To glean insights into these transitional states, we used manifold pyramidal cells. Our data suggest that epileptogenesis and the learning leveraged in nearest-neighbor information to automati- accompanying hippocampal sclerosis are evolving processes cally organize cells in trajectories along a principal tree, reflecting that affect the transcriptional state of neurons in a cell-type-spe- progression of associated biological processes (Qiu et al., cific manner. 2017a). Hippocampal sclerosis is a heterogeneous histopathological Low dimensional embedding of the automatically learned un- entity (Blu€mcke et al., 2013). The prognosis value of certain sub- derlying trajectory produced a spanning tree, revealing a topo- types has been debated, but emerging data suggest that it de- logical structure with four main branches (A, B, D, and E) and serves further consideration in light of cell type specificity (Prada two bifurcation points (1 and 2; Figure 6A). Most nuclei at Jardim et al., 2018). Superficial and deep CA1 pyramidal cells are branches A and B were from cells identified as deep, whereas innervated differently by local-circuit GABAergic interneurons, branches D and E were populated mostly by superficial cells express different neuromodulatory receptors, project differen- (Figure 6B). Interestingly, control and epileptic cells were distrib- tially to cortical and subcortical regions, and participate distinctly uted differently. Branches A and D were mostly populated by in hippocampal oscillations (Cembrowski et al., 2016a; Soltesz control pyramidal cells, whereas branches B and E contained and Losonczy, 2018; Valero and de la Prida, 2018). The more the vast majority of epileptic cells (Figures 6C and 6D). The tra- sustained excitability of superficial pyramidal cells (in terms of jectory topology suggested that the transcriptional state of single post-event depolarization and increased firing rate) and their cells progresses along sublayer-specific disease trajectories consistent stronger expression of the immediate-early gene- from a basal state (branches A and D) to the epileptic condition associated protein c-Fos after individual seizures suggest that (branches B and E) (Figure 6C). Notably, epilepsy-specific they have a larger activity burden compared with deep cells (Fig- Pyr_ES cells were retrieved at the end of the superficial epileptic ure 1). Although our snRNA-seq analysis identified many tran- E path, suggesting that they may be superficial neurons at a ter- scripts regulated differentially between deep and superficial minal pathological state (Figure 6B, green arrowhead). CA1 neurons that may underlie pro-epileptogenic changes in To identify mechanisms underlying single-cell transcriptional intrinsic excitability, microcircuit mechanisms are key contribu- changes from health to disease, we performed quantitative com- tors to these differences. First, perisomatic inhibition by parval- parison of gene expression kinetics across trajectories in each bumin (PV) basket cells is remarkably higher in deep compared cellular population (DDRTree method) (Qiu et al., 2017b). By with superficial cells (Soltesz and Losonczy, 2018; Valero et al., comparing the top-ranked 250 DEGs in deep (maximum false dis- 2015). Second, superficial cells are more likely to be driven by covery rate [FDR] adj. p = 8.94 3 104, likelihood ratio test) versus presynaptic CA3 activity than deep cells (Valero et al., 2015). superficial cells (maximum FDR adj. p = 5.8 3 1019), we found Similarly, direct inputs from layer III pyramidal neurons in the about 50% of genes (124 genes) associated with epilepsy car- medial entorhinal cortex are biased by deep CA1 cells, whereas rying significant cell-type-specific trajectories. Consistent with projections from the lateral entorhinal cortex are stronger in su- our previous results, we found that genes deregulated in epilepsy perficial cells (Masurkar et al., 2017). Given loss of medial layer in superficial CA1 neurons show tight functional association III inputs and sprouting of lateral inputs together with specific among them, leading to a high proportion of significantly enriched cannabinoid type 1 receptors pathways (Laurent et al., 2015; ontology terms compared with genes deregulated in deep CA1 Maroso et al., 2016), it is very likely that these changes contribute neurons (Figures 6E and 6F). Again, significantly enriched GO to more sustained activation of superficial cells in the epilepsy terms were related to structural and functional neuronal plasticity. context. Interestingly, pseudotemporal kinetics also revealed marked dif- Transcriptional dysregulation is a central feature of most ferences in the progression across sublayers of genes with impor- neurodegenerative diseases. Our transcriptional profiling in tant functions in regulating neuronal survival. For instance, the laser-microdissected CA1 deep and superficial sublayers (LCM gene transcript levels of Dapk1 showed a robust increase in su- RNA-seq) revealed strong regionalized responses to acquired perficial epileptic cells, whereas the gene expression levels of epilepsy. We found significantly larger transcriptional changes the pro-survival receptors TrkB (Ntrk2) and TrkB (Ntrk3), along of GO terms associated with immune and inflammatory pro- with other genes with key functions in neurons, such as Sorcs3 cesses, cytokine production, and programmed cell death in and Aff3, were upregulated specifically in deep pyramidal cells the superficial sublayer (Figure 2). This includes Cx3cr1, (Figures 6G and S6F). Other genes also related to neuronal vulner- P2ry12, and Il18, among other genes described previously, play- ability, such as Camk2a and Nrg3, displayed similar changes ing major roles in microglial activation and neurodegeneration along the disease trajectory in both cell types (Figures 6G and associated with epilepsy (Ali et al., 2015; Eyo et al., 2014; Vezzani S6F). Strikingly, the subset of apoptosis-related genes Spag5, and Baram, 2007). For the deep pyramidal sublayer, transcrip- Matk, and Chl1 reached their maximal expression level in Pyr_ES tional changes were milder in general and mostly involved dereg- at the end of the transcriptional trajectory in superficial cells (Fig- ulation of genes associated with changes in excitability, homeo- ures 6G, 6H, and S6F, green). static regulation, and synaptic signaling, such as Kcnh7, Kcnip4, 10 Cell Reports 35, 109229, June 8, 2021 ll Article OPEN ACCESS Figure 6. Pseudotemporal ordering of transcriptional changes unfolds separate disease trajectories across CA1 sublayers (A) UMAP plot of pseudotemporal trajectory analysis of Pyr subpopulations along with discriminative dimensionality reduction via learning a tree(DDRTree) coordinates uncovers four major branches (A, B, D, and E) and segregated from two branchpoints (1 and 2). (B) Distribution of Pyr subpopulations (deep, Sup, and Pyr_ES) across pseudotime trajectories. An arrowhead points to Pyr_ES cells located at the endpoint of branch E. (C) Distribution of Pyr subpopulations across pseudotime trajectories by condition (control and epileptic). (D) Distribution of Pyr subtypes across the trajectory topology. (E) Heatmap of the top-ranked branched expression analysis modeling (BEAM) significant changes through progression across the trajectory from branch A (basal state) to B (epilepsy) (deep cells). Each row indicates the standardized kinetic curves of a gene. (F) Same as (E) for the Sup-specific trajectory from basal state (branch D) to the epileptic condition (branch E). (G) Pseudotime kinetics across disease trajectory for selected DEGs in epilepsy. (H) Schematic of disease trajectories from control to epilepsy. Scn3b, and Gabrb3 (Dixit et al., 2016; Steinlein and Noebels, transcriptional changes affecting cell-cell communication (i.e., 2000). Many of these changes were confirmed by snRNA-seq neuron-microglia interactions). Such different latent trends might analysis (Figure 5), which disclosed significant upregulation of suggest that deep and superficial cells are at different stages in the Dapk1 specifically in superficial CA1 neurons (Bialik and the neurodegenerative process or being affected by pathways Kimchi, 2006), whereas genes associated with the pro-survival leading to differential vulnerability. For instance, although we receptors TrkB (Ntrk2) and TrkC (Ntrk3) and Sorcs3 were mostly found deregulated expression of genes associated with upregulated in deep cells (Huang and Reichardt, 2001; Pfisterer apoptotic pathways that were common to both cell types in and Khodosevich, 2017). response to epilepsy (i.e., Akt3, Nrg3, and Camk2a) (Ashpole By leveraging unsupervised branched expression learning, we and Hudmon, 2011; He et al., 2020; Kanekura et al., 2005; Ta- identified deviating transcriptional disease trajectories in deep warayama et al., 2018; Vest et al., 2010), genes related to cell and superficial pyramidal cells (Figure 6). Some of the GO terms death were found to be deregulated specifically in superficial involved suggest that transition from the normal to the epileptic or deep neurons. In superficial cells, epilepsy resulted in condition in deep cells is more likely associated with changes decreased expression of brain-derived neurotrophic factor in synaptic signaling and excitability (i.e., they reflect microcircuit (Bdnf) and increased transcript levels of the pro-apoptotic alterations), whereas superficial cells move differentially along pathway activator Dapk1 (Bialik and Kimchi, 2006). Conversely, Cell Reports 35, 109229, June 8, 2021 11 ll OPEN ACCESS Article in deep pyramidal neurons, epilepsy led to a robust increase in B LCM RNA-seq single-cell informed data analysis expression of the gene encoding the pro-survival receptors B Single-nuclei isolation TrkB (Ntrk2) and TrkC (Ntrk3) (Huang and Reichardt, 2001; Pfis- B Single-nucleus RNA sequencing terer and Khodosevich, 2017). Interestingly, pseudotime analysis B snRNA-seq analysis mapped some control cells along the epileptic branches in both B Cell trajectories and pseudotime analysis sublayers, suggesting that latent sublayer-specific transcrip- B Cell-type immunostaining and analysis tional processes might actually be running along life (Bishop B FluoroJade staining et al., 2010). B In situ hybridization analysis Our snRNA-seq analysis also disclosed an epilepsy-specific B RNAscope analysis pyramidal cell population, Pyr_ES, which accumulated at the d QUANTIFICATION AND STATISTICAL ANALYSIS end of the superficial trajectory branch (Figure 6B, arrowhead). This subset of pyramidal cells displayed remarkable expression SUPPLEMENTAL INFORMATION of neurodegeneration-related transcripts such as Cdk5, Ckb, Matk, Chl1, and Spag5 (Patrick et al., 1999). The presence of Supplemental information can be found online at https://doi.org/10.1016/j. celrep.2021.109229. cells with very high numbers of Spag5 molecules in the superficial sublayer of the epileptic CA1 was confirmed by com- ACKNOWLEDGMENTS bined immunofluorescence and single-molecule amplification methods. Our results indicate that genes dysregulated more in We thank Antonio Caler, Laura Dolón, Ester Lara, Beatriz Lázaro, and Juan epilepsy are expressed specifically by a vulnerable subset of py- Moriano for technical help. This work was supported by grants from MICINN ramidal neurons undergoing later stages of neurodegeneration (RTI2018-098581-B-I00 to L.M.P.), Fundación Tatiana Pérez de Guzman el by the time of sampling. We propose that Pyr_ES cells reflect Bueno, and the SynCogDis Network (SAF2014-52624-REDT and SAF2017- 90664-REDT to L.M.P. and A. Bayes). Collaboration between L.M.d.l.P. and the accumulated pro-epileptic transitional changes leading to Y.H. was supported by Human Frontiers Science Program (HFSP) grant epileptogenesis and neurodegeneration, as suggested by their RGP0022/2013. J.P.L.-A. was supported by grants from MICIU co-financed extreme location along disease trajectories. by ERDF (RYC-2015-18056 and RTI2018-102260-B-I00) and Severo Ochoa Our results identify previously unobserved heterogeneity in the grant SEV-2017-0723. R.R.-V. and A. Bayes were supported by MINECO neuronal patterns of activity of deep and superficial CA1 pyrami- BFU2015-69717-P and RTI2018-097037-B-100 and a Marie Curie career inte- dal neurons in epilepsy, show specific gene expression signa- gration grant (ref. 304111). A.V.M. was supported by MICINN (SAF2017- tures across CA1 deep and superficial sublayers that are associ- 85717-R) and Fundación Alicia Koplowitz. A. Barco was supported by grants SAF2017-87928-R from MICINN co-financed by ERDF and RGP0039/2017 ated with neuronal loss and hippocampal sclerosis, reveal from the Human Frontiers Science Program Organization. The Instituto de disease trajectories of deep and superficial CA1 pyramidal neu- Neurociencias is a ‘‘Centre of Excellence Severo Ochoa.’’ D.G.-D. and rons in epilepsy, and uncover the underlying transcriptional pro- C.M.N. hold PhD fellowships from MICINN (BES-2013-064171 and BES- grams. By dissecting the transcriptional landscape across CA1 2016-076281, respectively). sublayers in epilepsy, our work offers insights into the mecha- nisms regulating epileptogenesis and highlights the importance AUTHOR CONTRIBUTIONS of leveraging cell type specificity to better understand the Conceptualization, L.M.d.l.P. and J.P.L.-A.; methodology, L.M.d.l.P., J.P.L.- phenotypic complexities accompanying hippocampal sclerosis A., A. Bayés, and A. Barco; investigation, E.C., M.V., B.G., D.C.M., C.M.N., in epilepsy. L.B.-E., R.R.-V., A.V.M., I.F.-L., and D.G.-D.; formal analysis, E.C., A.M.-G., M.V., J.P.L.-A., and L.M.d.l.P.; software, A.M.-G., and J.P.L.-A.; resources, STAR+METHODS Y.H., and M.S.; writing – original draft, L.M.d.l.P. and J.P.L.-A.; writing – re- viewing & editing, L.M.d.l.P., J.P.L.-A., E.C., and A.M.-G.; funding acquisition, L.M.d.l.P. and J.P.L.-A.; supervision, L.M.d.l.P. and J.P.L.-A. Detailed methods are provided in the online version of this paper and include the following: DECLARATION OF INTERESTS d KEY RESOURCES TABLE The authors declare no competing interests. d RESOURCE AVAILABILITY B Lead contact Received: September 16, 2020 B Materials availability Revised: March 18, 2021 B Data and code availability Accepted: May 17, 2021 d EXPERIMENTAL MODEL AND SUBJECT DETAILS Published: June 8, 2021 d METHOD DETAILS REFERENCES B EEG recordings and analysis B In vivo recording and labeling of single cells Ali, I., Chugh, D., and Ekdahl, C.T. (2015). Role of fractalkine-CX3CR1 pathway B Analysis of intracellular single-cell recordings in seizure-induced microglial activation, neurodegeneration, and neuroblast B cFos immunostaining and analysis production in the adult rat brain. Neurobiol. Dis. 74, 194–203. B Laser capture microdissection (LCM) and RNA isola- Anders, S., Pyl, P.T., and Huber, W. (2015). HTSeq–a Python framework to tion work with high-throughput sequencing data. Bioinformatics 31, 166–169. B LCM RNA-seq library construction and sequencing Ashpole, N.M., and Hudmon, A. (2011). Excitotoxic neuroprotection and B LCM RNA-seq data analysis vulnerability with CaMKII inhibition. Mol. Cell. Neurosci. 46, 720–730. 12 Cell Reports 35, 109229, June 8, 2021 ll Article OPEN ACCESS Benjamini, Y., and Hochberg, Y. (1995). Controlling the False Discovery Rate: NMDA receptors and microglial P2Y12 receptors after status epilepticus. A Practical and Powerful Approach to Multiple Testing. J. R. Stat. Soc. B Meth- J. Neurosci. 34, 10528–10540. odol. 57, 289–300. Habib, N., Li, Y., Heidenreich, M., Swiech, L., Avraham-Davidi, I., Trombetta, Bialik, S., and Kimchi, A. (2006). The death-associated protein kinases: struc- J.J., Hession, C., Zhang, F., and Regev, A. (2016). Div-Seq: Single-nucleus ture, function, and beyond. Annu. Rev. Biochem. 75, 189–210. RNA-Seq reveals dynamics of rare adult newborn neurons. Science 353, Biber, K., Neumann, H., Inoue, K., and Boddeke, H.W.G.M. (2007). Neuronal 925–928. ‘On’ and ‘Off’ signals control microglia. Trends Neurosci. 30, 596–602. He, J., Green, A.R., Li, Y., Chan, S.Y.T., and Liu, D.X. (2020). SPAG5: An Bishop, N.A., Lu, T., and Yankner, B.A. (2010). Neural mechanisms of ageing Emerging Oncogene. Trends Cancer 6, 543–547. and cognitive decline. Nature 464, 529–535. Hickman, S.E., Kingery, N.D., Ohsumi, T.K., Borowsky, M.L., Wang, L.C., €mcke, I., Kistner, I., Clusmann, H., Schramm, J., Becker, A.J., Elger, C.E., Blu Means, T.K., and El Khoury, J. (2013). The microglial sensome revealed by Bien, C.G., Merschhemke, M., Meencke, H.-J., Lehmann, T., et al. (2009). To- direct RNA sequencing. Nat. Neurosci. 16, 1896–1905. wards a clinico-pathological classification of granule cell dispersion in human Hoek, R.H., Ruuls, S.R., Murphy, C.A., Wright, G.J., Goddard, R., Zurawski, mesial temporal lobe epilepsies. Acta Neuropathol. 117, 535–544. S.M., Blom, B., Homola, M.E., Streit, W.J., Brown, M.H., et al. (2000). Down- €mcke, I., Thom, M., Aronica, E., Armstrong, D.D., Bartolomei, F., Bernas- Blu regulation of the macrophage lineage through interaction with OX2 (CD200). coni, A., Bernasconi, N., Bien, C.G., Cendes, F., Coras, R., et al. (2013). Inter- Science 290, 1768–1771. national consensus classification of hippocampal sclerosis in temporal lobe Huang, E.J., and Reichardt, L.F. (2001). Neurotrophins: roles in neuronal devel- epilepsy: a Task Force report from the ILAE Commission on Diagnostic opment and function. Annu. Rev. Neurosci. 24, 677–736. Methods. Epilepsia 54, 1315–1329. Huang, W., Sherman, B.T., and Lempicki, R.A. (2009). Systematic and integra- Butler, A., Hoffman, P., Smibert, P., Papalexi, E., and Satija, R. (2018). Inte- tive analysis of large gene lists using DAVID bioinformatics resources. Nat. grating single-cell transcriptomic data across different conditions, technolo- Protoc. 4, 44–57. gies, and species. Nat. Biotechnol. 36, 411–420. Huber, W., Carey, V.J., Gentleman, R., Anders, S., Carlson, M., Carvalho, B.S., Cembrowski, M.S., Bachman, J.L., Wang, L., Sugino, K., Shields, B.C., and Bravo, H.C., Davis, S., Gatto, L., Girke, T., et al. (2015). Orchestrating high- Spruston, N. (2016a). Spatial Gene-Expression Gradients Underlie Prominent throughput genomic analysis with Bioconductor. Nat. Methods 12, 115–121. Heterogeneity of CA1 Pyramidal Neurons. Neuron 89, 351–368. Iacopino, A.M., and Christakos, S. (1990). Specific reduction of calcium-bind- Cembrowski, M.S., Wang, L., Sugino, K., Shields, B.C., and Spruston, N. ing protein (28-kilodalton calbindin-D) gene expression in aging and neurode- (2016b). Hipposeq: a comprehensive RNA-seq database of gene expression generative diseases. Proc. Natl. Acad. Sci. USA 87, 4078–4082. in hippocampal principal neurons. eLife 5, e14997. Chen, Z., Brodie, M.J., Liew, D., and Kwan, P. (2018). Treatment Outcomes in Inostroza, M., Cid, E., Brotons-Mas, J., Gal, B., Aivar, P., Uzcategui, Y.G., Patients With Newly Diagnosed Epilepsy Treated With Established and New Sandi, C., and Menendez de la Prida, L. (2011). Hippocampal-dependent Antiepileptic Drugs: A 30-Year Longitudinal Cohort Study. JAMA Neurol. 75, spatial memory in the water maze is preserved in an experimental model of 279–286. temporal lobe epilepsy in rats. PLoS ONE 6, e22372. Cid, E., and de la Prida, L.M. (2019). Methods for single-cell recording and la- Jefferys, J.G.R., Menendez de la Prida, L., Wendling, F., Bragin, A., Avoli, M., beling in vivo. J. Neurosci. Methods 325, 108354. Timofeev, I., and Lopes da Silva, F.H. (2012). Mechanisms of physiological and epileptic HFO generation. Prog. Neurobiol. 98, 250–264. Coras, R., Pauli, E., Li, J., Schwarz, M., Rössler, K., Buchfelder, M., Hamer, H., Stefan, H., and Blumcke, I. (2014). Differential influence of hippocampal sub- Kanekura, K., Hashimoto, Y., Kita, Y., Sasabe, J., Aiso, S., Nishimoto, I., and fields to memory formation: insights from patients with temporal lobe epilepsy. Matsuoka, M. (2005). A Rac1/phosphatidylinositol 3-kinase/Akt3 anti- Brain 137, 1945–1957. apoptotic pathway, triggered by AlsinLF, the product of the ALS2 gene, antag- onizes Cu/Zn-superoxide dismutase (SOD1) mutant-induced motoneuronal Curia, G., Longo, D., Biagini, G., Jones, R.S.G., and Avoli, M. (2008). The pilo- cell death. J. Biol. Chem. 280, 4532–4543. carpine model of temporal lobe epilepsy. J. Neurosci. Methods 172, 143–157. Katic, J., Loers, G., Tosic, J., Schachner, M., and Kleene, R. (2017). The cell Davies, K.G., Hermann, B.P., Dohan, F.C., Jr., Foley, K.T., Bush, A.J., and Wy- adhesion molecule CHL1 interacts with patched-1 to regulate apoptosis dur- ler, A.R. (1996). Relationship of hippocampal sclerosis to duration and age of ing postnatal cerebellar development. J. Cell Sci. 130, 2606–2619. onset of epilepsy, and childhood febrile seizures in temporal lobectomy pa- tients. Epilepsy Res. 24, 119–126. Kiselev, V.Y., Kirschner, K., Schaub, M.T., Andrews, T., Yiu, A., Chandra, T., de Lanerolle, N.C., Kim, J.H., Williamson, A., Spencer, S.S., Zaveri, H.P., Eid, Natarajan, K.N., Reik, W., Barahona, M., Green, A.R., and Hemberg, M. T., and Spencer, D.D. (2003). A retrospective analysis of hippocampal pathol- (2017). SC3: consensus clustering of single-cell RNA-seq data. Nat. Methods ogy in human temporal lobe epilepsy: evidence for distinctive patient subcat- 14, 483–486. egories. Epilepsia 44, 677–687. Laurent, F., Brotons-Mas, J.R., Cid, E., Lopez-Pigozzi, D., Valero, M., Gal, B., Dixit, A.B., Banerjee, J., Srivastava, A., Tripathi, M., Sarkar, C., Kakkar, A., and de la Prida, L.M. (2015). Proximodistal structure of theta coordination in Jain, M., and Chandra, P.S. (2016). RNA-seq analysis of hippocampal tissues the dorsal hippocampus of epileptic rats. J. Neurosci. 35, 4760–4775. reveals novel candidate genes for drug refractory epilepsy in patients with Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., Marth, G., MTLE-HS. Genomics 107, 178–188. Abecasis, G., and Durbin, R.; 1000 Genome Project Data Processing Sub- Dobin, A., Davis, C.A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., Batut, group (2009). The Sequence Alignment/Map format and SAMtools. Bioinfor- P., Chaisson, M., and Gingeras, T.R. (2013). STAR: ultrafast universal RNA-seq matics 25, 2078–2079. aligner. Bioinformatics 29, 15–21. Liao, Y., Wang, J., Jaehnig, E.J., Shi, Z., and Zhang, B. (2019). WebGestalt Dong, H.-W., Swanson, L.W., Chen, L., Fanselow, M.S., and Toga, A.W. 2019: gene set analysis toolkit with revamped UIs and APIs. Nucleic Acids (2009). Genomic-anatomic evidence for distinct functional domains in hippo- Res. 47, W199–W205. campal field CA1. Proc. Natl. Acad. Sci. USA 106, 11794–11799. Love, M.I., Huber, W., and Anders, S. (2014). Moderated estimation of fold Du, F., Eid, T., Lothman, E.W., Köhler, C., and Schwarcz, R. (1995). Preferential change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. neuronal loss in layer III of the medial entorhinal cortex in rat models of tempo- Luo, J., Elwood, F., Britschgi, M., Villeda, S., Zhang, H., Ding, Z., Zhu, L., ral lobe epilepsy. J. Neurosci. 15, 6301–6313. Alabsi, H., Getachew, R., Narasimhan, R., et al. (2013). Colony-stimulating fac- Eyo, U.B., Peng, J., Swiatkowski, P., Mukherjee, A., Bispo, A., and Wu, L.J. tor 1 receptor (CSF1R) signaling in injured neurons facilitates protection and (2014). Neuronal hyperactivity recruits microglial processes via neuronal survival. J. Exp. Med. 210, 157–172. Cell Reports 35, 109229, June 8, 2021 13 ll OPEN ACCESS Article Maroso, M., Szabo, G.G., Kim, H.K., Alexander, A., Bui, A.D., Lee, S.-H., Lutz, Soltesz, I., and Losonczy, A. (2018). CA1 pyramidal cell diversity enabling par- B., and Soltesz, I. (2016). Cannabinoid Control of Learning and Memory allel information processing in the hippocampus. Nat. Neurosci. 21, 484–493. through HCN Channels. Neuron 89, 1059–1073. Steinlein, O.K., and Noebels, J.L. (2000). Ion channels and epilepsy in man and Masurkar, A.V., Srinivas, K.V., Brann, D.H., Warren, R., Lowes, D.C., and Sie- mouse. Curr. Opin. Genet. Dev. 10, 286–291. gelbaum, S.A. (2017). Medial and Lateral Entorhinal Cortex Differentially Excite Stuart, T., Butler, A., Hoffman, P., Hafemeister, C., Papalexi, E., Mauck, W.M., Deep versus Superficial CA1 Pyramidal Neurons. Cell Rep. 18, 148–160. 3rd, Hao, Y., Stoeckius, M., Smibert, P., and Satija, R. (2019). Comprehensive Mathern, G.W., Kuhlman, P.A., Mendoza, D., and Pretorius, J.K. (1997). Hu- Integration of Single-Cell Data. Cell 177, 1888–1902.e21. man fascia dentata anatomy and hippocampal neuron densities differ depend- Tawarayama, H., Yamada, H., Amin, R., Morita-Fujimura, Y., Cooper, H.M., ing on the epileptic syndrome and age at first seizure. J. Neuropathol. Exp. Shinmyo, Y., Kawata, M., Ikawa, S., and Tanaka, H. (2018). Draxin regulates Neurol. 56, 199–212. hippocampal neurogenesis in the postnatal dentate gyrus by inhibiting DCC- McInnes, L., Healy, J., and Melville, J. (2018). UMAP: Uniform Manifold induced apoptosis. Sci. Rep. 8, 840. Approximation and Projection for Dimension Reduction. arXiv, ar- Thom, M., Liagkouras, I., Martinian, L., Liu, J., Catarino, C.B., and Sisodiya, Xiv:1802.03426. https://arxiv.org/abs/1802.03426. S.M. (2012). Variability of sclerosis along the longitudinal hippocampal axis McKenzie, B.A., Dixit, V.M., and Power, C. (2020). Fiery Cell Death: Pyroptosis in epilepsy: a post mortem study. Epilepsy Res. 102, 45–59. in the Central Nervous System. Trends Neurosci. 43, 55–73. Thorvaldsdóttir, H., Robinson, J.T., and Mesirov, J.P. (2013). Integrative Geno- Morin-Brureau, M., Milior, G., Royer, J., Chali, F., Le Duigou, C., Savary, E., mics Viewer (IGV): high-performance genomics data visualization and explora- Blugeon, C., Jourdren, L., Akbar, D., Dupont, S., et al. (2018). Microglial phe- tion. Brief. Bioinform. 14, 178–192. notypes in the human epileptic temporal lobe. Brain 141, 3343–3360. Valero, M., and de la Prida, L.M. (2018). The hippocampus in depth: a sublayer- Morris, M.E., Baimbridge, K.G., el-Beheiry, H., Obrocea, G.V., and Rosen, A.S. specific perspective of entorhinal-hippocampal function. Curr. Opin. Neuro- (1995). Correlation of anoxic neuronal responses and calbindin-D28k localiza- biol. 52, 107–114. tion in stratum pyramidale of rat hippocampus. Hippocampus 5, 25–39. Valero, M., Cid, E., Averkin, R.G., Aguilar, J., Sanchez-Aguilera, A., Viney, T.J., Patrick, G.N., Zukerberg, L., Nikolic, M., de la Monte, S., Dikkes, P., and Tsai, Gomez-Dominguez, D., Bellistri, E., and de la Prida, L.M. (2015). Determinants L.H. (1999). Conversion of p35 to p25 deregulates Cdk5 activity and promotes of different deep and superficial CA1 pyramidal cell dynamics during sharp- neurodegeneration. Nature 402, 615–622. wave ripples. Nat. Neurosci. 18, 1281–1290. Pfisterer, U., and Khodosevich, K. (2017). Neuronal survival in the brain: neuron Valero, M., Averkin, R.G., Fernandez-Lamo, I., Aguilar, J., Lopez-Pigozzi, D., type-specific mechanisms. Cell Death Dis. 8, e2643. Brotons-Mas, J.R., Cid, E., Tamas, G., and Menendez de la Prida, L. (2017). Pfisterer, U., Petukhov, V., Demharter, S., Meichsner, J., Thompson, J.J., Ba- Mechanisms for Selective Single-Cell Reactivation during Offline Sharp- tiuk, M.Y., Martinez, A.A., Vasistha, N.A., Thakur, A., Mikkelsen, J., et al. Wave Ripples and Their Distortion by Fast Ripples. Neuron 94, 1234–1247.e7. (2020). Identification of epilepsy-associated neuronal subtypes and gene van der Maaten, L., and Hinton, G. (2009). Visualizing Data using t-SNE. expression underlying epileptogenesis. Nat. Commun. 11, 5038. J. Mach. Learn. Res. 9, 2579–2605. Prada Jardim, A., Liu, J., Baber, J., Michalak, Z., Reeves, C., Ellis, M., Novy, J., Vest, R.S., O’Leary, H., Coultrap, S.J., Kindy, M.S., and Bayer, K.U. (2010). de Tisi, J., McEvoy, A., Miserocchi, A., et al. (2018). Characterising subtypes of Effective post-insult neuroprotection by a novel Ca(2+)/ calmodulin-depen- hippocampal sclerosis and reorganization: correlation with pre and postoper- dent protein kinase II (CaMKII) inhibitor. J. Biol. Chem. 285, 20675–20682. ative memory deficit. Brain Pathol. 28, 143–154. Vezzani, A., and Baram, T.Z. (2007). New roles for interleukin-1 Beta in the Qiu, X., Mao, Q., Tang, Y., Wang, L., Chawla, R., Pliner, H.A., and Trapnell, C. mechanisms of epilepsy. Epilepsy Curr. 7, 45–50. (2017a). Reversed graph embedding resolves complex single-cell trajectories. Winden, K.D., Bragin, A., Engel, J., and Geschwind, D.H. (2015). Molecular al- Nat. Methods 14, 979–982. terations in areas generating fast ripples in an animal model of temporal lobe Qiu, X., Hill, A., Packer, J., Lin, D., Ma, Y.A., and Trapnell, C. (2017b). Single- epilepsy. Neurobiol. Dis. 78, 35–44. cell mRNA quantification and differential analysis with Census. Nat. Methods Wittner, L., Eross, L., Szabó, Z., Tóth, S., Czirják, S., Halász, P., Freund, T.F., 14, 309–315. and Maglóczky, Z.S. (2002). Synaptic reorganization of calbindin-positive neu- Rodrigues, G.R., Kandratavicius, L., Peixoto-Santos, J.E., Monteiro, M.R., rons in the human hippocampal CA1 region in temporal lobe epilepsy. Neuro- Gargaro, A.C., Geraldi, C. de V.L., Velasco, T.R., and Leite, J.P. (2015). science 115, 961–978. Increased frequency of hippocampal sclerosis ILAE type 2 in patients with Wyler, A.R., Curtis Dohan, F., Schweitzer, J.B., and Berry, A.D. (1992). A mesial temporal lobe epilepsy with normal episodic memory. Brain 138, e359. grading system for mesial temporal pathology (hippocampal sclerosis) from Saghafi, S., Ferguson, L., Hogue, O., Gales, J.M., Prayson, R., and Busch, anterior temporal lobectomy. J. Epilepsy 5, 220–225. R.M. (2018). Histopathologic subtype of hippocampal sclerosis and episodic You, J.C., Muralidharan, K., Park, J.W., Petrof, I., Pyfer, M.S., Corbett, B.F., memory performance before and after temporal lobectomy for epilepsy. Epi- LaFrancois, J.J., Zheng, Y., Zhang, X., Mohila, C.A., et al. (2017). Epigenetic lepsia 59, 825–833. suppression of hippocampal calbindin-D28k by DFosB drives seizure-related Scandaglia, M., Lopez-Atalaya, J.P., Medrano-Fernandez, A., Lopez-Cas- cognitive deficits. Nat. Med. 23, 1377–1383. cales, M.T., Del Blanco, B., Lipinski, M., Benito, E., Olivares, R., Iwase, S., Zeisel, A., Muñoz-Manchado, A.B., Codeluppi, S., Lönnerberg, P., Manno, G., Shi, Y., and Barco, A. (2017). Loss of Kdm5c Causes Spurious Transcription La Juréus, A., Marques, S., Munguba, H., He, L., Betsholtz, C., et al. (2015). and Prevents the Fine-Tuning of Activity-Regulated Enhancers in Neurons. Cell types in the mouse cortex and hippocampus revealed by single-cell Cell Rep. 21, 47–59. RNA-seq. Science 347, 1138–1142. Slomianka, L., Amrein, I., Knuesel, I., Sørensen, J.C., and Wolfer, D.P. (2011). Zeisel, A., Hochgerner, H., Lönnerberg, P., Johnsson, A., Memic, F., van der Hippocampal pyramidal cells: the reemergence of cortical lamination. Brain Zwan, J., Häring, M., Braun, E., Borm, L.E., La Manno, G., et al. (2018). Molec- Struct. Funct. 216, 301–317. ular Architecture of the Mouse Nervous System. Cell 174, 999–1014.e22. 14 Cell Reports 35, 109229, June 8, 2021 ll Article OPEN ACCESS STAR+METHODS KEY RESOURCES TABLE REAGENT or RESOURCE SOURCE IDENTIFIER Antibodies Rabbit anti-calbindin D-28k Swant Cat# CB38; RRID: AB_10000340 Mouse anti-calbindin D-28k Swant Cat# 300; RRID: AB_10000347 Rabbit anti-Wfs1 Proteintech Cat# 11558-1-AP; RRID: AB_2216046 Rabbit anti-Iba1 Wako Cat# 019-19741; RRID: AB_839504 Rabbit anti-GFAP Sigma Cat# G9269; RRID: AB_477035 Rabbit anti-Olig2 Millipore Cat# AB9610; RRID: AB_570666 Rabbit anti-cFos Santa Cruz Biotechnology Cat# sc-52; RRID: AB_2106783 Alexa Fluor633 goat anti-rabbit IgG Thermo Fisher Cat# A21070; RRID: AB_2535731 Alexa Fluor488 donkey anti-mouse IgG Thermo Fisher Cat# A21202; RRID: AB_141607 Alexa Fluor488 goat anti-mouse IgG Jackson Immunoresearch Cat# 115-545-003; RRID: AB_2338840 Rhodamine Red goat anti-mouse IgG Jackson Immunoresearch Cat# 115-295-003; RRID: AB_2338756 Alexa Fluor488-conjugated streptavidin Jackson Immunoresearch Cat# 016-540-084; RRID: AB_2337249 Chemicals, peptides, and recombinant proteins BSA Sigma A7906 EDTA Sigma 607-429-00-8 PBS VWR E404-200 TABS Tris-HCl Sigma T5941 NaCl VWR 7647-14-5 MgCl2 Sigma M2670 IGEPAL Sigma I8896 ARCTURUS PicoPure RNA isolation kit Thermo Fisher Scientific KIT0204 RNase-free DNase set QIAGEN 79254 Bisbenzimide H33258 Sigma-Aldrich Cat# B2883; CAS: 23491-45-4 Neurobiotin tracer Vector Labs Cat# SP-1120 FluoroJade C dye Sigma AG325 DPX mountant VWR 360294H ProLong Gold antifade mountant ThermoFisher Scientific P10144 Critical commercial assays Chromium single cell 30 library & gel bead kit v2 10x Genomics cat#: 120237 (v2) Chromium i7 multiplex kit 10x Genomics N/A Chromium single cell A chip kit 10x Genomics cat#: 120236 TruSeq stranded mRNA library prep kit Illumina cat#: RS-122-2101/2 HiSeq sequencing kit TruSeq v4 Illumina cat#: FC-401-4003 RNAscope fluorescent multiplex reagent kit ACD 320850 Deposited data LCM-RNA-seq Gene Expression Omnibus (NCBI) GSE143555 snRNA-seq Gene Expression Omnibus (NCBI) GSE143560 Experimental models: organisms/strains Mouse: Thy1.2-G-CaMP7-DsRed2 RIKEN Bioresources Center stock RBRCO6579 Mouse: Calb1-2A-dgCre-D Jackson Labs N/A Mouse: TdTomato reporter line Ai9 Jackson Labs Stock No:007905 Rat: Wistar Instituto Cajal Animal facility N/A Mouse: C57BL/6 Instituto Cajal Animal facility N/A (Continued on next page) Cell Reports 35, 109229, June 8, 2021 e1 ll OPEN ACCESS Article Continued REAGENT or RESOURCE SOURCE IDENTIFIER Oligonucleotides Wfs1-forw (50 CATCCCCAACAACCTGTTCC) Sigma NM_031823.1, bp 783-1631 Wfs1-rev (50 ACACCAGGTAGGGCACCACC) Sigma N/A Ndst4-forw (50 ATACATCCAAAACTGATCCCAC) Sigma XM_006233274.2, bp 536-1024 Ndst4-rev (50 AAAAGCACTGGCTGGTAGGTAG) Sigma N/A Syt17-forw (50 ACACATCCAAGTCCACATACAG) Sigma NM_138849.1, bp 378-1118 Syt17-rev (50 GTGTACCAGCTGGATTTTCACA) Sigma N/A Hrt1a-forw (50 CTTTCTATATCCCGCTGTTGCT) Sigma J05276.1, bp 730-1663 Hrt1a-rev (50 CACACAGATACTAGTGGTTCTC) Sigma N/A Scn7a-for (50 GATTATTTCTTCCCTTGTGTATGC) Sigma NM_131912.1, bp 2812-3550 Scn7a-rev (50 CTAAAGTCATCTTCGCTCAAGG) Sigma N/A Recombinant DNA Mm-Spag5 ACD 505691 Mm-Dcc-C3 ACD 427491-C3 Mm-Dapk1-C3 ACD 880221-C3 RNAscope 3-plex positive control probe-Mm ACD 320881 RNAscope 3-plex negative control probe ACD 320871 Software and algorithms Sequencer Software HiSeq Control Software Illumina, Inc https://support.illumina.com/sequencing/ sequencing_instruments/hiseq_2500/ downloads.html?langsel=/us R R Foundation for Statistical https://www.R-project.org/; RRID: SCR_001905 Computing RStudio RStudio, PBC https://www.rstudio.com/; RRID: SCR_000432 FastQC Simon Andrews https://www.bioinformatics.babraham.ac.uk/ projects/fastqc; RRID: SCR_014583 STAR Dobin et al., 2013 https://code.google.com/archive/p/rna-star; RRID: SCR_004463 Samtools Li et al., 2009 http://samtools.sourceforge.net; RRID: SCR_002105 HTSeq Anders et al., 2015 https://htseq.readthedocs.io/en/master/; RRID: SCR_005514 DESeq2 Love et al., 2014 http://www.bioconductor.org/packages/release/ bioc/html/DESeq2.html; RRID: SCR_015687 Cell Ranger Huang et al., 2009 https://david.ncifcrf.gov; RRID: SCR_017344 Seurat Butler et al., 2018; Stuart et al., 2019 http://satijalab.org/seurat/; RRID: SCR_016341 Monocle 2 Qiu et al., 2017a, 2017b http://cole-trapnell-lab.github.io/monocle-release; RRID: SCR_016339 DAVID Bioinformatics Huang et al., 2009 https://david.ncifcrf.gov; RRID: SCR_001881 Webgestalt Liao et al., 2019 http://www.webgestalt.org; RRID: SCR_006786 MATLAB 2016b Mathworks https://www.mathworks.com; RRID: SCR_001622 ImageJ NIH Image https://imagej.net/ImageJ; RRID: SCR_003070 SPSS Statistics IBM https://www.ibm.com/products/spss-statistics; RRID: SCR_002865 Other Resource website for the LCM RNA-seq data This paper http://lopezatalayalab.in.umh-csic.es/CA1_ Sublayers_&_Epilepsy Resource website for the snRNA-seq data This paper http://lopezatalayalab.in.umh-csic.es/CA1_ SingleNuclei_&_Epilepsy Silicon probes: 16-channel linear; 100 mm Neuronexus A1x16-5mm-100-413 inter-spacing; 413 mm2 electrode area e2 Cell Reports 35, 109229, June 8, 2021 ll Article OPEN ACCESS RESOURCE AVAILABILITY Lead contact Further information and requests for reagents and resources may be directed to the lead contact, Dr. Liset M de la Prida (lmprida@ cajal.csic.es). Materials availability This study did not generate any new unique reagents. Data and code availability No unpublished custom code, software or algorithm was used in this study. Freely available software and algorithms used for analysis are listed in the resource table. Some analyses were specifically designed for the purpose of this paper using routines written in MAT- LAB 7.10 (MathWorks). All custom scripts and data contained in this manuscript are available upon request from the Lead Contact. The accession number for the bulk-tissue transcriptome (LCM RNA-seq) and single-nuclei RNA-seq (snRNA-seq) data reported in this paper is Gene Expression Omnibus (GEO) database: GSE143555 (LCM RNA-seq) https://www.ncbi.nlm.nih.gov/geo/query/ acc.cgi?acc=GSE143555; GSE143560 (snRNA-seq) https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE143560. Two open web applications provide visualization of transcriptomic data: LCM RNA-seq data: http://lopezatalayalab.in.umh-csic.es/ CA1_Sublayers_&_Epilepsy. snRNA-seq data: http://lopezatalayalab.in.umh-csic.es/CA1_SingleNuclei_&_Epilepsy. EXPERIMENTAL MODEL AND SUBJECT DETAILS All experimental protocols and procedures were performed according to the Spanish legislation (R.D. 1201/2005 and L.32/2007), the European Communities Council Directives of 1986 (86/609/EEC) and 2003 (2003/65/CE) for animal research, and were approved by the Ethics Committee of the Instituto Cajal. Adult male Wistar rats (180–200 g; 6-8 weeks postnatal), as well as wild-type C57 and transgenic mice (20-25 g; 8-9 weeks pos- natal), were treated with multiple intraperitoneal injections of kainate (5 mg/kg) at hourly intervals until they reached status epilepticus. Transgenic mouse lines included the Thy1.2-G-CaMP7-DsRed2 (c57BL/6J-Tg(Thy1-G-CaMP7-DsRed2)492Bsi, stock RBRCO6579, RIKEN) and the Calb1-Cre (Jackson Lab, stock No:023531; Calb1-2A-dgCre-D) crossed with the tdTomato reporter line (Jackson Lab, stock No:007905; Ai9). The status was defined as a condition of continuous seizures lasting longer than 30 min. In a subset of rats, the lithium-pilocarpine model was used. These rats were i.p. injected with pilocarpine hydrochloride 12–24 h after the injection of lithium chloride (127 mg/kg, i.p.). Between one and four doses of 10 mg/kg pilocarpine were injected every 30 min until the status epilepticus was reached. Diazepam (4 mg/kg) was injected 1 h post-status to stop convulsions in all animals. They received intraperitoneal injections of 5% dextrose in saline (2.5 ml) and their diet was supplemented with fruit and powder milk during the following 2–3 days. After 3 days, animals behaved normally and were housed individually. Control animals were injected with saline and received treatments similar to epileptic animals. Animals used in this study were examined in the chronic phase of the status epilepticus model of TLE (2-23 weeks post-status) when they already exhibited spontaneous seizures and interictal discharges (IID) (Figures S1A and S1B). In epileptic rats and mice, these events were typically associated with high-frequency oscillations (HFOs; Figure S1B), which are considered biomarkers of epileptogenesis (Jefferys et al., 2012; Valero et al., 2017). We focused on the dorsal hippocampus given the major role in associ- ated cognitive comorbidities of epilepsy and more consistent neuronal loss as compared to ventral (Thom et al., 2012; Valero et al., 2017). In some epileptic animals, we aimed to induce convulsive seizures using high-pitched sounds (pulses of 95-100 dB and 1-20 s duration at random intervals of 1-20 s during 10 min). A convulsive motor seizure was defined according to clinical criteria, as sus- tained and repeated forelimbs automatisms with or without falling. Control animals did never exhibit convulsive seizures using this protocol. METHOD DETAILS EEG recordings and analysis To evaluate epileptogenesis, some rats and mice were implanted with either intracranial 16-channel silicon probes or skull EEG-grid of 32-channels (Neuronexus) under isoflurane anesthesia (1.5%–2% mixed in oxygen 400–800 ml/min). Jeweler’s screws were in- serted into the skull for providing additional anchoring and reference/ground connections (over the cerebellum). The implant was secured with dental cement. Animals were recovered from anesthesia and returned to home cages. For recordings, EEG signals were pre-amplified (4x gain) and recorded with a 32-channel AC amplifier (Axona), further amplified by 100, filtered by analog means at 1 Hz to 5 kHz, and sampled at 20 kHz/channel with 12 bits precision. EEG recordings were synchro- nized with a ceiling video camera (30 frames/s) to track the animal position in space. Analysis of electrophysiological signals was implemented in MATLAB 9.3 (MathWorks). EEG signals from the frontal and parietal electrodes were used to identify theta periods during running (band-pass 4-12 Hz) and periods of immobility characterized by low frequency delta activity (0.5-4 Hz). Forward-backward-zero-phase finite impulse response (FIR) filters of order 512 were used to Cell Reports 35, 109229, June 8, 2021 e3 ll OPEN ACCESS Article preserve temporal relationships between channels and signals. Spectral values fitted to 1/f were similar between groups for fre- quencies > 150 Hz. HFO events were defined from the bandpass filtered signal (80-120 Hz) by thresholding (> 4 SDs). The power spectra were evaluated in a window of ± 0.2 ms around each detected event. Time-frequency analysis was performed by applying the multitaper spectral estimation in sliding windows with 97.7% overlap and frequency resolution of 10 Hz in the 90-600 Hz fre- quency range. HFO activity was evaluated as the power integral in the 80-120 Hz band. We used a combination of features to identify HFO events automatically and to classify them in different categories (Valero et al., 2017) (Figure S1). First, we identified large amplitude transient (< 100 ms) discharges using LFP and current-source-density signals (CSD, i.e., the second spatial derivative) at the SR and SLM. Second, we identified HFO events at the SP by frequency thresholding over 100 Hz. Then we used spectral indices such as entropy and fast ripple indices from candidate HFO events at SP together with amplitude information from LFP and CSD signals at SR to classify events as: a) SPW-ripples (low amplitude; 100-150 Hz); b) SPW-fast ripples (medium amplitude, > 150 Hz) and IID-HFO (larger amplitude, > 100 Hz). Events not meeting criteria were left unclassified. Surface EEG recordings from mice were analyzed similarly, by using channels over the dorsal hippocampus or at the frontal cortex to identify HFO events and channels over the parietal cortex using selected segments of the EEG. In vivo recording and labeling of single cells Rats were anesthetized with urethane (1.2 g/kg, i.p.), fastened to the stereotaxic frame and warmed to keep their body temperature at 37 C. Two bilateral craniotomies of 1 mm diameter were performed for CA3 stimulation (AP: 1.2 mm, ML: 2.9 mm) and CA1 re- cordings (AP: 3.7 mm; ML: 3 mm). The dura was gently removed, the cisterna magna was drained and the craniotomy covered with warm agar to reinforce stability. A 16-channel silicon probes (NeuroNexus Tech; 100 mm interspaced, 413 mm2 contact) was advanced perpendicular along the CA1-DG-CA3c axis guided by extracellular stimulation and electrophysiological hallmarks. Concentric bipolar electrodes were advanced 3.5 mm with 30 in the coronal plane to stimulate CA3. Stimulation consisted of biphasic square pulses (0.2 ms duration, 0.05-1.2 mA every 5 s). Extracellular signals were pre-amplified (4x gain) and recorded with a 16(32)-channel AC amplifier (Multi- channel Systems), further amplified by 100, analogically filtered at 1Hz to 5 kHz, and sampled at 20 kHz/channel with 12 bits precision with a Digidata 1440. A subcutaneous Ag/AgCl wire in the neck served as reference. Recording and stimulus position was confirmed by post hoc histological analysis. Intracellular recording and labeling were obtained in current-clamp mode using sharp pipettes (1.5 mm/0.86 mm outer/inner diam- eter borosilicate glass; A-M Systems, Inc) filled with 1.5 M potassium acetate and 2% Neurobiotin (Vector Labs, Inc; 50-100 MU) (Cid and de la Prida, 2019). Signals were acquired with an intracellular amplifier (Axoclamp 900A) at 100x gain. The resting potential, input resistance and amplitude of action potentials was monitored all over the course of experiments. After data collection, Neurobiotin was ejected using 500 ms depolarizing pulses at 1-3 nA at 1 Hz for 10-45 min. Animals were perfused with 4% paraformaldehyde (PFA) and 15% saturated picric acid in 0.1 M, pH 7.4 phosphate buffered saline (PBS). Brains were postfixed overnight at room temperature (RT), washed in PBS and serially cut in 70 mm coronal sections (Leica VT 1000S vibratome). Sections containing the stimulus and probe tracks were identified with a stereomicroscope (S8APO, Leica). Sections containing Neurobiotin-labeled cells were localized by incubating them in 1:400 Alexa Fluor 488-conjugated streptavidin (Jackson ImmunoResearch 016-540-084) with 1% Triton X-100 in PBS (PBS-Tx) for 2 hours at room temperature (RT). To evaluate morpho- logical features of single cells, sections containing the somata of recorded cells were processed with Triton 1% in PBS, blocked with 10% fetal bovine serum (FBS) in PBS-Tx and incubated overnight at RT with the primary antibody solution containing rabbit anti-Cal- bindin (1:500, CB D-28k, Swant CB-38) or mouse anti-Calbindin (1:1000, CB D-28k, Swant 300) antibodies with 1% FBS in PBS-Tx. After three washes in PBS-Tx, sections were incubated for 2 hours at RT with appropriate secondary antibodies: goat anti-rabbit Alexa Fluor 633 (1:500, Invitrogen, A21070), and goat anti-mouse Alexa Fluor488 (Jackson Immunoresearch, 115-545-003) or goat anti-mouse Rhodamine Red (1:200, Jackson ImmunoResearch, 115-295-003) in PBS-Tx-1% FBS. Following 10 min incubation with bisbenzimide H33258 (1:10000 in PBS, Sigma, B2883) for nuclei labeling, sections were washed and mounted on glass slides in Mowiol (17% polyvinyl alcohol 4-88, 33% glycerin and 2% thimerosal in PBS). All morphological analyses were performed blindly to electrophysiological data. The distance from the cell soma to radiatum was measured from confocal images using information from Calbindin and bisbenzimide staining. All pyramidal cells included in this study were localized within the CA1 region. Calbindin immunostaining was used to estimate the width of the superficial sub-layer from the border to the stratum radiatum. Superficial cells were defined based on the location of the soma within the calbindin sublayer, inde- pendently on their immunoreactivity (Valero et al., 2015). The border with radiatum was estimated for each section and the distance from the recorded cell somata was measured using ImageJ (NIH Image). Analysis of intracellular single-cell recordings Interictal discharges and sharp-waves associated with HFO events (ripples or fast ripples) recorded simultaneously to the intracel- lular membrane potential were detected as explained before. The timing of these events was used to identify the corresponding membrane potential deflections. Time-frequency analysis of HFO events was performed by applying the multitaper spectral estima- tion in sliding windows with 97.7% overlap and frequency resolution of 10 Hz in the 90-600 Hz frequency range (only the 100-600 Hz range is shown) to data sweeps aligned by sharp-wave ripple events (±1 s). Membrane potential responses of single-cells were eval- uated in peri-event plots before (200 to 150 ms), during (±50 ms) and after (150 to 200 ms) HFO events. e4 Cell Reports 35, 109229, June 8, 2021 ll Article OPEN ACCESS Passive electrophysiological properties (input resistance, membrane decay and capacitance) of neurons recorded intracellularly in vivo were measured using 500 ms currents step in current-clamp mode. Cells with intracellular action potential amplitude smaller than 40 mV were excluded. Resting membrane potential and input resistance were estimated by linear regression between baseline potential data and the associated holding current. Intrinsic firing properties, including action potential threshold, half-width duration and AHP were estimated from the first spike in response to depolarizing current pulses of 0.2 nA amplitude and 500 ms duration. The sag and maximal firing rate were calculated from current pulses of ± 0.3 nA amplitude. A bursting index was defined as the ratio of the number of complex spikes (minimum of 3 spikes < 8 ms inter-spike interval) over the total number of spikes recorded during theta activity. cFos immunostaining and analysis To evaluate immediate-early gene expression associated to sound-induced convulsive seizures, animals were perfused 1 hour after and their brains cut in 70 mm coronal sections. Selected sections were stained against c-Fos using a polyclonal antibody at 1:250 (Santa Cruz Biotechnology sC-52) and bisbenzimide. Using one 20x confocal mosaic per animal, we quantified cFos intensity at CA1 pyramidal cells by delineating single-cell nuclei stained with bisbenzimide in one confocal plane (ImageJ). The mean intensity of cFos signal from each cell was then normalized by subtracting the background (set at 0) and dividing by the maximal positive signal in the mosaic, which was always at granule cell level (set at 1). No significant differences of background were observed across sec- tions. Delineated cells were ranked by their distance to radiatum to classify them as deep or superficial, according to standard mea- surements of Calb1-layer thickness. Laser capture microdissection (LCM) and RNA isolation Brains from 3 control and 3 epileptic rats were dissected, longitudinally cut in half (separating both hemispheres), wrapped in aluminum foil and immediately frozen in liquid nitrogen. To avoid circadian effects on gene expression all samples were collected in the morning before noon and conserved at 80 C until use. The hippocampal region of each hemisphere was cut in 20 mm slices in a cryostat (Leica) (chamber temperature: 20 C; block temperature: 30 C) and placed on 1.0 mm PEN-membrane covered slides (Carl Zeiss). Slides were conserved at 20 C until use. Right before microdissection, slides were dried with vapor of liquid ni- trogen. The CA1 cell layer was microdissected with a Leica 6000 laser microdissector through a 40x non-oil immersion objective to obtain cell bodies of superficial and deep sublayers separately (Figure 2A). Microdissected deep and superficial areas were collected in different empty caps of 0.5 mL Eppendorf tubes. After microdissection, samples were processed following ARCTURUS PicoPure RNA Isolation Kit (Thermo Fisher Scientific) instructions in order to extract and isolate total RNA. Briefly, 50 ml of extraction buffer was added into the cap, incubated at 42 C for 30 min, centrifuged at 800 x G for 2 min and stored at 80 C. The same volume of 70% ethanol was added to the cell extract and the mixture was pipetted into a pre-conditioned RNA purification column. The column was centrifuged 2 min at 100 x G and 30 s at 16000 x G, and washed with 100 ml of Wash Buffer 1. To completely eliminate DNA, the pu- rification column was treated with 40 ml of DNase (diluted 1/8 in RDD Buffer) (QIAGEN), incubated 15 min and centrifuged at 8000 x G for 15 s. Then the column was washed twice with 100 ml of Wash Buffer 2, and centrifuged at 8000 x G after the first wash and at 16000 x G after the second one. Finally, RNA was eluted into a new 0.5 ml Eppendorf tube by adding 11 ml of elution buffer onto the column membrane, incubating the column for 1 min at room temperature, and centrifuging the column for 1 min at 1000 x G and at 16000 x G immediately after. Total RNA samples were stored at 80 C. RNA integrity number (RIN) was similar in control (4.7 ± 0.8) and epileptic rats (5.7 ± 0.5; p = 0.07; 3 replicates x 2 sublayers per group), as well as for deep (5.2 ± 0.9) and superficial sublayers (5.3 ± 0.1; p = 0.81; n = 6 replicates per sublayers). LCM RNA-seq library construction and sequencing RNA preparation for sequencing deep and superficial CA1 sublayers from control and epileptic rats was performed as described in Scandaglia et al. (2017). The twelve samples were sequenced according to manufacturer instructions in a HiSeq2500 sequencer (Illumina, Inc). Libraries were strand specific (reverse) and sequenciation was performed in paired-end configuration with a read length of 75 bp. Library size of read pairs for the different samples analyzed was between 47 and 59 Million reads. RNA-seq data can be accessed at the GEO repository (GSE143555). LCM RNA-seq data analysis Alignment quality control of sequenced samples (LCM RNA-seq) was assessed with FastQC (v.0.11.3) (Babraham Institute) and RNA-seq tracks were visualized using IGV (v.2.3.57) (Thorvaldsdóttir et al., 2013). LCM RNA-seq reads were mapped to the rat genome (Rnor_6.0.83) using STAR (v.2.5.0c) (Dobin et al., 2013), and files were further processed with Samtools (v.0.1.19). Aligned reads were counted to gene transcripts using HTSeq (v.0.6.1) (Anders et al., 2015). Differential expression analysis was performed using DESeq2 (v.1.10.0) (Love et al., 2014) of the bioconductor suite (Huber et al., 2015) in the R (v.3.2.2) statistical computing platform. The experimental design consisted in two factors (treatment and anatomical area) and there was also grouped samples (samples from different anatomical areas (minus (deep) and plus (superficial) that were obtained from individual mice). Genes were considered differentially expressed at Benjamini-Hochberg (BH) adj. p < 0.05 and FC > 0.3 (Benjamini and Hochberg, 1995), except otherwise specified. GO analysis was performed using DAVID (v.6.8) bioinformatics platform (Huang et al., 2009). Cell Reports 35, 109229, June 8, 2021 e5 ll OPEN ACCESS Article LCM RNA-seq single-cell informed data analysis Single-cell RNA-seq data from Zeisel et al. (2015) was reanalyzed with consensus cluster SC3 algorithm (Kiselev et al., 2017). From the original 3,005 cells, pyramidal and interneurons were removed from somatosensory cortex, resulting in a total of 2,442 cells. Re- maining cells were re-analyzed downstream with SC3. Clustering stability was optimal for 6 clusters. One cluster presented high het- erogeneity (mixed population cluster), and was reanalyzed with the optimal clustering stability (5 clusters). Marker genes were tested for every cluster with Wilcoxon signed rank test. Top 10 genes with the area under the receiver operating characteristic (ROC) curve (AUC) power > 0.85 and with the adj. p < 0.01 from both cluster analyses were selected, with a total of 69 bona fide population markers. Based on the markers, populations were fused/splitted and 3 populations were isolated in the first clustering (pyramidal neurons, interneurons and oligodendrocytes) and 4 more in the second (astrocytes, endothelial cells, microglia and mural cells). Sixteen outlier cells were removed by total_counts, total_features or pct_counts_spike criteria, remaining 2,426 cells with a high cor- respondence with the original classification: Astrocytes (155 cells), Endothelial (177 cells), Interneurons (174 cells), Microglia (85 cells), Mural (56 cells), Oligodendrocytes (804 cells) and Pyramidal (975 cells). Next, to capture convoluted gene signatures on LCM RNA-seq we used gene sets of equal size that were differentially regulated between superficial and deep sublayers of control and epileptic rats. To obtain these gene sets, we first filtered the genes by significance (adj. p < 0.1), then significantly differentially regulated genes were ranked by fold-change, and top 250 differentially regulated genes for each sublayer and condition, were re- tained for further analysis. We then, obtained normalized gene expression from scRNA-seq data for these gene sets and performed linear dimensionality reduction by principal component analysis to summarize and visualize the cells in the low-dimensional space. We also performed Pearson pairwise correlation and hierarchical clustering for these gene sets across all cells, to capture cell-type specific gene signatures that were mapped to major cell types in CA1 by using previously identified cell markers. Single-cell RNA-seq data from Allen Brain Map portal - Mouse Whole Cortex and Hippocampus SMART-seq (2019) with 10x-SMART- seq taxonomy (2020) (http://portal.brain-map.org/atlases-and-data/rnaseq/mouse-whole-cortex-and-hippocampus-smart-seq) was downloaded and subsetted as follows: from the original 76,533 cells, cells from hippocampal region were retained, and subsequently, cells from CA2, CA3, DG or empty subclass were removed. Then, cells of subclasses Astro, Micro-PVM, Endo, Oligo or VLMC were added. This initial filtering led a total of 5,506 cells. Of these, cells tagged by subclass as CA1-ProS (1,592 cells), Astro (976 cells), Lamp5 (864 cells), Vip (462 cells), Sncg (416 cells), Sst (266 cells), Oligo (236 cells), Endo (213 cells), Micro-PVM (176 cells), VLMC (159 cells), Pvalb (69 cells) were retained for further analysis, while cells tagged by subclass as SUB-ProS (29 cells), L2 IT RHP (21), NP SUB (9), L2/3 IT CTX-1 (3), L5 IT TPE-ENT (3), Sst Chodl (3), L2/3 IT ENTl (2), L6 CT CTX (2), Meis2 (2), L2/3 IT CTX-2 (1), L2/3 IT PPP (1), L6b CTX (1) were removed because of their low cell number and/or mismatch with CA1 or glial cells. The remaining 5,429 cells were grouped in their corresponding major cell populations as follows: Astrocytes (Astro) (976 cells), Endothelial (Endo) (213 cells), In- terneurons (Lamp5, Pvalb, Sncg, Sst, Vip) (2077 cells), Microglia (Micro-PVM) (176 cells), Mural (VLMC) (159 cells), Oligodendrocytes (Oligo) (236 cells) and Pyramidal (CA1-ProS) (1,592 cells). Next, we proceed as indicated before to capture convoluted gene signatures on LCM RNA-seq by using normalized expression data from this scRNA-seq dataset, for differentially regulated between superficial and deep sublayers of control and epileptic rats. Single-nuclei isolation We accurately isolated single nuclei from 2 control and 2 epileptic Calb1::CrexTdTomato young adult mice. We used these animals to facilitate identification of the CA1 region under a fluorescent scope. Animals were sacrificed 12 weeks after the kainate/saline admin- istration by cervical dislocation and brains were dissected and cut in 300 um thick slices in a vibratome covered by ice cold HBBS 1x. As in LCM studies, samples were collected before noon to avoid circadian effects. The dorsal CA1 region were manually dissected from 4 consecutive slices and put altogether into 400 mL of ice-cold MACS buffer (0,5% BSA, 2 mM EDTA, PBS 1x). CA1 portions were transferred to a dounce homogenizer (20404 Lab Unlimited) containing 400 mL of MACS buffer and were homogenized 12-15 times each with the pestle. The cell suspension was transferred to a 2 mL Eppendorf tube and centrifuge 15 min 500 x G 4 C. Cell pellets were resuspended in 500 mL of lysis buffer (10 mM Tris-HCL, 10 mM NaCl, 3 mM MgCl2, 0,1% IGEPAL) and kept 5 min on ice. Samples were then spun down at 500 x G for 30 min in a pre-chilled centrifuge. The nuclei pellet was resuspended in PBS 1x 1% BSA and then 15,000 nuclei were purified by flow cytometry in a BD FACS Aria III. The whole process was carried out at 4 C. Single-nucleus RNA sequencing Purified intact nuclei from mouse hippocampal CA1 area were processed through all steps to generate stable cDNA libraries. For every sample, 15,000 nuclei were loaded into a Chromium Single Cell A Chip (10x Genomics) and processed following the manufac- turer’s instructions. Single-nuclei RNA-seq libraries were prepared using the Chromium Single Cell 30 Library & Gel Bead kit v2 and i7 Multiplex kit (10x Genomics). Pooled libraries were then loaded on a HiSeq2500 instrument (Illumina) and sequenced to obtain 75 bp paired-end reads following manufacturer instructions. On sequencing depth, 262 million fragments were generated for the control condition and 296 for the epileptic dataset. Libraries reached a sequencing saturation of 86.9% for control and 91.2% for epilepsy condition. snRNA-seq data can be accessed at the GEO repository (GSE143560). snRNA-seq analysis Quality control of sequenced reads was performed using FastQC (Babraham Institute). Sequenced samples were processed using the Cell Ranger (v.2.2.0) pipeline (10x Genomics) and aligned to the GRCm38 (mm10) mouse reference genome customized to count e6 Cell Reports 35, 109229, June 8, 2021 ll Article OPEN ACCESS reads in introns (pre-mRNA) over the Ensembl gene annotation version 94. Barcodes with total unique molecular identifier (UMI) count > 10% of the 99th percentile of the expected recovered cells were selected for further analysis. Using this criterion, we retrieved 3,661 (control), 3,078 (epileptic) high quality nuclei per sample. Mean reads per nucleus were 71,449 (control) and 96,150 (epileptic). Median genes per nucleus were 71,449 (control) and 96,150 (epileptic). Minimum UMI count per nucleus were 710 (control), 540 (epileptic), well above the typical quality standards in single-cell/nucleus sequencing. Single-nucleus RNA-seq data were subsequently pre-pro- cessed and further analyzed in R (v.3.4.4) using Seurat (v.2.3.4) (Butler et al., 2018; Stuart et al., 2019). Filtering parameters were as follows: genes, nCell < 5; cells, nGene < 200. Data were then normalized using global-scaling normalization (method: LogNormalize, scale.factor = 10,000). An initial exploratory analysis was performed on each dataset separately. This analysis retrieved a similar number of clusters in each dataset that were approximately equal in size. We next combined both datasets using the function MergeSeurat. Highly variable genes (HVGs) were detected using FindVariableGenes function with default parameters. Then, normalized counts on HVGs were scaled and centered using ScaleData function with default parameters. Principal component analysis (PCA) was performed over the first ranked 1,000 HVGs. Plots of the two principal components of the PCA, where cells were colored by dataset, excluded the presence of batch effects. Cluster detection was carried out with Louvain algorithm in FindClusters function, using 20 first PCA dimensions and resolution of 0.6 (the default number in Seurat and the optimal according to cell number, data dispersion and co-expression of previously reported cell markers). Plots of the two principal components of the PCA where cells were colored by dataset of origin excluded the presence of batch effects. This analysis identified 13 clusters. CA1 pyramidal neurons populated three of these clusters: one cluster was enriched in bona fide gene markers of deep cells (Ndst4, Coll11a1) whereas a second one was enriched in canonical markers of superficial cells (Calb1, Epha3). The third cluster showed a mixed identity. An additional round of clustering segmented this population in three additional clusters that were enriched in deep and superficial markers, respectively, and a third cluster that could not be annotated based on the presence of gene markers of CA1 sublayer neurons. The vast majority of cells within this cluster (66/67) were from epileptic mice (Pyr_ES). Next, the FindMarkers function was used to identify gene markers and to determine the cell populations represented by each cluster. Finally, cell subtypes were manually aggregated based on the presence of canonical markers of known cell types into six distinct major cell types: excitatory neurons (Excit); inhibitory neurons (Inter); oligodendrocytes (ODC); oligodendrocyte progenitor cells (OPCs); microglia (Microglia); astrocytes (Astro). Visualization and embedding were performed using stochastic nearest neighbors (tSNE) (van der Maaten and Hinton, 2009) and uniform manifold approximation and projection (UMAP) (McInnes et al., 2018) methods over PCA using the 20 first PCA dimensions. UMAP plots of gene expression show normalized count (UMIs) per nucleus. The equalized expression between fixed percentiles was plotted according to the following criterion: the minimum expression was adjusted to 5% and the maximum expression was adjusted to 95% in all UMAP expression plots. To evaluate effects of epilepsy, datasets from both conditions were merged and HVGs were identified for each dataset as above indicated. Only HVGs that were detected in all datasets were used to perform visualization and embedding as described above. Clustering was performed on merged dataset from both conditions and populations were identified combining these results with clustering information obtained in control and epileptic datasets separately, together with co-expres- sion of population markers. Differential expression analysis (DEA) was used to identify population gene markers. For DEA, the nuclei of each population were contrasted against all the other nuclei in the merged dataset using Wilcoxon Rank Sum test on normalized counts. For epilepsy effect analysis, in the merged dataset, the nuclei of each population from the epileptic dataset were contrasted against all the other nuclei of the same population in control using Wilcoxon Rank Sum test on normalized counts. GO functional enrichment analyses were per- formed using DAVID (v.6.8) bioinformatics platform (Huang et al., 2009). All DEA and functional enrichment analyses applies False Discovery Rate (FDR) method by BH post hoc p adjustment, except otherwise specified. Cell trajectories and pseudotime analysis The disease pseudotime analysis was performed using Monocle 2 (v.2.8.0) (Qiu et al., 2017a, 2017b). First, the Seurat merged dataset was transformed to Monocle object and cells from Pyr_CA1 were subset. The size factor and dispersion of the subset was estimated, and data was normalized and preprocessed. Genes under the minimum level detection threshold of 0.1 and detected in less than 10 cells were filtered with the function setOrderingFilter. Genes defining how a cell progress through a pseudo-time disease trajectory were selected with the function differentialGeneTest (Monocle’s main differential analysis routine). 2,579 genes (64,68% of a total of 3,987 genes considered as expressed) were significative with FDR adj. p < 1% for the combination of factors: SeuratCluster+Condition, and thus, defined the high dimensional space for pseudotemporal trajectory analysis. Discriminative dimensionality reduction with trees (DDRTree) reduction algorithm learns the principal graph and specifies the tra- jectory. DDRTree was applied inside the function reduceDimension, and got the default parameters: norm_method = ‘‘log,’’ pseu- do_expr = 1, relative_expr = TRUE, auto_param_selection = TRUE (automatically calculate the proper value for the ncenter (number of centroids)) and scaling = TRUE (scale each gene before running trajectory reconstruction). Prior the dimensional reduction, the function reduceDimension also performed a variance-stabilization of the data (because the expressionFamily of the data was neg- binomial.size). Finally, the cells were ordered according to pseudo-time with the function orderCells, which added a pseudo-time value and state for each cell; together encode where each cell maps to the trajectory. Cell Reports 35, 109229, June 8, 2021 e7 ll OPEN ACCESS Article For enrichment analysis on pseudotime trajectories, top-ranked 250 branched expression analysis modeling (BEAM) significant changes through the progression across the disease trajectory for each sublayer were clusterized and GO enrichment analyses on upregulated and downregulated clusters were performed using DAVID (v.6.8) bioinformatics platform (Huang et al., 2009). Cell-type immunostaining and analysis To evaluate the contribution of different cell-types to hippocampal sclerosis, control and epileptic rats and mice were perfused with 4% paraformaldehyde (PFA) and 15% saturated picric acid in 0.1 M PBS, pH 7.4. Brains were postfixed overnight and cut in 70 mm coronal sections (Leica VT 1000S vibratome). Sections containing the dorsal-intermediate hippocampus were processed with Triton 1% in PBS and blocked with 10% fetal bovine serum (FBS) in PBS-Tx. Sections were incubated overnight at RT with 1% FBS PBS-Tx solution containing primary antibodies against a battery of cell-type specific markers. The list of antibodies includes: rabbit anti-cal- bindin (1:1000, CB D-28k, Swant CB-38) or mouse anti-calbindin (1:500, CB D-28k, Swant 300) to identify superficial CA1 pyramidal cells; rabbit anti-Wfs1 (1:500, Protein Tech, 11558-1-AP) for CA1 pyramidal cells; rabbit anti-Iba1 (1:1000, Wako, 019-19741) for mi- croglia; rabbit anti-GFAP (1:1000, Sigma, G9269) for astrocytes; rabbit anti-Olig2 (1:200, Millipore, AB9610) for oligodendrocytes. After three washes in PBS-Tx, sections were incubated for 2 hours at RT with appropriate secondary antibodies: goat anti-rabbit Alexa Fluor633 (1:200, ThermoFisher, A-21070), and donkey anti-mouse Alexa Fluor488 (1:200, ThermoFisher, A-21202) or goat anti-mouse Rhodamine Red (1:200, Jackson ImmunoResearch, 115-295-003) in PBS-Tx-1%FBS. Following 10 min incubation with bisbenzimide H33258 (1:10000 in PBS, Sigma, B2883) for nuclei labeling, sections were washed and mounted on glass slides in Mowiol (17% polyvinyl alcohol 4-88, 33% glycerin and 2% thimerosal in PBS). To acquire multichannel fluorescence stacks from recorded cells, a confocal microscope (Leica SP5) with LAS AF software v2.6.0 build 7266 (Leica) was used. The following channels (fluorophore, laser and excitation wavelength, emission spectral filter) were used: a) bisbenzimide, Diode 405 nm, 415–485 nm; b) Alexa Fluor 488 nm, Argon 488 nm, 499–535 nm; c) Rhodamine Red / Alexa Fluor 568 / Texas Red, DPSS 561 nm, 571–620 nm; d) Alexa Fluor 633, HeNe 633 nm, 652–738 nm; and objectives HC PL APO CS 10.0 3 0.40 DRY UV, HCX PL APO lambda blue 20.0 3 0.70 IMM UV and HCX PL APO CS 40.0 3 1.25 OIL UV were used. FluoroJade staining To evaluate neurodegenerating neurons we used coronal sections from epileptic rats perfused at different time points post-status. Selected sections were immunostained against Wfs1 followed by FluoroJade staining. To this purpose, sections were pretreated for 5 min with 1% sodium hydroxide in 80% ethanol, followed by 70% ethanol (2 min) and distilled water (2 min). Sections were then incubated 10 min in 0.06% potassium permanganate, rinsed in distilled water and immersed into 0.0001% solution of FluoroJade C dye (Sigma AG325) in 0.1% acetic acid (pH 3.5) for 10 min. After a brief wash in distilled water, they were mounted on gelatin-coated slides, air-dried, coverslipped with DPX and examined under a confocal microscope as described above. FluoroJade positive cells exhibited bright green fluorescence. In situ hybridization analysis Selected sections from control and epileptic rats were processed for in situ hybridization using standard methods. Briefly, riboprobes were prepared from Rat Enpp2 cDNA (Image clone ID 7115236) and using RT-PCR from rat adult hippocampus to prepare Wfs1 (NM_031823.1, from bp 783 to 1631), Ndst4 (XM_006233274.2, from bp 536 to 1024), Syt17 (NM_138849.1, from bp 378 to 1118), Hrt1a (J05276.1, from bp 730 to 1663) and Scn7a (NM_131912.1, from bp 2812 to 3550) cDNAs. Similar riboprobes were used in sections from normal and epileptic mice. Riboprobe hybridization was detected using alkaline phosphatase–coupled anti- digoxigenin Fab fragments (Sigma). Hybridized sections were mounted in glycerol and photographed using a Nikon stereoscope and a DCC Nikon camera. RNAscope analysis Control and epileptic mice were perfused with 4% paraformaldehyde (PFA) and 15% saturated picric acid in 0.1 M PBS, pH 7.4. Brains were post-fixed overnight and cut in 50 mm coronal sections (Leica VT 1000S vibratome) which were kept at 20 C in a so- lution of 30% glycerol and 30% etylenglycol in PBS. RNAscope Fluorescent Multiplex Assay (Advanced Cell Diagnostics) was carried out according to the manufacturer’s protocols. Briefly, sections containing the dorsal-intermediate hippocampus (around 2 mm AP from Bregma) were washed in PBS-Tx, and mounted onto SuperFrost Plus microscope slides (10149870, ThermoFisher Scientific). They were then dehydrated at 60 C follow by ethanol, pretreated with a target retrieval solution (322000, ACD) and protease III (322340, ACD), and co-hybridized with Spag5 (Mm-Spag5, 505691, ACD) and Dcc (Mm-Dcc-C3, 427491, ACD) or Dapk1 (Mm- Dapk1-C3, 880221-C3, ACD) probes. Finally, the amplification steps (RNAscope Fluorescent Multiplex Detection reagents, 320851, ACD) were followed, using Atto 550 for Spag5 and Atto 647 for Dcc or Dapk1 as fluorescent labels. The RNAscope 3- plex positive control probe set (320881, ACD), with probes to Polr2a, PPIB and UBC, was used to confirm preservation of sample RNA. The negative control probe to bacterial DapB (320871) was used to establish non-specific labeling. Following the RNAscope protocol, sections were blocked for 30 min with 10% FBS in PBS-Tx, and incubated overnight at RT with rabbit anti-Wfs1 (1:500, Protein Tech, 11558-1-AP) in 1% FBS PBS-Tx. After three washes in PBS, sections were incubated for 2 h with donkey anti-rabbit Alexa Fluor488 (1:200, ThermoFisher Scientific, A-21206) in 1% FBS PBS-Tx, washed twice in PBS and mounted using ProLong Gold Antifade mountant (ThermoFisher Scientific, P10144). e8 Cell Reports 35, 109229, June 8, 2021 ll Article OPEN ACCESS Multichannel fluorescence stacks were achieved in a confocal microscope (Leica SP5) with LAS AF software v2.6.0 build 7266 (Leica), with a 40x objective (HCX PL APO CS 40.0 3 1.25 OIL UV), at 0.5 mm z-interval. Pinhole was set at 1 Airy, and following chan- nels settings were applied for each fluorophore (excitation wavelength, emission spectral filter): Argon (488 nm, 499–553 nm) for Alexa Fluor488, DPSS (561 nm, 570–630 nm) for Atto 550, and HeNe (633 nm, 645–740 nm) for Atto 647. Diode (405 nm, 415– 485 nm) was used as an unstained channel to identify autofluorescent material, which can be abundant in epileptic tissue. Images were acquired at 8 bits, 387.5 3 387.5 mm (1024 3 1024 pixels) in size. Brightness and contrast were adjusted with the ImageJ Fiji software (NIH Image). To estimate the amount of Dapk1, Spag5 and Dcc transcripts per cell, a total of 5 control and 5 epileptic animals (12 cells per an- imal) were analyzed. Each ROI was drawn in one confocal plane as the outline of a pyramidal cell soma at their maximum diameter based on Wfs1 immunoreactivity, using ImageJ. Only those cells focused within the first 10 mm were selected to avoid reduction of the signal due to depth in the section. In either Dapk1, Spag5 or Dcc channels, single signal dots were counted, and clusters were converted to dot number as suggested by the manufacturer (ACD): (total intensity - average background intensity x area)/single dot average intensity. QUANTIFICATION AND STATISTICAL ANALYSIS Statistical analyses were performed with MATLAB, SPPSS and using the computing environment R (R Development Core Team, 2005, http://www.R-project.org). No statistical method was used to predetermine sample sizes. Normality and homoscedasticity were evaluated with the Kolmogorov–Smirnov and Levene’s tests, respectively. The exact number of replications for each experi- ment is detailed in text and figures. Several ways ANOVAs or Friedman tests were applied. Post hoc comparisons were evaluated with the Tukey-Kramer, Student or Wilcoxon tests. Deep-superficial trends were evaluated using Spearman correlation and tested against 0 (i.e., no correlation was the null hypothesis) at p < 0.05 (two sided). Cell Reports 35, 109229, June 8, 2021 e9