1 Introduction

The neural operations supporting memory require the involvement of various brain systems interconnected through neural networks (Kitamura et al., 2017). The activity of these networks is fine-tuned by experience through the selective recruitment of ensembles of neurons ascribed to contain the bits of information associated with memory (Josselyn and Tonegawa, 2020). Neuronal ensemble activity is shaped by synaptic plasticity mechanisms that modulate the weight and efficacy of the ensemble’s synaptic connections (Takeuchi et al., 2014; Malenka and Bear, 2004). Changes in gene expression are a key underlying mechanism in synaptic plasticity (Mayford et al., 2012), and recent studies have identified multiple profiles of gene expression associated with memory (Marco et al., 2020; Rao-Ruiz et al., 2019; Reshetnikov et al., 2020).

In the rodent hippocampus, encoding of spatial memory information is supported by the diversity of synaptic computations within and across its sub-fields (Poo et al., 2016; Mayford, 2014). Each hippocampal subregion (i.e., Dentate Gyrus, CA1, CA2, CA3) shows distinct patterns of activity to process spatial information and to form or retrieve memories (Middleton and McHugh, 2020; Chevaleyre and Siegelbaum, 2010; Senzai, 2019; Basu and Siegelbaum, 2015; Amaral and Witter, 1989; Sloviter and Lomo, 2012). At the regional level, patterns of activity are conferred through the combination of cyto-architecture, local synaptic circuitry, neuronal cell types and synaptic inputs unique to each region (van Strien et al., 2009; Senzai, 2019; Kesner et al., 2004; Lisman, 1999; Ryan et al., 2015; Hunsaker and Kesner, 2013; Molitor et al., 2021; Rao-Ruiz et al., 2019; Shpokayte et al., 2022; Sullivan et al., 2021; Shah et al., 2016; Cho et al., 2016). At the single-cell level, the functional state of hippocampal neurons is defined by connectivity, hippocampal subregional location, and experience-dependent recruitment into a memory ensemble (Josselyn and Tonegawa, 2020; Tonegawa et al., 2018). The investigation of memory across molecular and functional levels necessitates an approach that connects the scales of hippocampal organization between cellular mechanisms and networks of cells. Though hundreds of individual genes supporting synaptic plasticity at the cellular level have been identified (Sanes and Lichtman, 1999; Alberini and Kandel, 2014; Asok et al., 2019; Magee and Grienberger, 2020; Malenka and Bear, 2004), the transcriptomic profiles which characterize memory across hippocampal subregions and memory-associated neuronal ensembles are understudied.

In this study, we present an exploratory approach to expand our understanding of the spatial organization of transcriptomic profiles by providing evidence that hippocampal subregions undergo distinct differential gene expression of biological mechanisms to support memory. We characterized changes in gene expression across the major subregions of the hippocampus using bulk RNA sequencing of micro-dissected dorsal hippocampal subregions from mice that learned an active place avoidance memory task. The regionalized changes in gene expression with memory in the hippocampus is a developing area of investigation often using behavioral paradigms like fear conditioning or the Morris Water Maze (Cho et al., 2016; Sardoo et al., 2022; Reshetnikov et al., 2020; Harris et al., 2020). Our approach involved training mice in the active place avoidance paradigm (Cimadevilla et al., 2000; Stuchlik et al., 2013), a spatial memory task that distinguishes itself from these other aversive paradigms because it requires proper utilization of cognitive control to solve the spatial task (Chung et al., 2021; Kelemen and Fenton, 2016). Results from our bulk RNA sequencing data suggest that the consolidation of an active place avoidance memory is supported by regionalized transcriptomic profiles separating the CA1 and CA3 from the DG.

Following the establishment of the memory-associated transcriptomic profiles of bulk hippocampal subregions we deepened our exploration of these transcriptomic changes using an innovative RNA sequencing methodology known as spatial transcriptomics. Spatial transcriptomics provides localized whole transcriptome sequencing across tissue sections at near single cell resolution (a.k.a. spatial transcriptomic spots) (Stahl et al., 2016). With this approach, we investigated the putative memory-associated neuronal ensemble by characterizing the spatial distribution of immediate early gene (IEG)-associated gene expression across coronal dorsal hippocampal sections. Memory-associated neuronal ensembles are a sparsely distributed population of neurons which are strongly activated together during the acquisition and recall of a memory, and are thought to encode the traces of information relevant to that memory (Ryan et al., 2015; Guan et al., 2016; Chen et al., 2019). Numerous studies have identified IEGs such as Arc, c-Fos, Egr1, c-Jun, Npas4, Bdnf, and Fmrp to characterize memory-associated neuronal ensembles (Sauvage et al., 2019; Josselyn and Tonegawa, 2020; Meenakshi et al., 2021; Erwin et al., 2020). Moreover, the expression of Arc mRNA has been used to tag and characterize the properties of memory-associated neuronal ensembles (Guzowski et al., 1999; Ramirez-Amaya et al., 2005; Lee et al., 2022).

Through our combined understanding of IEG activation and memory-associated neuronal ensembles, we used the detection of IEG expression in spatial transcriptomic spots to infer the location of these ensembles and explore the transcriptomic changes supporting spatial memory in the hippocampus. Consistent with our bulk RNA sequencing findings, spatial transcriptomics also showed differential gene expression between the hippocampal regions from an active place avoidance memory trained compared to an untrained mouse. We then investigated the difference between memory-associated neuronal populations through the analysis of spatial transcriptomic spots expressing the IEGs Arc, Egr1, and c-Jun. Comparison of the gene expression profiles from these IEG-expressing spatial transcriptomic spots between the memory trained and the untrained mouse revealed differential expression of genes involved in energy production and cytoplasmic translation. Novel analyses comparing IEG-expressing spatial transcriptomic spots within each mouse hippocampal sample (i.e., within the memory trained mouse or within the untrained mouse) revealed that IEG-expressing spatial spots were transcriptomically and spatially distinct from one another. With Arc, only the comparison between Arc-expressing and non-Arc-expressing spatial transcriptomic spots in the memory trained mouse –but not the untrained mouse– detected gene expression profiles linked with synaptic plasticity and memory. Additionally, similar comparisons amongst Egr1 and c-Jun expressing versus non-expressing spots revealed collections of differentially expressed distinguishing each of the investigated populations of IEG-expressing spots within the trained sample.

While still at an exploratory level, our spatial transcriptomics results highlight the role of Arc in shaping memory-associated ensembles and of the detection of gene expression profiles unique to the comparison Arc-expressing versus non-Arc-expressing spatial transcriptomic spots in an animal trained in the active place avoidance memory task. We speculate that a neuronal network being actively recruited in the consolidation of memory (such as the hippocampal network in our study) manifests transcriptomically divergent cellular microenvironments (such as the Arc-expressing and non-Arc-expressing spatial transcriptomic spots in our study) that could facilitate the processing of memory information.

2 Methods 2.1 Animals

A total of 18 adult ArcCreERT2::eYFPflx mice (Denny et al., 2014) with C57BL/6 genetic background aged 3–4 months were utilized across all experimental cohorts. Bulk RNA sequencing experiments utilized 3 male and 3 female mice, 2 females and 1 male were assigned to the trained behavioral condition, and 1 female and 2 males were assigned to the untrained behavioral condition (see section “Behavior” for details). Spatial transcriptomics experiments utilized 2 male mice, 1 mouse was assigned to each behavioral condition. RT-qPCR experiments utilized 12 male mice, with 6 mice assigned to each behavioral condition. Prior to experimental onset mice were bred in-house at the SUNY Downstate Health Sciences University vivarium (Brooklyn, NY, USA). Mice were housed in groups of two to five per cage. Beginning on the first day of training and testing, mice were single housed in shoebox cages in a sound attenuation cubicle (Med Associates). For every experimental run, a minimum of two mice —each undergoing the same behavioral conditioning— were simultaneously housed in the sound attenuation cubicle. Ad libitum food and water was provided. Mice were randomly assigned to behavioral cohorts before the start of the experiment. Mice were handled daily for 3 days prior to the start of the experiment to reduce anxiety and improve voluntary approach. All animal procedures proposed are approved by, and will be performed following, the Institutional Animal Care and Use Committee guidelines at SUNY Downstate Health Sciences University.

2.2 Behavior

All procedures were performed in compliance with the Institutional Animal Care and Use Committee of the State University of New York, Downstate Health Sciences University. ArcCreERT2::eYFPflx male mice were trained in a hippocampus-dependent two-frame active place avoidance task. The place avoidance system consisted of a 40-cm diameter arena with a parallel rod floor that could rotate at 1 rpm. The position of the animal was tracked using PC-based software (Tracker, Bio-Signal Group Corp., Brooklyn, NY) that analyzed 30-Hz digital video images from an overhead camera. Mice in the trained condition learned the active place avoidance task. Place avoidance of a 60° zone was reinforced by a constant current foot shock (60 Hz, 500 ms, 0.2 mA) that was scrambled (5-poles) across pairs of the floor rods. Rotation of the arena would carry the mouse into the shock zone unless the animal actively avoided the zone. Entering the shock zone for more than 500 ms triggered shock. Additional shocks occurred every 0.5 s until the animal left the shock zone. Measures of place avoidance were computed by TrackAnalysis software (Bio-Signal Group Corp., Brooklyn, NY, USA). The number of detections of the animal in the shock zone was tracked as the performance metric for the training trials.

On day one, mice assigned to the trained behavioral condition (trained mice) received a 30-min trial with the shock off to habituate to the rotating arena. Across the next 2 days the animals experienced four training trials, with two 30-min trials a day with the activated shock with a 40 min inter-trial-interval. Control (untrained mice) experienced identical training conditions, except for the shock always being off. The number of detections of the animal in the shock zone was tracked as the training performance metric for the training trials. Memory retention performance was assessed 24 h after the final training session in a 10-min retention test with the shock off. Latency to the animal’s second entrance into the shock zone was tracked as the memory retention performance metric for the retention test. Previous research indicates the latency to second entrance as a valid behavioral metric for retention test performance because it represents place avoidance memory recall that is less affected by to the animal’s potential errors in self-localization once placed in the rotating arena (Cimadevilla et al., 2000; Stuchlik et al., 2013; Cimadevilla et al., 2001).

2.3 Microdissected bulk RNA sequencing

60 min post retention test, mice were euthanized (5% isoflurane vaporized in 100% oxygen), and brains were extracted, washed in ice cold artificial cerebrospinal fluid, blocked, and mounted on a vibratome stage. Hippocampal subregions (DG, CA3, CA1) were microdissected from 400 μm thick coronal live tissue sections. Subregions from the dorsal hippocampus were extracted and microdissected using microsurgical tools. Microdissected pieces of tissue from each animal were pooled by subregion in 500 μL of prechilled TRIzol in a 1.5 mL microcentrifuge tubes and stored at −80°C.

Total RNA was extracted using the Direct-zol RNA miniprep kit (Zymo Research, Irvine, CA, United States) according to the manufacturer’s protocol. RNA quality was assessed on an Agilent 2200 TapeStation (Agilent Technologies, Palo Alto, CA, United States). Samples with a RIN greater than 7 were processed for library preparation with NEBnext rRNA depletion kit and ULTRAII FS RNA-seq Library Preparation Kit for Illumina (New England Biolabs, Ipswich, MA, United States). This procedure involved steps for mRNA enrichment, fragmentation, random primed cDNA synthesis, second strand synthesis, end repair, A-tailing, adaptor ligation, and PCR amplification. Size selection and cleanup was performed with SPRI select beads (Beckman Coulter, Indianapolis, IN, United States). Sequencing was performed on the NovaSeq 6000 (Illumina, Inc., San Diego, CA, United States) to obtain paired end sequencing reads.

2.4 Bulk RNA sequencing data preprocessing

RNA sequencing reads obtained from NovaSeq 6000 were systemically processed to ensure high-quality outputs for downstream analyses. To evaluate the quality of the reads, quality control analysis was performed using FASTQC v0.11.9. RNA sequencing reads in Binary Base Call (BCL) format were subsequently converted using bcl2fastq2 v2.20 to FASTQ format, and simultaneously demultiplexing to assign sequences to their respective samples based on index sequences. Demultiplexed samples that did not pass quality control were excluded from future analyses in the pipeline. Next, FASTQ files were then mapped to mouse MM10 (GRCm39) reference genome using STAR 2.7.10a aligner. Lastly, the transcript quantification of the reads was done using Salmon v1.2.0 to create a count matrix for downstream analyses.

2.5 Bulk RNA sequencing analyses

We conducted several differential expression analyses using DESeq2 v1.1.0 in the Illumina BaseSpace, comparing RNA counts between trained and untrained sample groups. Genes with a false discovery rate (FDR) below 0.05 were deemed differentially expressed. To identify the variances in our count data, we performed principal component analysis (PCA) on the normalized count data, encompassing all hippocampal regions and within each specific region.

To visualize our data, we generated a heatmap in ComplexHeatmap v2.18.0 using Z-scores, which standardized the expression levels of genes across samples (Gu et al., 2016). The heatmap, organized with dendrogram and sorted using Euclidean distance, allowed us to identify patterns of similarity in gene expression across samples.

To explore deeper into biological significance, we conducted Gene Ontology (GO) enrichment analysis using clusterProfiler v4.10.0 in R (Wu et al., 2021), focusing on biological processes, cellular components, and molecular functions. Only enriched GO terms with an FDR less than 0.01 were identified as significant, and the top 10 GO terms were displayed in a ranked dot plot. Redundant terms were filtered out using simplify with a cutoff of 0.7. Biological processes were ordered along the y-axis based on their ordered significance within each column from left to right. Additionally, if a biological process in one analysis is amongst the top 10 of another analysis, the dot signifying its enrichment will be shown at the corresponding row in the plot.

We also utilized Vennplex v1.0.0.2 software to compare and visualize the overlaps among differentially expressed genes (DEGs) from our differential expression analyses (Cai et al., 2013). The gene lists were categorized into upregulated and downregulated, allowing the visualizing of the gene overlaps and identification of counter-regulated genes (genes exhibiting opposite expression trends in different regions). We employed the SuperExactTest v1.1.0 in R to statistically assess whether the observed overlapping genes across different hippocampal regions significantly exceeded the expected overlap by random coincidence (Wang et al., 2015). The assessment was performed using hypergeometric distribution and Fisher’s Exact test, and the statistical significance of these gene overlaps was shown in upset-style plots.

2.6 Tissue preparation for spatial transcriptomics

60 min after retention test, mice were euthanized and brains were extracted and immediately prepared for snap freezing in cuvettes of Tissue-Tek Optimal Cutting Temperature compound (Sakura Finetek USA, Torrance, CA, United States) floating in a bath of methyl-butane chilled by liquid nitrogen. Brain blocks were cryosectioned at −20°C (10X Genomics, 2023) and 10 μm thick coronal sections containing the dorsal portion of the hippocampus were obtained (bregma −1.7 and − 2.2 mm, Paxinos and Franklin, 2019). Collected tissue sections are trimmed to fit within the 5 mm × 5 mm capture area on the Visium gene expression slide. Sections were mounted on a prechilled Visium gene expression slide or a Tissue Optimization slide.

For Optimization, all tissue sections were collected from a single tissue. For Visium gene expression, slide capture areas contained one section per biological sample. Slides were fixed in prechilled methanol before staining using the Visium spatial tissue optimization protocol (Genomics, 2019) or the Visium spatial gene expression protocol (Genomics, 2023).

Tissues permeabilization time on the gene expression slide was set to 18 min based on our tissue optimization results. Images were taken according to the Visium Spatial Gene Expression Imaging Guidelines (Genomics, 2021). Brightfield histology images for the Visium Tissue Optimization and Gene Expression slides were taken using the Leica Aperio CS2 Slide Scanner at 20x magnification. Tissue optimization fluorescent images were taken on a Zeiss LSM 800 (555 nm LED, 75% intensity, and 200 ms exposure).

mRNA was extracted and libraries were prepared following the Visium Spatial Gene Expression User Guide (Genomics, 2023). Libraries were indexed using the Dual Index Kit TT Set A (PN-1000215). Indexed libraries were loaded at 300 pM and sequenced on a NovaSeq 6000 System (Illumina) using a NovaSeq S4 Reagent Kit (200 cycles, catalog no. 20027466, Illumina), targeting a sequencing depth of approximately 2.0 × 108 read-pairs per sample. Sequencing depth was determined through a calculation of a suggestion 50,000 read pairs per spatial capture spot (Genomics, 2023). Our tissue samples covered roughly 80% of the 5,000 spatial capture spots within the Visium fiducial frame [(0.8 × 5,000) × 50,000 = 2.0 × 108]. Sequencing was performed using the following read protocol: read 1: 28 cycles; i7 index read: 10 cycles; i5 index read: 10 cycles; and read 2: 91 cycles.

2.7 Spatial transcriptomics preprocessing

Raw RNA sequencing data obtained from NovaSeq 6000 were systemically processed to ensure high-quality outputs for downstream analyses. Raw sequencing data in BCL format were converted using bcl2fastq2 v2.20 to FASTQ format, and simultaneously demultiplexing to assign sequences to their respective samples based on both i5 and i7 index sequences. Following the conversion, these FASTQ files were subsequently processed with TrimGalore v.0.6.5 to automate quality, adapter trimming and perform quality control.

For the alignment and the quantification of gene expression, Space Ranger v2.0.0 (10x Genomics) was used to generate and quantify raw gene-barcode matrices from RNA sequencing data. Alignment and quantification of Visium spatial gene expression data utilizes the same tools implemented in the analysis of single cell RNA-seq data sets. Data were aligned to refdata-gex-mm10-2020-A, a mouse genome index provided by 10x Genomics, which is an annotated version of mm10 genome assembly. The gene-barcode matrices identified the number of unique molecular identifiers (UMIs) associated with each gene for each cellular barcode, allowing the quantification of transcript abundance.

After the alignment and quantification, Space Ranger toolset was utilized to normalize the gene expression data within the raw gene-barcode matrices to account for any technical variations. Finally, we also utilized Space Ranger to align gene expression data with spatial coordinates derived from histological images of each sample, facilitating the visualization of transcriptional activity to specific tissue location.

2.8 Spatial transcriptomics analysis

Gene expression outputs from Space Ranger were analyzed using the software package Seurat v5 in the R coding environment (Hao et al., 2023; Butler et al., 2018). Spatial gene expression data were cropped to only include the capture spots within the dorsal hippocampus. Using default Seurat commands, data from each sample were normalized, transformed (using SCTransform v2) and integrated. The integrated data set were further normalized and clustered (resolution = 0.8) in PCA space. Clustered data can be visualized through UMAP projections of the capture spots or in the 2D space of the tissue section. Cell type annotation was performed by integrating our data with a random subset of 10,000 hippocampal cells from the Allen Brain Atlas Cortex and Hippocampus Single Cell Taxonomy (Yao et al., 2021). Predicted cell-type identities were calculated and optimized based on gene expression similarity and anatomical location.

Capture spots were grouped along two different criteria to calculate differential gene expression in our samples. First, capture spots belonging to the whole hippocampus and each anatomical cell layer in the DG, CA3 and CA1 regions were compared across training conditions. Next, within each training condition, hippocampal capture spots were separated into two groups based on the detectable (>0) expression of the immediate early gene (e.g., Arc + and Arc− groups). IEG expressing and not expressing groups of spots were compared both within training conditions and across training conditions. Differential gene expression was calculated using the Wilcoxon rank sum test (the default method for Seurat). Genes with an of FDR below 0.05 were considered differentially expressed genes/transcripts. GO term enrichment and overlap analyses on detected DEGs were performed as described above for bulk RNA sequencing DEGs.

Expression data for 23 IEGs from all hippocampal spots evaluated through pairwise Pearson’s Correlations using the R-language base stats package. Pearson’s correlation coefficients between all IEG pairs were calculated within samples, alongside p values to assess significance. Heatmaps of the correlation matrix were plotted using the R-language package, ComplexHeatmap.

2.9 RT-qPCR

Dorsal hippocampal subregions (CA1, CA3, and DG) were microdissected as described for bulk sequencing. Tissue was homogenized using TRIzol (Thermo Fisher Scientific, Waltham, MA, United States) before RNA extraction by guanidinium-phenol-chloroform extraction with ethanol precipitation was performed and quantified using spectrophotometry. cDNA was reverse transcribed from sample bulk RNA using Thermo Fisher’s SuperScript IV Kit (Thermo Fisher Scientific, Waltham, MA, United States). Relative gene expression was quantified using the SYBR green master mix system (Ponchel et al., 2003). 10 μL reactions were prepared in the wells of a 384 well plate, hippocampal regions were each prepared on separated plates. Each well contained 2.5 ng of template cDNA and 0.5 μM each of forward and reverse primers for one gene per well for a total reaction volume of 10 μL (sequences and melting temperatures can be found in Table 1). Each gene was loaded in triplicate in each plate to account for pipetting errors. Plates were incubated in the Bio-Rad CFX384 Real-Time PCR Detection System (Bio-Rad Life Sciences, Hercules, CA, United States) using the following schedule: 50°C for 5 min (1 cycle), 95°C for 10 min (1 cycle), 60°C for 1 min (1 cycle), and 95°C for 15 s followed by 60°C for 1 min (40 cycles).

Table 1. RT-qPCR primer sequences.

Primers were selected to validate the most significant genes with the highest fold change of expression for each regional comparison of trained and untrained animals. The genes Crym and Kit were selected from the DEGs detected in the DG comparison. Rasal1 and Pcsk1n were selected from the DEGs detected in the CA3 comparison. And Bok, Homer3, Serpina3n, and Shank3 were selected from the DEGs detected in the CA1 comparison.

Ct values for each reaction are automatically determined (roughly 20% of the plateau value). Data were analyzed for relative gene expression quantification utilizing the ∆∆Ct method, using the expression of Gapdh as the internal control. Relative expression values were normalized to the expression value of the untrained group. The Levene’s test was uses to assess the homogeneity of variances for Log2 transformed relative expression. Upon confirmation of homogeneity of variances, Student’s t-tests were performed on each primer target across training conditions using log2 transformed relative expression value for each comparison.

3 Results 3.1 APA trained mice exhibit place avoidance memory

To investigate gene expression changes induced by active place avoidance (APA) training in the hippocampal network, we conducted bulk RNA-sequencing and spatial transcriptomics (see Figure 1A for experimental timeline). In the APA, mice learn to avoid a 60° shock zone on a rotating circular arena (Cimadevilla et al., 2000; Stuchlik et al., 2013). Mice trained in the APA task with the shock zone on (n = 4, 2 females, 2 male) exhibited a decrease in the number of detections in the shock zone during the acquisition trials (Figure 1B) and longer re-entry times to the shock zone during the memory retention test trials (Figure 1C) relative to the mice in the untrained cohort (n = 4, 1 female, 3 males). Re-entry times to the shock zone (i.e., time to second entrance) was used as the retention test performance metric because it represents performance that is more resistant to the animal’s potential errors in self-localization once placed in the rotating arena (Cimadevilla et al., 2000; Stuchlik et al., 2013; Cimadevilla et al., 2001). From this pool of animals 3 trained (2 females, 1 male) and 3 untrained (1 female, 2 male) were utilized in Bulk RNA-sequencing experiments and 1 trained and 1 untrained animal (2 males) were utilized in spatial transcriptomic experiments.

Figure 1. Mice trained in the APA learned to avoid the location of an unmarked shock zone. (A) Experimental timeline. Mice were trained in the active place avoidance paradigm over the course of 4 days following 3 days of handling. On day 0 mice were habituated (HT) to a rotating arena with no active shock zone. On day 1 and 2 mice assigned to the trained (TR) behavioral condition (n = 4) received two 30-min trials with an active shock zone. Mice assigned to the untrained (UT) behavioral condition (n = 4) were not exposed to an active shock zone. On day 3 all mice received a 10-min retention test trial (exposed to the rotating arena with no active shock zone). Mice were sacrificed 60 min following the end of the retention test, brains were collected for further processing. Samples from three of the four trained mice (2 females, 1 male) and untrained mice (1 female, 2 males) were processed with bulk RNA sequencing, while the remaining mice (1 trained male, 1 untrained male) had their samples processed with spatial transcriptomics. (B) Training performance measured as number detections in the shock zone. Mice assigned to the trained behavioral condition learned to avoid the location of the shock zone indicated by the decreased number of shocked triggered. Exponential fit learning curves included with shaded error bars (SE). (C) Memory performance measured as time to second entrance in the retention test trial. Trained mice demonstrated higher time to second entrance during the retention test. Trained mean 375 s ± 70.4 s (SE), Untrained mean 58.7 s ± 8.51 s (SE). Line and point color reflect the assigned behavioral condition (untrained = blue, trained = red). Point shape reflects the experimental pipeline a given animal was assigned to (square = bulk RNA-seq, circle = spatial transcriptomics).

3.2 Hippocampal subregional transcriptomes are stratified by training condition

To explore the molecular mechanisms following memory retention test across the dorsal hippocampus, bulk RNA-sequencing was performed on the Dentate Gyrus (DG), CA3, and CA1 hippocampal subregions collected from trained and untrained mice (trained: 2 females, 1 male; untrained: 1 female, 2 male). Principal component analysis (PCA) from all 3 subregions combined reveals a separation of samples by behavioral cohort (Figure 2A), indicating that differences in training contribute substantially to the overall variability of gene expression in our samples. Separation between samples is further accentuated when PCA is performed separately on data from each subregion (Figures 2B–D). Consistently, hierarchical clustering of gene expression profiles reliably grouped samples by their behavioral cohort (Figure 2E).

Figure 2. APA training results in substantial changes in hippocampal subregion gene expression profiles. (A–D) Left, principal component analyses of all hippocampal subregions (A), CA1 (B), CA3 (C), DG (D) illustrate a clear separation between samples collected from trained (red) and untrained (blue) mice along the axis that explains the most variance. Hippocampal regions are also indicated by the point shape, squares refer to the CA1, circles to the CA3, and triangles to the DG. Right, heatmap with top 20 contributing genes to the first two principal components of each are reported in the heat map to the right of each plot. Sample size for the CA1 was 3 trained and 3 untrained. Sample size for the CA3 was 3 trained and 2 untrained. And Sample size for the DG was 3 trained and 2 untrained. Two samples from separate untrained animals (1 DG, 1 CA3) did not pass sequencing quality control and were excluded from PCA analyses. (E) Normalized gene expression of the top 2000 most significant DEGs. Hierarchical clustering of gene expression profiles plotted as a dendrogram clusters samples by training condition.

3.3 Hippocampal subregional enrichment of learning and memory-related biological processes

We carried out differential gene expression analysis comparing APA-trained to untrained mice. We detected 1,663, 606, 315, and 350 differentially expressed genes (DEGs) (adjusted p-value <0.05) in the combined hippocampal subfields, CA1, CA3 and DG respectively, when contrasting APA-trained and untrained mice (Figures 3A–D and Supplementary Figures 1A–D). Despite high statistical significance (FDR < 0.005), most DEG in each trained versus untrained comparison exhibited modest fold change (|log2 fold| < 1.0) difference. Next, we investigated the overlap of regionally detected DEGs (FDR < 0.05) to elucidate the distribution of DEGs across hippocampal subregions. We found the largest overlap of upregulated DEGs between the CA3 and CA1 subregions with a total of 127 intersecting genes (Figure 3E and Supplementary Figure 1E).

Figure 3. Bulk RNA-seq of hippocampal subregions reveals similar expression profiles in CA1 and CA3 regions in trained samples. (A–D) Volcano plot of DEGs between trained and untrained animals in all hippocampal subregions combined (A), and in the CA1 (B), CA3 (C), and DG (D) subregions individually. Sample size for the CA1 was 3 trained and 3 untrained. Sample size for the CA3 was 3 trained and 2 untrained. And Sample size for the DG was 3 trained and 2 untrained. Two samples from separate untrained animals (one DG, one CA3) did not pass sequencing quality control and were excluded from differential gene expression analyses. In all hippocampal subregions combined, 966 genes were upregulated and 697 were downregulated. 393 genes were upregulated and 213 were downregulated in the CA1, 273 genes were upregulated and 42 were downregulated in the CA3, and 117 genes were upregulated and 223 were downregulated in the DG. See Supplementary Figures 1A–D for full sized plots. (E) Overlaps in regional DEGs demonstrate greatest similarity between CA1 and CA3 subregions. DEGs from each regional analysis were stratified by direction of fold change. Overlaps in DEGs between the CA3 and CA1 subregions were significant (Fisher’s Exact Test). See Supplementary Figure 1E for further details. (F) Regional enrichment of biological processes detected amongst all DEGs (left) and stratified by up- (middle) and down-regulated (right) DEGs. The top 10 biological processes detected in each analysis are shown on the y-axis. Biological processes are ordered based on their statistical significance in the left-most column in which they are detected. Dot color reflects the statistical significance [−log10(FDR)] of the biological process enrichment. Dot size reflects the number of detected DEGs mapped to the genes involved in a given biological process.

GO enrichment analysis was conducted to identify overrepresented biological process among DEGs from the trained versus untrained comparison and identified within each hippocampal subregion (The Gene Ontology Consortium, 2019). Among the DEGs detected in the combined hippocampal subregions, we found overrepresentation of genes involved in the modulation of neurotransmitter signaling (Figure 3F). These DEGs included Ywhag, Bsn, and Nup153, which have been previously implicated in these cellular functions (Wachi et al., 2017; Annamneedi et al., 2021; Leone et al., 2019). DEGs detected in both the CA3 and CA1 were also enriched with genes involved in the modulation of neurotransmitter signaling. The genes Snap25, Shank1, Shank 3, Bsn and Ywhag were identified in the overlap of upregulated DEGs between these regions and have been studied for their involvement in synaptic transmission (Irfan et al., 2019; Mao et al., 2015; Shi et al., 2017; Kouser et al., 2013). DEGs detected in the DG were enriched with genes involved in neurogenesis and transcriptional regulation. Among these DEGs were Kit, H2afz, and Kmt2e, which have been studied for their relationship to neurogenesis and transcriptional regulation (Katafuchi et al., 2000; Ferreira et al., 2016; Barbieri et al., 2018; Zovkic et al., 2014; Collins et al., 2019). When enrichment analyses were separated by up- and down-regulated DEGs, the differences separating the CA1 and CA3 from the DG became more striking with no significant enrichment of biological processes among the downregulated DEGs of the CA1 and CA3. These data reveal a spatial distribution of gene expression across the dorsal hippocampus of trained animals, defined by a greater degree of similarity between the CA1 and CA3 than between either of these subregions and the DG.

3.4 RT-qPCR validation of regionally detected DEGs

We performed targeted gene expression analyses using RT-qPCR on DG, CA3, and CA1 subregion samples collected from similarly behaviorally conditioned (i.e., trained and untrained) animals (Supplementary Figures 3A,B). We probed a total of 8 genes (with Gapdh as an internal control), each identified for their high statistical significance and relatively large positive fold change of expression in one of the region-specific comparisons between trained and untrained samples. In the DG, we identified significant differential expression of the gene Kit between trained and untrained animals, which replicates our results from bulk RNA sequencing comparisons of in this region (Supplementary Figure 3C). In the CA3 we identified significant differential expression of the gene Rasal1 between trained and untrained animals, replicating our results from bulk RNA sequencing (Supplementary Figure 3D). In the CA1, Homer3 and Crym were significantly differentially expressed (Supplementary Figure 3E). Differential expression of Homer3 in the CA1 was an expected validation of our bulk RNA sequencing findings, but the detection of Crym differential expression was unexpected as it was selected as a candidate DEG from the DG. In sum, RT-qPCR allowed us to validate the differential expression of select target genes with the largest expression fold changes in our bulk RNA-seq dataset. Consistent with prior reports, some of our target genes which exhibited comparatively lower fold changes of expression escaped detection because of the limitation in RT-qPCR relative to RNA-sequencing approaches (Coenye, 2021; Everaert et al., 2017; Reshetnikov et al., 2020).

3.5 Differences in the spatial distribution of hippocampal gene expression between trained and untrained animals

To further investigate the spatial distribution of DEGs associated with APA memory training across the hippocampus, we performed spatial transcriptomics on coronal sections containing the dorsal hippocampus from one trained (male) and one untrained (male) mouse. Computational analyses of integrated capture spots in the hippocampal slices from trained and untrained animals revealed distinct clusters (seen in a Uniform Manifold Approximation and Projection (UMAP) plot) which map along anatomical boundaries in a spatial transcriptomic plot (Figures 4A,B). Spatial transcriptomics data were integrated with the Allen Brain Atlas (Yao et al., 2021) to computationally annotate each spot with the cell-type most prominently detected (Figures 4C,D).

Figure 4. DEGs detected between trained and untrained animals within hippocampal subregions using spatial transcriptomics are enriched with genes involved in similar biological processes to those seen in bulk RNA-seq. (A,B) UMAP plot (A) and spatial transcriptomic maps of hippocampal spots (B) from trained and untrained mouse samples (n = 1 per group). Discrete color scale displays unsupervised graph-based clustering of hippocampal spots. Clusters are labeled 1 through 7 from greatest (1) to least (7) number of spots belonging to each unbiased cluster. Unbiased clustering groups spots along anatomical boundaries reflecting natural transcriptomic differences between morphologically distinct regions in the hippocampus. (C,D) UMAP plot (C) and cell-type annotated spatial transcriptomic maps of hippocampal spots (D) in samples from the trained and untrained mouse. Discrete color scale reflects the predominant cell type surveyed by each capture spot which was determined by cell-type annotation of hippocampal capture spots through automated integration of our data with the single cell Allen Brain Atlas for mice. (E–H) Volcano plot of DEGs between all hippocampal (E), CA1 (F), CA3 (G), and DG (H) cell-layer spots from trained and untrained mice (n = 1 per group). In the analysis of all hippocampal combined, 205 genes were upregulated and 147 were downregulated. 55 genes were upregulated and 20 were downregulated in the CA1 cell layer, 54 genes were upregulated and 18 were downregulated in the CA3 cell layer, and 129 genes were upregulated and 58 were downregulated in the DG cell layer. (I) Overlap of regionally detected DEGs stratified by direction of fold change. Overlaps were tested with Fisher’s exact test. See Supplementary Figure 2 for further details. (J) Regional enrichment of biological processes detected amongst all DEGs (left) and stratified by up- (middle) and down-regulated (right) DEGs. The top 10 biological processes detected in each analysis are shown on the y-axis. Biological processes are ordered based on their statistical significance in the left-most column in which they are detected. Dot color reflects the statistical significance [−log10(FDR)] of the biological process enrichment. Dot size reflects the number of detected DEGs mapped to the genes involved in a given biological process. See Supplementary Figure 2 for enrichment of Cellular Component and Molecular Function GO terms.

Differential gene expression analysis between trained and untrained mice identified a total of 352 DEGs in all hippocampal spots, with further stratification by subregional cell layers detecting 75 DEGs in the CA1, 72 DEGs in the CA3, and 187 DEGs in the DG (Figures 4E–I and Supplementary Figure 2). GO enrichment analyses of the DEGs detected between all hippocampal spots identified biological processes related to energy production and protein expression (Figure 4J). DEGs detected in the analysis of the CA1 cell layer were enriched with biological processes related to energy production and synaptic function. In the CA3, DEGs were enriched with genes involved in neurotransmitter release and energy production. And in the DG, DEGs were enriched with genes involved in energy production and translational machinery. Notably, the biological process “regulation of synaptic plasticity” (GO: 0048167) was enriched among the DEGs detected between all subregional comparisons between trained and untrained animals.

These spatial transcriptomics data are concordant with our bulk RNA sequencing data. The overlap of DEGs detected in each hippocampal region between the two techniques was lower than expected (10, 10, 22, in the CA1, CA3, and DG, respectively). Notwithstanding, there is stronger concordance between the biological process enrichment detected with the two techniques. We found upregulation of genes involved in the regulation of synaptic plasticity across all three hippocampal subregions in both experiments. These findings reinforce a subregion-specific model of memory-associated transcriptional activation following memory recall.

3.6 Arc-expressing hippocampal spots exhibit changes in the expression of memory-associated genes and biological processes

Our data suggest that differences in behavioral conditioning contributed to variations in the spatial distribution of synaptic plasticity related gene expression across the hippocampal network. It is well established that a surge of Arc expression accompanies memory associated processes such as consolidation and retrieval, and differences in behavioral conditioning affect the recruitment of Arc-expressing hippocampal neurons which form the memory-associated neuronal ensemble (Gouty-Colomer et al., 2016; Gallo et al., 2018; Carrillo-Reid, 2022; Gouty-Colomer et al., 2016; Hsieh et al., 2021; Chia and Otto, 2013). As such, we studied the Arc-expressing spatial transcriptomics spots to assess changes in biological processes within the memory associated neuronal ensemble following APA memory recall. Spatial transcriptomic spots with detectable expression of Arc mRNA (>0 copies) were determined to be Arc + in this experiment. Given the resolution and capture efficiency of spatial transcriptomics (Asp et al., 2020), these spots represent regions of space with high expression of Arc mRNA across one or more cells/cellular compartments. Hence, we utilized these Arc + spots as a proxy for the location of the putative memory-associated neuronal populations with induced expression of Arc following the retrieval of the active place avoidance memory.

We observed a sparse distribution of Arc + spots in our samples consistent with the characteristic expression of Arc protein in the hippocampus from studies investigating memory associated ensembles (Mayford, 2014; Chia and Otto, 2013; Gouty-Colomer et al., 2016). Out of the 542 spots in the hippocampus from the trained animal, 195 spots had detectable expression of Arc mRNA (Arc+). The untrained animal, 169 Arc + spots out of 568 hippocampal spots (Supplementary Figure 4C). In both samples, spots with the highest level of Arc expression concentrated primarily in the CA1 cell layer (Figures 5A,B).

Figure 5. Arc-positive spots in the trained sample are enriched with genes involved in synaptic plasticity. (A,B) UMAP plot (A) and spatial transcriptomic maps of hippocampal spots (B) in the trained and untrained samples colored to reflect normalized expression of Arc. Arc positivity is determined by detectable expression (>0 copies) of Arc mRNA in a spatial transcriptomic spot. (C,D) Volcano plots of differential gene expression of Arc positive (C) and negative (D) spots across behavioral conditions. In the analysis of Arc + spots across behavioral conditions 85 upregulated and 89 downregulated genes were detected. In the analysis of Arc- spots across behavioral conditions 113 upregulated and 240 downregulated genes were detected. (E,F) Volcano plots of differential gene expression between Arc positive and negative spots within trained (E) and untrained (F) samples. In the analysis of Arc spots within the trained animal 95 upregulated and 6 downregulated genes were detected. In the analysis of Arc spots within the untrained animal 47 upregulated and 31 downregulated genes were detected. (G) Biological processes overrepresented in the analysis of Arc-expressing spots within trained and untrained samples amongst all DEGs (left) and stratified by up- (middle) and down-regulated (right) DEGs. The top 10 biological processes detected in each analysis are shown on the y-axis. Biological processes are ordered based on their statistical significance in the left-most column in which they are detected. Dot color reflects the statistical significance [−log10(FDR)] of the biological process enrichment. Dot size reflects the number of detected DEGs mapped to the genes involved in a given biological process. (H) Biological processes overrepresented in the analysis of Arc-expressing spots across trained and untrained samples amongst all DEGs (left) and stratified by up- (middle) and down-regulated (right) DEGs. The top 10 biological processes detected in each analysis are shown on the y-axis. Biological processes are ordered based on their statistical significance in the left-most column in which they are detected.

To assess synaptic plasticity-related transcriptional changes in these Arc + spots, we conducted multiple differential expression comparisons: (1) across sets of Arc + or Arc− spots between the trained and untrained animal (Figures 5C,D), and (2) within the trained or untrained animal between sets of Arc + and Arc− (Figures 5E,F). The comparison of Arc + spots across trained and untrained animals detected 174 significant DEGs, and comparison of Arc− spots across trained and untrained animals detected 353 significant DEGs. GO term enrichment analysis of DEGs detected across Arc + spots in the trained and untrained animal revealed enrichment with biological processes related to ribosomal function and energy production (Figure 5H). A finding which is consistent with the regionalized differences observed between the trained and untrained animals through spatial transcriptomics (Figure 4J). Rpl4 and Rpl5 were differentially expressed in the analysis across Arc + spots in the trained and untrained, and are known for their roles in ribosomal biogenesis (Robledo et al., 2008). Enrichment analyses of DEGs detected in Arc− spots across training conditions revealed a similar set biological processes related to energy production as those detected in the analysis of Arc + spots. Differences between the two comparisons are accentuated when enrichment analyses are performed separately on up- and down-regulated DEGs.

101 DEGs were detected within Arc + and Arc− spots in the trained animal, and 78 DEGs were detected within Arc + and Arc− spots in the untrained animal. Notably, substantially more biological processes related to synaptic plasticity and neuronal excitability were enriched in the DEGs Arc + spots of the trained animal as compared to Arc + spots of the untrained animal. This difference is enhanced when the enrichment analyses are separated by up- and down-regulated DEGs, where only the upregulated DEGs within the Arc + spots of the trained animal demonstrate enrichment of any biological processes (Figure 5G). Prkcb and Dkk3, known for their role in synaptic plasticity (Weeber et al., 2000; Martin-Flores et al., 2024), were detected amongst upregulated genes in the comparison of Arc + and Arc− spots in the trained animal. In summary, enrichment analyses in Arc + spots across training conditions detected biological processes related to energy metabolism and ribosomal biogenesis, while Arc + spots within the trained sample revealed biological processes related to synaptic plasticity and transmission.

3.7 IEG-expressing hippocampal spots comprise distinct gene expression profiles

Our findings demonstrate that behavioral training induces discrete transcriptional changes in the population of spots marked by Arc expression. However, memory traces are likely comprise overlapping yet molecularly distinct neuronal populations tagged by IEGs (Sun et al., 2020; Nambu et al., 2022; Okuno, 2011; Gallo et al., 2018). We sought to examine spots marked by the expression of the IEG most and least correlated with the expression of Arc to evaluate the functional relationship between these subsets of the memory-associated ensemble. Pairwise correlations identified Egr1 and c-Jun as the IEGs most and least correlated with Arc, respectively (Supplementary Figures 4A,B, 5A–F).

Unlike Arc, which codes for a cytoskeleton regulating protein, both Egr1 and c-Jun code for transcription factors. IEGs which function as transcription factors are thought to promote the expression of activity related genes which execute the cellular level changes necessary for memory storage (Yap and Greenberg, 2018; Sun et al., 2020; Kim et al., 2018; Yelhekar et al., 2024). Genes downstream of Egr1-driven transcription are involved in vesicular release, neurotransmitter metabolism, receptor expression, and synaptic plasticity (Duclot and Kabbaj, 2017; Sekiya et al., 2019). Genes downstream of c-Jun-driven transcription in the brain are involved in neurite growth, axonal regeneration, and synaptic long term depression (Raivich and Behrens, 2006;