Abstract
A fundamental, evolutionarily conserved biological mechanism required for long-term memory formation is rapid induction of gene transcription upon learning in relevant brain areas. For episodic types of memories, two regions undergoing this transcription are the dorsal hippocampus (dHC) and prelimbic (PL) cortex. Whether and to what extent these regions regulate similar or distinct transcriptomic profiles upon learning remain to be understood. Here, we used RNA sequencing in the dHC and PL cortex of male rats to profile their transcriptomes in untrained conditions (baseline) and at 1 h and 6 d after inhibitory avoidance learning. We found that, of 33,713 transcripts, >14,000 were significantly expressed at baseline in both regions and ∼3000 were selectively enriched in each region. Gene Ontology biological pathway analyses indicated that commonly expressed pathways included synapse organization, regulation of membrane potential, and vesicle localization. The enriched pathways in the dHC were gliogenesis, axon development, and lipid modification, while in the PL cortex included vesicle localization and synaptic vesicle cycle. At 1 h after learning, 135 transcripts changed significantly in the dHC and 478 in the PL cortex; of these, only 34 were shared. Biological pathways most significantly regulated by learning in the dHC were protein dephosphorylation, glycogen and glucan metabolism, while in the PL cortex were axon development and axonogenesis. The transcriptome profiles returned to baseline by 6 d after training. Thus, a significant portion of dHC and PL cortex transcriptomic profiles is divergent, and their regulation upon learning is largely distinct and transient.
SIGNIFICANCE STATEMENT Long-term episodic memory formation requires gene transcription in several brain regions, including the hippocampus and PFC. The comprehensive profiles of the dynamic mRNA changes that occur in these regions following learning are not well understood. Here, we performed RNA sequencing in the dorsal hippocampus and prelimbic cortex, a PFC subregion, at baseline, 1 h, and 6 d after episodic learning in rats. We found that, at baseline, dorsal hippocampus and prelimbic cortex differentially express a significant portion of mRNAs. Moreover, learning produces a transient regulation of region-specific profiles of mRNA, indicating that unique biological programs in different brain regions underlie memory formation.
Introduction
A fundamental and evolutionarily conserved mechanism required for the formation of long-term memory is an initial phase of de novo transcription and translation in the relevant nervous system regions (Klann and Dever, 2004; Costa-Mattioli et al., 2009; Alberini and Kandel, 2014). These processes are necessary for stabilizing a labile memory representation, a process known as memory consolidation (McGaugh, 2000; Dudai, 2012; Squire et al., 2015).
For episodic types of memories, which store spatial and contextual information, de novo transcription required for memory consolidation occurs in the medial-temporal lobe regions, which include the hippocampus (HC) as well as functionally linked cortical regions, such as the prefrontal cortex (PFC) (Dash et al., 2004; Wiltgen et al., 2004; Alberini, 2009). Over time, the HC appears to become dispensable for long-term memory storage, while cortical regions remain functionally engaged, suggesting that distinctive regional biological mechanisms must be responsible for memory consolidation. This temporal redistribution of the memory representation, which is critical for maintaining and storing the learned information long-term, is known as systems-level consolidation (Frankland and Bontempi, 2005; Goshen et al., 2011; Dudai et al., 2015; Squire et al., 2015).
To fully understand how long-term memories are formed and stored, we need a comprehensive understanding of the transcriptomic profiles that change in response to learning across the brain regions that constitute the activated memory system, as well as of their progression over time. In other words, we need to obtain a system(s) biology understanding of memory, which is still in large part lacking. Furthermore, for decades, the molecular studies on long-term memory consolidation focused mostly on the characterization of mechanisms in single brain regions, particularly the HC. Only more recently has some attention been turned to cortical regions. Comprehensive, unbiased biological assessments, such as RNA sequencing (RNAseq), have been conducted but have been mostly limited to single brain regions. For example, Cho et al. (Cho et al., 2016; Mathew et al., 2016) used RNAseq together with ribosomal profiling (active translational profiles) in the mouse HC after contextual fear conditioning and found an initial translational wave at 5 min after training followed by multiple transcriptional waves, beginning 10–30 min after training and continuing for 4 h. Bero et al. (2014) performed RNAseq in the ACC, a subregion of the PFC involved in memory storage, at 1 h after contextual fear conditioning in mice and reported a significant upregulation of synaptic plasticity transcriptional profiles paired with dendritic spine remodeling. Rizzo et al. (2017), using RNAseq on polyribosomal-associated mRNAs, reported two distinct waves of mRNA recruitment to actively translating ribosomes in the PFC (ACC, prelimbic [PL], and infralimbic cortices) at 1 and 6 h after contextual fear learning in mice. These and other previous studies emphasized the general concept that brain plasticity mechanisms underlie memory formation. In most cases, these studies focused their attention on the fundamental pathways shared across brain regions and systems; however, whether distinct mechanisms are also engaged in each brain region remains to be determined. Halder et al. (2016) used RNAseq, ChIP-seq, and methylated DNA immunoprecipitation-seq in the hippocampal CA1 subregion and ACC after contextual fear conditioning in mice at 1 h and 4 weeks after training. They reported changes in transcription and methylation profiles of plasticity genes in both CA1 and ACC at 1 h after training and differential methylation changes only in the ACC at 4 weeks after training (Halder et al., 2016). Other studies compared cell type-specific transcriptional profiles in different regions, such as the HC and PFC in mice, but only at baseline (Zhang et al., 2014; Zeisel et al., 2015; Cembrowski et al., 2016). These results provided important cell type-specific transcriptional databases, but with the caveat that cell-sorting techniques damage cytoplasmic and peripheral cell compartments (van den Brink et al., 2017; Nguyen et al., 2018) where rapid mRNA changes occur in response to learning (Miyashiro et al., 1994; Zhong et al., 2006; Ostroff et al., 2019).
In sum, most RNAseq studies have been limited to single brain regions and have investigated different learning paradigms and time points after training, hindering cross-study comparative analyses. To better understand the systems-level biology of memory formation, it is important to compare the transcriptional changes that take place across brain regions to identify common and/or distinct biological processes, and to use the same learning paradigm and species. Furthermore, as transcription required for memory consolidation is known to be transient (Alberini and Kandel, 2014; Santini et al., 2014), the temporal progression of the transcriptional changes induced by learning also needs to be elucidated. Finally, as most RNAseq studies have been conducted using mice, little is known about transcriptional profiles across brain regions in other species.
To begin addressing these gaps in knowledge, we used RNAseq to compare dorsal HC (dHC) and PL cortex transcriptomics at baseline, and at 1 h and 6 d following inhibitory avoidance (IA) learning in rats, an episodic aversive event that evokes long-term memory after a single trial.
Materials and Methods
Animals
Adult male Long-Evans rats (200-250 g) were used for all experiments. Animals were housed in pairs and maintained on a 12 h light–dark cycle. Experiments were performed during the light cycle. All rats were allowed ad libitum access to food and water and were handled for 3 min per day for 5 d before behavioral manipulations. All protocols complied with the National Institutes of Health's Guide for the care and use of laboratory animals and were approved by the New York University Animal Welfare Committee.
IA behavior
IA was conducted as previously described (Taubenfeld et al., 2001). The IA chamber (Med Associates) consisted of a rectangular Perspex box divided into a safe compartment and a shock compartment. The safe compartment was white and illuminated by a light fixture fastened to a wall of the compartment. The shock compartment was dark and made of black Perspex. Footshocks were delivered to the grid floor of the shock compartment via a constant current scrambler circuit. The two compartments were separated by an automatically operated sliding door. The apparatus was located in a sound-attenuated, nonilluminated room. During training sessions, each rat was placed in the safe compartment facing away from the door. After 10 s, the door separating the compartments was automatically opened, allowing the rat access to the shock compartment. The door closed 1 s after the rat entered the shock compartment, and a 2 s 0.6 mA footshock was delivered. Ten seconds after delivery of the footshock, the rat was returned to its home cage. Untrained rats were handled but otherwise remained in the home cage. Animals were killed 1 h or 6 d after IA training, or at matched time points for untrained animals. Unpaired-shock control rats were exposed to the IA chamber and allowed to enter the shock compartment; however, a footshock was not administered upon the sliding door being closed. One hour later, the animals were placed in a different box, subjected to a footshock (2 s, 0.6 mA), immediately removed and returned to the homecage. This novel context had different dimensions and was located in a different behavioral room entirely to avoid familiarity of external cues. Animals were killed 1 h after receiving the unpaired footshock.
Tissue dissection, isolation of RNA, RNA sequencing alignment, and read counting
Rats were killed by decapitation. Their brains were quickly dissected and snap-frozen in prechilled isopentane on dry ice. To obtain individual PL cortex and dHC samples for each brain, we combined the bilateral PL cortices or dHC regions using a Neuro Punch (19 gauge, Fine Science Tools), as described previously (Milekic et al., 2007). The RNA of each sample was isolated using the RNAeasy Plus Universal Kit following the manufacturer's protocol (QIAGEN). All experiments were validated by confirming a training-dependent increase of the immediate early gene c-fos, a well-established marker of neural activation (Tischmeyer and Grimm, 1999) using RT-PCR (primer information in Extended Data Fig. 4-1). Experimental conditions, n per group, and workflow are described in Figure 1. The samples were submitted to the New York University School of Medicine Genome Technology Center for further quality control, quantification of mRNA read counts, and normalization. RNA integrity numbers (RINs) for each pool of RNA (see Fig. 1) were determined using an Agilent Technologies 2100 Bioanalyzer. RINs are generated by an algorithm that measures values on a scale of 1 to 10, based on the integrity of RNA with 10 being the highest. For all pools, the minimum RIN was 8.6, the maximum RIN was 9.9, with an average across all pools of 9.5. RNAseq libraries were obtained using the Illumina TruSeq Stranded mRNA LT kit (catalog #RS- RS-122-2101 or RS-122-2102) on a Beckman Biomek FX instrument, using 200 ng of total RNA as input, amplified by 12 cycles of PCR, and sequenced by an Illumina 4000 as paired end 150 reads. The reads were mapped to the rat reference genome (Rnor_6.0.84) using the STAR aligner (version 2.5.0c) (Dobin et al., 2013). Alignments were guided by a Gene Transfer File (Ensembl GTF Rnor_6.0.84), and the mean read insert sizes and their SDs were calculated using Picard tools (version 1.126) (http://broadinstitute.github.io/picard/). The read per million normalized BigWig files were generated using BEDTools (version 2.17.0) (Quinlan and Hall, 2010) and bedGraphToBigWig tool (version 4). The read count tables were generated using HTSeq (version 0.6.0) (Anders et al., 2015), normalized based on their library size factors using DEseq2, and differential expression analysis was performed (Love et al., 2014). Expression of individual transcripts were deemed significantly different from untrained at false discovery rate (FDR) < 0.1.
RNAseq data validation using RT-PCR
RNAseq data were validated by RT-PCR on mRNA extracted from distinct cohorts of rats. Ten transcripts were chosen from the list of genes found enriched in each region at baseline. Twenty to 22 transcripts were chosen from the list of differentially expressed genes (DEGs) at 1 h after training, and 10 transcripts were chosen from the data obtained at 6 d after training. All primer sequences were selected either from previous literature on brain tissue in rats or designed using Primer3 (http://bioinfo.ut.ee/primer3/). Primer specificity against the target sequence was confirmed using NCBI Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) and tested for primer dimers using a DNA gel to demonstrate that their PCR product generated a single band of the expected size. Finally, the denaturation curve performed after RT-PCR cycling with an intercalating dye was analyzed to ensure a single distinct peak in the plot of the negative derivative of fluorescence versus temperature. This indicated that the amplified double-stranded DNA products were not because of nonspecific amplification. RT-PCR analyses were performed using CFX96 Touch Real-Time PCR Detection System (Bio-Rad) with iQ SYBR Green Supermix (Bio-Rad). For each sample, 500 ng of cDNA was amplified using one initial denaturation step at 95°C for 5 min, followed by 40 cycles at 95°C for 15 s, 59°C for 30 s, and 72°C for 20 s. Triplicates of each cDNA amplification were analyzed by RT-PCR, and the averaged cycle threshold (Ct) value was used for relative quantifications (see Critical Factors for Successful Real-Time RT-PCR by QIAGEN). Samples that contained Ct values >36 were excluded from analysis, considering the loss in RT-PCR calculation accuracy at this level (VanGuilder et al., 2008; Caraguel et al., 2011). Transcripts enriched in the PL cortex and dHC at baseline were normalized to the level of ube2d2 of each respective sample (because the levels of gapdh were significantly different between regions), whereas genes differentially expressed within each region following IA training and in unpaired conditions were normalized to the level of gapdh of each respective sample (all primer sequences listed in Extended Data Figs. 2-5, 4-1).
Gene expression analysis, Gene Ontology (GO) biological process analysis, and gene set enrichment analysis (GSEA)
Heatmaps, volcano plots, and gene-concept maps were generated using clusterProfiler (Yu et al., 2012) packages in R/Bioconductor (Gentleman et al., 2004; Huber et al., 2015). The Database for Annotation, Visualization and Integrated Discovery (DAVID, version 6.8) (Huang et al., 2009) was used for the functional description of gene expression changes using GO biological process terms. Only GO biological processes with an FDR of <5% were considered. Dotplots were generated using the “compareCluster” function in the clusterProfiler package, and gene-concept networks were generated using the “cnetplot” function in the clusterProfiler package. GSEA was performed using GSEA (Broad Institute, version 4.1.0) (Mootha et al., 2003; Subramanian et al., 2005). A preranked analysis was performed using log2 fold change as the ranking metric. Gene sets for GSEA were derived from the Molecular Signatures Database (MSigDB, version 7.0) using the C5 collection of GO biological processes. Only gene sets with an FDR of <25% were considered, following guidelines set by the Broad Institute (https://www.gsea-msigdb.org/gsea/doc/GSEAUserGuideFrame.html).
Cell type- and synaptic compartment-specific analyses
To identify cell type-specific gene expression in baseline and learning-induced transcriptomic profiles, we compared our lists of transcripts to a database of cell type-specific mRNA expression published in Zeisel et al. (2015), which established selectively enriched transcripts in astrocytes, ependymal cells, endothelial cells, interneurons, microglia, mural cells, and oligodendrocytes. To determine the transcripts expressed in dendritic compartments, we used gene lists compiled by Middleton et al. (2019). Finally, to determine the transcripts related to synaptic specialization, we used the SynGO portal (syngoportal.org) published by Koopmans et al. (2019).
Statistical analyses
RNAseq data statistical analyses and GO biological process analyses are described in the above sections, Tissue dissection, isolation of RNA, RNA sequencing alignment, and read counting, and Gene expression analysis, Gene Ontology (GO) biological process analysis, and gene set enrichment analysis (GSEA). Pairwise comparisons for RT-PCR validations were performed using Student's t tests in GraphPad Prism for MacOS (version 9.0.0). p < 0.05 was considered statistically significant for pairwise comparisons. p values, t values, and degrees of freedom are shown in Extended Data Figs. 2-6, 4-2, 4-3, 4-4, 4-5, 9-3, and 9-4.
Availability of data
All datasets are available in the GEO repository under accession number GSE150950.
Results
RNAseq was conducted in dHC and PL cortex of untrained rats (also referred to as baseline conditions), and rats killed at 1 h and 6 d following IA training. In IA, the animals learn to avoid a context paired with a footshock, hence forming a long-term episodic aversive memory after a single learning trial. This discrete learning event allows for well-defined temporal investigations. The RNA preparations from all groups of rats that underwent RNAseq were processed in parallel. All samples extracted 1 h after IA learning were validated using RT-PCR to establish the induction of the immediate early gene (IEG) c-fos in each region at 1 h after training relative to untrained conditions. This c-fos induction is a well-established and commonly used proxy of neural activation following learning (Tischmeyer and Grimm, 1999; Katche et al., 2010; Zhang et al., 2011) (Fig. 1).
Experimental design and workflow. 1) Rats were trained in IA and killed 1 h or 6 d after training. Untrained rats were killed at matched time points and used as controls. 2) Brains were snap frozen in isopentane on dry ice immediately after death. Bilateral dHC and PL cortex tissue was punched on a cryostat and RNA was isolated. 3) 500 ng of RNA from each brain region from individual animals (n = 15) was reverse-transcribed to perform RT-PCR to test for c-fos induction 1 h after IA training, as a proxy for learning-induced cellular activity. Data are shown as mean percentages ± s.e.m. of untrained control (U) mean values. 4) RNA samples for each behavioral condition were pooled into 3 groups of 5 rats (3 biological replicates of n = 5). All replicates underwent RNA sequencing and bioinformatics analysis (alignment, read counting, normalization). 5) GO pathway analysis as well as cell type- and compartment-specific transcript analyses were conducted.
We conducted genome-wide RNAseq on three independent pools of bilateral dHC and PL cortex (5 rats per pool) in three behavioral conditions: untrained rats, trained rats killed at 1 h after IA training, and trained rats killed at 6 d after IA training (Fig. 1). The time points of 1 h and 6 d after training were chosen based on previous studies of contextual and spatial types of memories, including IA, indicating that gene expression following learning may persist for days. Indeed, significant transcriptional changes, such as induction of IEGs, are detected between 1 and 2 h after training (Tischmeyer and Grimm, 1999; Cavallaro et al., 2002; Levenson et al., 2004; Mei et al., 2005; Katche et al., 2010; Zhang et al., 2011; Bero et al., 2014; Poplawski et al., 2014; Peixoto et al., 2015; Cho et al., 2016; Halder et al., 2016). Furthermore, extended time course assessments following IA training in rats indicated that the dHC induction of specific genes required for memory consolidation, including CCAAT enhancer-binding protein β and δ, insulin-like growth factor 2, muscle-specific tyrosine kinase receptor, and BDNF persist for ∼1 d after training, and return to baseline by 2 d after training (Garcia-Osta et al., 2006; Bekinschtein et al., 2008; D. Y. Chen et al., 2011; Arguello et al., 2013; Bambah-Mukku et al., 2014). Finally, in the PL cortex, where much less has been explored about the temporal progression of learning-induced gene expression, BDNF protein levels were found to be significantly elevated 1 week after IA training (Ye et al., 2017).
dHC and PL cortex have partially distinct transcriptomic profiles
A comparison of the transcriptomic profiles of the dHC and PL cortex at baseline revealed that, out of a reference library comprised of 33,713 transcripts (Extended Data Fig. 2-1), 3062 were significantly enriched in the dHC, whereas 3157 were significantly enriched in the PL cortex (FDR < 0.1), and 8153 transcripts were similarly expressed in both regions (FDR > 0.1). Finally, 19,329 transcripts did not reach a threshold of 10 normalized read counts after averaging across biological replicates in either region. Therefore, these transcripts were considered not expressed (Fig. 2A,B).
Gene expression and GO analyses of dHC versus PL cortex. A, Venn diagram of transcript enrichment in the dHC versus PL cortex. Full dataset from differential gene expression analysis found in Extended Data Figure 2-1. B, Heatmap of 14,372 transcripts expressed in the dHC and PL cortex above threshold of 10 averaged normalized read counts. Data are expressed as –log10(read count). C, Transcript validation: 10 genes from the pool of transcripts selectively enriched in either the dHC or the PL cortex underwent single-gene validation via RT-PCR. RT-PCR validations were conducted on a new cohort of animals. Primer sequences used for RT-PCRs are listed in Extended Data Figure 2-5. n = 6/group. Data are shown as mean percentages ± s.e.m. of PL cortex untrained control mean values. *p < 0.05; **p < 0.01; ***p < 0.001; two-tailed Student's t test. Statistical values are shown in Extended Data Figure 2-6. D, GO biological processes analysis (top 5 most statistically significant processes in each condition shown) was performed for the transcripts enriched in dHC or PL cortex, as well as transcripts whose expression was similar in both regions. The size of the circle represents the number of genes comprising each GO biological process. Color represents statistical significance, indicated by the legend on the right (–log10 adjusted p value). Full dataset from GO enrichment analysis found in Extended Data Figure 2-2. E, Number of enriched transcripts per cell type in the DHC and PL cortex, as determined from a database of cell type-specific gene expression generated from a large-scale single-cell RNAseq conducted in the mouse cortex and HC (Zeisel et al., 2015). Lists of cell type-specific DEGs found in Extended Data Figures 2-3 and 2-4.
Figure 2-1
List of enriched transcripts in untrained rats in dHC versus PL cortex. Download Figure 2-1, XLSX file.
Figure 2-2
GO enrichment for dHC and PL cortex enriched and shared expression transcriptomic profiles in untrained rats. Download Figure 2-2, XLSX file.
Figure 2-3
Cell type-specific expression in untrained rats in dHC enriched transcripts, using gene lists from Zeisel et. al. (2015). Download Figure 2-3, XLSX file.
Figure 2-4
Cell type-specific expression in untrained rats in PL enriched transcripts, using gene lists from Zeisel et. al. (2015). Download Figure 2-4, XLSX file.
Figure 2-5
Primer sequences used for RT-PCR validation 1 h after IA training and in unpaired conditions. Download Figure 2-5, DOCX file.
Figure 2-6
Statistical values for validations shown in Figure 2C for regional enrichment in untrained rats. Download Figure 2-6, DOCX file.
Single-gene RT-PCR analysis on independent sets of dHC and PL cortex samples validated the RNAseq data. Ten transcripts from the subsets of high, medium, and low expression levels were randomly selected and assessed by RT-PCR, which confirmed their differential expression between the two regions (Fig. 2C).
We then performed GO enrichment analyses using the R package clusterProfiler (Yu et al., 2012) to compare the biological processes derived from either the enriched or the commonly expressed transcripts in the dHC and PL cortex. We found that among the biologically distinct pathways enriched in the dHC were gliogenesis, axon development, and lipid modification. In contrast, the pathways enriched in the PL cortex included synaptic vesicle cycle and vesicle-mediated transport in synapse (Fig. 2D). Biological processes significantly enriched in both brain regions, but comprising different gene sets, included synapse organization, regulation of membrane potential, as well as vesicle localization. Finally, significantly enriched biological processes derived from transcripts with shared expression levels between the dHC and PL cortex included noncoding RNA metabolic processes, ribonucleoprotein complex biogenesis, and DNA repair (Fig. 2D; Extended Data Fig. 2-2).
In order to assess the cell type specificity of the enriched transcripts at baseline in the dHC and PL cortex, we compared our findings with a database generated from single-cell RNAseq conducted in the mouse hippocampal CA1 subregion and somatosensory cortex (Zeisel et al., 2015). We reasoned that this analysis would be informative despite the cross-species comparison because ∼94%-95% of all known protein-coding and regulatory elements are conserved between the rat, mouse, and human genomes (Rat Genome Sequencing Project Consortium, 2004). Therefore, we inferred that cell-specific transcriptomes are likely to be similar between the rat and mouse genomes. This analysis revealed that, compared with the PL cortex, the dHC had an approximately threefold increase in oligodendrocyte-enriched transcripts, and an approximately fivefold increase in astrocyte-enriched transcripts. Conversely, the PL cortex had a higher number of interneuron-enriched transcripts compared with the dHC (Fig. 2E; Extended Data Figs. 2-3, 2-4). Both regions expressed a similar number of transcripts enriched in ependymal, endothelial, and mural cells (pericytes and vascular smooth muscle cells), suggesting that their vasculature transcriptomic profiles are similar (Fig. 2E).
Differential transcriptome profiles induced in dHC and PL cortex 1 h after IA training. A, Venn diagram displaying the number of DEGs at 1 h after IA training in the dHC and PL cortex as well as the number of common DEGs regulated in both regions. The 34 transcripts that were commonly regulated in both regions are listed below the Venn diagram. Full datasets are found in Extended Data Figures 3-1 and 3-2. B, Heatmap of all DEGs in the dHC and PL cortex, expressed as log2 fold change from average untrained read counts. C, D, Volcano plots for DEGs in the dHC (C) and PL cortex (D), with a horizontal line at the –log10(FDR = 0.1) threshold, and a vertical line at log2 fold change = 0, delineating positive versus negative regulation. Blue represents negatively regulated genes. Red represents positively regulated genes.
Figure 3-1
List of DEGs at 1 h after IA training in the dHC. Download Figure 3-1, XLSX file.
Figure 3-2
List of DEGs at 1 h after IA training in the PL cortex. Download Figure 3-2, XLSX file.
Distinct transcriptomic profiles are regulated 1 h after learning in the dHC and PL cortex
We next determined the transcriptomic profiles of the dHC and PL cortex at 1 h after IA training. We found 135 DEGs in the dHC and 478 DEGs in the PL cortex of trained rats compared with the respective regions of untrained animals (FDR < 0.1). Only 34 DEGs were commonly regulated in both regions (Fig. 3A; Extended Data Figs. 3-1, 3-2). A greater proportion of transcripts were downregulated in the PL cortex (55% downregulated/45% upregulated DEGs) compared with the dHC (21% downregulated/79% upregulated DEGs) (Fig. 3B–D). Of the 34 DEGs commonly regulated, 22 transcripts significantly increased in both regions. These included 15 IEGs, such as arc, egr1, egr4, nr4a1, nr4a2, and nr4a3 (Fig. 3A). IEGs are rapidly and transiently regulated in response to cell activation and are required for activity-dependent processes, including synaptic plasticity and memory (Morgan and Curran, 1991; Alberini, 2009; Minatohara et al., 2015). Other noteworthy transcripts commonly upregulated by learning in both dHC and PL cortex were slc2a1, which encodes the glucose transporter GLUT1, and s100β. s100β encodes a calcium-binding protein primarily expressed by astrocytes and oligodendrocytes and is reported to be involved in LTP and long-term memory in the mouse HC (Nishiyama et al., 2002). The remaining 12 common DEGs were significantly downregulated. These included mag, an integral myelin-associated glycoprotein gene, fcho1, which encodes an endocytic adaptor protein, and adam11, which encodes for an A Disintegrin and Metalloprotease (ADAM) family protein member (Fig. 3A).
These transcriptomic data were validated by using RT-PCR on independent sets of individual rat RNA samples. Twenty-two dHC DEGs and 20 PL cortex DEGs ranging from low to high read count expression were assessed. In the dHC, these transcripts included the glucose transporters slc2a1 and slc2a2; the IEGs nr4a1, nr4a2, and egr1; the dual specificity phosphatases dusp1, dusp4, dusp5, and dusp6; the phosphatase regulatory subunits ppp1r12a, ppp1cb, ppp3r1, and ppp4r2; and the axon guidance regulator plxna3 (Fig. 4A). In the PL cortex, they included the IEG arc, the structural molecules cntn2, map1b, and reln; the amyloid precursor protein encoding gene app; and the axon guidance regulators plxna1, plxna3, plxnb3, sema3c, and sema4d (Fig. 4B). Twenty-one of 22 transcripts in dHC and 19 of 20 in the PL cortex confirmed the RNAseq results.
IA training-induced gene expression validation and their assessment in context-shock unpaired controls. A, B, Validation of transcripts induced by training. Single-gene RT-PCR was performed on 22 DEGs (FDR < 0.1) in the dHC (A) and 20 DEGs in the PL cortex (B). These RT-PCRs were conducted using new cohorts of rats (n = 4-11/group, individual data points overlayed on bar graphs). C, D, RT-PCRs for the same genes that underwent the validation in A and B were conducted in the dHC (C) and PL cortex (D) of rats that were exposed to unpaired context and shock (n = 5 or 6/group, individual data points overlayed on bar graphs). Primers used for all RT-PCRs are listed in Extended Data Figure 4-1. Data are shown as mean percentages ± s.e.m. of untrained control mean values. *p < 0.05; **p < 0.01; ***p < 0.001; two-tailed Student's t test. Statistical values are shown in Extended Data Figures 4-2, 4-3, 4-4, and 4-5.
Figure 4-1
Primer sequences used for RT-PCR validation 1 h after IA training and in unpaired conditions. Download Figure 4-1, DOCX file.
Figure 4-2
Statistical values for validations shown in Figure 4A for transcripts induced 1 h after IA training in the dHC. Download Figure 4-2, DOCX file.
Figure 4-3
Statistical values for validations shown in Figure 4B for transcripts induced 1 h after IA training in the PL cortex. Download Figure 4-3, DOCX file.
Figure 4-4
Statistical values for validations shown in Figure 4C for transcripts induced in unpaired conditions in the dHC. Download Figure 4-4, DOCX file.
Figure 4-5
Statistical values for validations shown in Figure 4D for transcripts induced in unpaired conditions in the PL cortex. Download Figure 4-5, DOCX file.
To assess whether the DEGs were related to associative learning as opposed to context or footshock exposures alone, we performed RT-PCRs of the dHC and PL cortex DEGs on mRNA extracted from rats that underwent an unpaired context-shock protocol. In this paradigm, the rats are exposed to the context and 1 h later to the footshock. This unpaired paradigm does not evoke IA memory (Garcia-Osta et al., 2006; D. Y. Chen et al., 2011). The RT-PCR analyses revealed that 19 of 22 genes in the dHC (Fig. 4C) and 19 of 20 genes in the PL cortex (Fig. 4D) were regulated only following training but not in the unpaired condition, suggesting that the majority of transcriptional regulations found by the RNAseq analyses accompany associative learning. Two transcripts, plxna3 in the dHC and app in the PL cortex, which were significantly downregulated at 1 h after IA training, were also found significantly downregulated with the unpaired condition protocol, indicating that some transcriptional changes are linked to either context or footshock experiences. Conversely, slc2a1 and slc1a3, which were upregulated 1 h after IA training in the dHC, were significantly decreased in the dHC in unpaired conditions, indicating that unpaired context and footshock experiences are also accompanied by distinctive transcriptional regulations. In sum, our validations suggested that the transcriptomic profile at 1 h after IA training is mostly linked to associative context-footshock learning.
We then used R/clusterProfiler to perform GO enrichment analysis on the transcriptomic profiles to identify the biological processes regulated by learning. The most significantly regulated pathways in the dHC were regulation of cell projection organization, axon development, protein dephosphorylation, and glycogen and glucan biosynthetic processes (Fig. 5A; Extended Data Fig. 5-1). In the PL cortex, the pathways that emerged as most significantly regulated were related to axon development, guidance, and axonogenesis (Fig. 5B; Extended Data Fig. 5-1). Although these pathways were also significantly regulated in the dHC, the number of transcripts per pathway and the degree of significance of these pathways were greater in the PL cortex transcriptional profile. Protein dephosphorylation and metabolic process transcriptional programs were uniquely regulated in the dHC.
Most significantly regulated biological processes in the dHC and PL cortex 1 h after IA training. A, B, Most significant biological processes obtained from GO enrichment analysis of DEGs regulated in the dHC (A) and PL cortex (B) 1 h after IA training. Gene count indicated by size of dot, and statistical significance of each biological pathway indicated by color of dot (–log10 adjusted p value) as per the legend shown on the right. Full dataset from GO enrichment analysis found in Extended Data Figure 5-1.
Figure 5-1
GO enrichment analysis of dHC and PL cortex DEGs at 1 h after IA training. Download Figure 5-1, XLSX file.
Notably, GO analyses do not account for the directionality of the changes in the biological pathway; hence, to consider the contribution of increase versus decrease in transcription, we visualized the DEGs of the top 5 most significantly regulated biological processes from each region in gene-concept networks, which depict linkages between gene expression and their associated GO terms. In this visualization, the relative coordinates of the GO terms are determined by the number and identity of the underlying DEGs. Therefore, clusters of biological processes will form among GO terms that have many overlapping DEGs. Based on these parameters, two groups of functionally clustered transcripts emerged in the dHC: (1) regulation of cell projection organization and axon development; and (2) glycogen/glucan biosynthetic processes, and protein dephosphorylation (Fig. 6A). In the PL cortex, the most significant GO terms clustered tightly together and were all related to axon guidance and development. We noted that, in the PL cortex, the vast majority of transcripts related to axon guidance and development were downregulated at 1 h after training (Fig. 6B).
Gene expression networks of the most significant GO biological processes induced 1 h after training. A, B, Gene concept networks depicting the gene composition of the top 5 most significantly regulated biological processes determined from GO enrichment analysis of DEGs in the dHC (A) and PL (B) 1 h after IA training. The relative coordinates of the GO terms are determined by the interconnectedness of their underlying DEG sets; therefore, biological processes with largely overlapping DEGs will be in closer proximity to one another. The color of the gene's node relative to untrained animals (log2 fold change) represents expression levels of individual transcripts.
One limitation of these GO analyses is that they evaluate only transcripts that are selected based on the imposed statistical cutoffs. These types of analyses may underestimate or omit subtle, nonstatistically significant transcriptional regulations across biological pathways, which, however, may be functionally relevant. Hence, we next performed GSEA (Subramanian et al., 2005) to determine whether gene sets defined by GO biological processes show statistical differences between untrained and trained transcriptomic profiles. GSEA in the dHC 1 h after IA revealed that the most significantly downregulated gene sets were amino acid activation and regulation of translational fidelity (Fig. 7A; Extended Data Fig. 7-1) and the most significantly upregulated gene sets were negative regulation of ERK1 and ERK2 cascades and cotranslational protein targeting to membrane (Fig. 7B; Extended Data Fig. 7-2). Parallel GSEA in the PL cortex 1 h after IA revealed that the most significantly downregulated gene sets were the semaphorin-plexin signaling pathway and regulation of axon guidance (Fig. 7C; Extended Data Fig. 7-3) and the most significantly upregulated gene sets were establishment of protein localization to endoplasmic reticulum and cotranslational protein targeting to membrane (Fig. 7D; Extended Data Fig. 7-4).
Most significantly regulated gene sets according to GSEA in dHC and PL cortex 1 h after IA training. GSEA in the dHC (A,B) and PL cortex (B,C) at 1 h after IA training compared with untrained conditions. Most significant GO biological process gene sets downregulated (A,C) and upregulated (B,D) are depicted in bar graphs on the top of each panel, with graphs for the two most significantly regulated gene sets and their respective normalized enrichment score (NES) and adjusted p values (Adj. p) depicted below. A positive NES indicates upregulation at 1 h after training compared with untrained controls, and a negative NES indicates downregulation. Green lines indicate the enrichment profiles by visualizing the running enrichment score, which reflects the degree to which a gene set is overrepresented at the top or bottom of the list of genes ranked by log2 fold change. Vertical black lines indicate hits in the transcriptomic profile in the gene sets to visualize their location in the ranked ordered dataset, depicted in gray. Heatmaps are overlaid on top of the hits to highlight the upregulation/downregulation of the transcript compared with untrained animals. Full GSEA datasets are found in Extended Data Figure 7-1 (downregulated), Extended Data Figure 7-2 (upregulated) for the dHC, and Extended Data Figure 7-3 (downregulated), and Extended Data Figure 7-4 (upregulated) for the PL cortex.
Figure 7-1
Downregulation of dHC gene sets from GSEA of transcriptional profiles at 1 h after IA training compared with untrained conditions. Download Figure 7-1, XLSX file.
Figure 7-2
Upregulation of dHC gene sets from GSEA of transcriptional profiles at 1 h after IA training compared with untrained conditions. Download Figure 7-2, XLSX file.
Figure 7-3
Downregulation of PL cortex gene sets from GSEA of transcriptional profiles at 1 h after IA training compared with untrained conditions. Download Figure 7-3, XLSX file.
Figure 7-4
Upregulation of PL cortex gene sets from GSEA of transcriptional profiles at 1 h after IA training compared with untrained conditions. Download Figure 7-4, XLSX file.
Collectively, these data led us to conclude that the transcriptional programs induced at 1 h after IA training in the dHC and PL cortex were highly distinct.
Learning-induced DEGs in specific cell types and in dendritic and synaptic compartments
We then determined cell type specificity of DEGs regulated 1 h after IA training in the dHC and PL cortex by analyzing our data against a cell type-specific RNAseq repository generated by cell sorting in mouse HC and cortex (Zeisel et al., 2015). Although this is again a cross-species comparison, we believed that this analysis could be informative given that the vast majority of protein-coding and regulatory elements are functionally conserved between rat and mouse (Rat Genome Sequencing Project Consortium, 2004).
These analyses revealed that, in the dHC at 1 h after training, there is an upregulation of genes enriched in astrocytes, ependymal cells, endothelial cells, interneurons, microglia, and mural cells and a downregulation of genes enriched in oligodendrocytes (Fig. 8A; Extended Data Fig. 8-1). In the PL cortex at 1 h after training, genes enriched in astrocytes, endothelial cells, microglia, and mural cells were upregulated, whereas genes enriched in ependymal cells, oligodendrocytes, and interneurons were significantly downregulated (Fig. 8B; Extended Data Fig. 8-2).
Cell type- and compartment-specific expression of DEGs 1 h after IA training in the dHC and PL cortex. A, B, Percent of DEGs regulated 1 h after training in the dHC (A) and PL cortex (B) in cell type-specific populations obtained from the database in Zeisel et al. (2015). Populations included astrocytes, ependymal cells, endothelial cells, interneurons, microglia, mural cells, and oligodendrocytes. A full list of cell type-specific DEGs is found in Extended Data Figures 8-1 and 8-2. C, Pie charts represent the percent of dendritically and nondendritically localized DEGs found at 1 h after training, according to the database of Middleton et al. (2019). A list of DEGs expressed in dendrites is found in Extended Data Figures 8-3 and 8-4. D, Pie charts represent the percent of DEGs at 1 h after IA training in the dHC and PL cortex that relate to synaptic functions according to the SynGO database published in Koopmans et al. (2019). A list of DEGs related to synaptic function and full SynGO analysis dataset is found in Extended Data Figures 8-5 and 8-6.
Figure 8-1
List of cell type-specific DEGs in the dHC at 1 h after IA training according to the gene list from Zeisel et al. (2015). Download Figure 8-1, XLSX file.
Figure 8-2
List of cell type-specific DEGs in the PL cortex at 1 h after IA training according to the gene list from Zeisel et al. (2015). Download Figure 8-2, XLSX file.
Figure 8-3
List of DEGs in the dHC localized in dendrites, according to the gene list from Middleton et al. (2019). Download Figure 8-3, XLSX file.
Figure 8-4
List of DEGs in the PL cortex localized in dendrites, according to the gene list from Middleton et al. (2019). Download Figure 8-4, XLSX file.
Figure 8-5
List of synaptic function-related DEGs in the dHC at 1 h after IA training according to the SynGO database (Koopmans et al., 2019). Download Figure 8-5, XLSX file.
Figure 8-6
List of synaptic function-related DEGs in the PL cortex at 1 h after IA training according to the SynGO database (Koopmans et al., 2019). Download Figure 8-6, XLSX file.
To determine which DEG transcripts are expressed in dendritic compartments, we used a recently published comprehensive catalog of dendritically localized mRNAs generated from 8 high-throughput transcriptional studies (Middleton et al., 2019). These analyses revealed that 49.6% (62 of 125) of dHC DEGs and 45.4% (217 of 478) of PL cortex DEGs 1 h after training are localized to dendrites (Fig. 8C; Extended Data Figs. 8-3, 8-4).
Finally, we compared the DEGs of both the dHC and PL cortex 1 h after training with the SynGO database, a synapse gene and protein repository exclusively based on published, expert-curated evidence. This database accumulates research about synapse biology using GO annotations and novel ontology terms (Koopmans et al., 2019). We found that only 11.9% (16 of 125) of dHC DEGs versus 21.1% (101 of 478) of PL cortex DEGs were related to synaptic function (Fig. 8D; Extended Data Figs. 8-5, 8-6). These data led us to conclude that the PL cortex more significantly regulated transcripts encoding for synaptic proteins.
Transcriptomic profiles in dHC and PL cortex return to baseline by 6 d
The DEGs significantly regulated 1 h after training returned to baseline levels 6 d after IA training (Fig. 9A,B; Extended Data Figs. 9-1, 9-2). Indeed, no transcripts were found to be significantly different at 6 d after training compared with untrained controls in either brain region. RT-PCR analyses of independent, individual RNA samples extracted from the dHC and PL cortex 6 d after IA training confirmed these results (Fig. 9C,D). One exception was plxna3, which was found significantly downregulated in the dHC at 6 d after IA training. This gene encodes for a Class 3 semaphorin receptor that regulates the development of hippocampal axonal projections in mice (Cheng et al., 2001), and likely contributes to cytoskeletal remodeling via semaphorin-plexin signal transduction pathways (Goshima et al., 2016). Studies of a similar gene in zebrafish suggested that this protein is important for axon pathfinding in the developing nervous system (Tanaka et al., 2007).
Training-induced mRNA changes return to baseline 6 d after training. A, B, Heatmaps of DEGs in the dHC (A) and PL cortex (B) at 1 h and 6 d after IA training with log2 fold change of each gene compared with the average normalized read count of untrained animals. All transcripts that were differentially regulated at 1 h after IA training returned to baseline levels 6 d after IA training. No significant changes were found at 6 d after training compared with the untrained group throughout the transcriptome (using a threshold of FDR < 0.1). Full dataset from differential gene expression analysis found in Extended Data Figures 9-1 and 9-2. C, D, Ten genes used for validation at 1 h after IA training (see Fig. 1A,C) underwent single-gene RT-PCRs on RNAs obtained from dHC (C) and PL (D) cortex extracted 6 d after training. RT-PCR validations were conducted on a new cohort of rats. n = 5 or 6/group. Data are shown as mean percentages ± s.e.m. of untrained control mean values. *p < 0.05 (two-tailed Student's t test). Statistical values are shown in Extended Data Figures 9-3 and 9-4.
Figure 9-1
List of DEGs at 6 d after IA training in the dHC. Download Figure 9-1, XLSX file.
Figure 9-2
List of DEGs at 6 d after IA training in the PL cortex. Download Figure 9-2, XLSX file.
Figure 9-3
Statistical values for validations shown in Figure 9C for transcripts induced 6 d after IA training in the dHC. Download Figure 9-3, DOCX file.
Figure 9-4
Statistical values for validations shown in Figure 9D for transcripts induced 6 d after IA training in the PL cortex. Download Figure 9-4, DOCX file.
Together, these data suggested that the learning-induced transcriptomic profiles of the dHC and PL cortex regulated at 1 h after training return to baseline 6 d later.
Discussion
Obtaining transcriptomic profiles of multiple brain regions in response to learning is key for understanding the systems-level biology underlying long-term memory formation. In this study, we showed that the baseline transcriptional profiles of the rat dHC and PL cortex, two regions essential for the formation of episodic contextual memories, are only partially overlapping, as approximately one-fifth of their transcripts were found to be selectively enriched in each region. In addition, we found that the transcriptional profiles regulated in each region at 1 h after IA learning are distinct, and, finally, that these changes return to baseline by 6 d after training.
At baseline, out of a reference library of 33,713 transcripts, 42.6% of transcripts were detected in the dHC and PL cortex, a size similar to what was found previously by the Allen Brain Institute in the adult mouse brain using high-throughput ISH (Lein et al., 2007) (http://www.brain-map.org). GO analyses of our data revealed that the biological processes most significantly enriched in the dHC include gliogenesis, axon development, and lipid modification, a result also confirmed by our additional analyses using the cell type-specific database generated by single-cell RNAseq in mouse brain (Zeisel et al., 2015). The biological processes most significantly enriched in the PL cortex were related to synaptic function and synaptic vesicle localization, whereas those significantly enriched in both brain regions included general brain cellular mechanisms, such as synapse organization, regulation of membrane potential, and vesicle localization. In sum, at baseline, while the HC expresses more transcripts related to glia and axon mechanisms, the PL cortex expresses more transcripts related to synaptic mechanisms. Notably, GO biological pathways that were commonly expressed in both regions comprised different gene sets, suggesting that a much higher biological distinction exists at the transcriptional level between dHC and PL cortex. Finally, the transcripts expressed in both regions reflected conserved, critical cellular mechanisms, such as noncoding RNA metabolic processes, ribonucleoprotein complex biogenesis, and DNA repair.
Our analyses were conducted on whole dHC and PL cortex extracts; therefore, they did not directly provide information about transcriptomic distinctions of the numerous cell populations present in the two regions. Our rationale for using whole-tissue mRNA extracts was to avoid the exclusion of important mRNAs that are localized in cellular processes, such as neuronal axons and dendrites, which is a confound in most studies that aim at obtaining cell type-specific sequencing analyses using single-cell or cell population separation procedures (van den Brink et al., 2017; Nguyen et al., 2018). Because one limitation of RNAseq is that it requires read count cutoffs, we cannot exclude that more (or less) overlap in the transcriptional profiles between brain regions exists if very low expressing transcripts could be included.
The transcriptomic profiles regulated in response to learning in the dHC and PL cortex were mostly distinct. First, only 34 of 613 total DEGs were shared. Second, the PL cortex regulated 4 times as many transcripts compared with the dHC. These remarkable differences remain to be understood and should be explored in future studies; here we speculate that, between these regions, the learning-evoked changes taking place in the PL cortex may involve distinct and/or a higher number of cell types. Future studies should also determine whether the learning-dependent gene expression occurs in discrete populations of cells that already express those transcripts, or whether it recruits new cells. It is likely that both of these types of regulation may occur based on the nature of the transcript, cell type, and brain area.
Our data showed that some biological pathways, such as axon development and axonogenesis, were more significantly regulated in the PL cortex 1 h after training compared with the dHC. These results imply that mechanisms different from the classically studied postsynaptic changes and dendritic remodeling (Leuner et al., 2003; Segal, 2017) are rapidly regulated on learning. More specifically, we found that learning induces significant downregulation in the PL cortex of the following: (1) Semaphorin-Plexin signaling via sema4d, sema3c, plxna1, plxna3, and plxnb3; (2) Slit-Robo signaling via dpysl2, dpysl4, dpysl5, and slit1; and (3) Netrin-1-DCC signaling via regulation of DCC netrin-1 receptor (dcc), all of which are pathways critical for axon guidance and axonogenesis. GSEA confirmed these results, demonstrating that the semaphorin-plexin signaling pathway and regulation of axon guidance pathways are downregulated. These data are in line with recent studies demonstrating that axon guidance pathways, which were primarily studied in synaptogenesis in the developing CNS (Tessier-Lavigne and Goodman, 1996), are also expressed in the mature brain and functionally involved in learning-induced brain plasticity in adult systems (Horn et al., 2013; Wong et al., 2019).
Conversely, the GO terms most significantly changed in the dHC, but not in the PL cortex, were related to metabolic processes, including glycogen and glucan biosynthesis via transcriptional upregulation of ptges3, epm2aip1, per2, and dyrk2. These results are in line with previous data indicating a key role for glycogen metabolism in the dHC following IA learning to support the high metabolic demands of active neuronal networks (Suzuki et al., 2011; Descalzi et al., 2019). To our knowledge, thus far, no studies have demonstrated a requirement of glycogen metabolism for long-term memory consolidation in the adult PL cortex. Hence, future investigations could determine whether these mechanisms may represent region-specific distinctive biological signatures. Another biological process found highly regulated in the dHC but not the PL cortex 1 h after training is protein dephosphorylation. This process included the upregulation of dual specificity phosphatases (DUSPs), including dusp1, dusp4, dusp5, and dusp6, which are Class I classical cysteine-based protein phosphatases that have the dual ability to dephosphorylate phospho-serine/threonine and phospho-tyrosine residues, and have been implicated in hippocampal memory as well as neurodegeneration (Abdul Rahman et al., 2016; Bhore et al., 2017). The protein dephosphorylation process also included several members of the protein serine/threonine phosphatase family, both catalytic and regulatory subunits (ppp1cb, ppp3r1, ppp1r12a, ppp4r2), which have a broad range of functions in the brain, including synaptic transmission, LTP, and hippocampal-dependent memories (Mansuy and Shenolikar, 2006).
It should be noted that one limitation of GO analysis is that the pathways are not independently defined; indeed, the same transcripts are included in multiple GO biological processes, which are functionally linked. For example, 5 DEGs upregulated in the dHC (dyrk2, epm2aip1, per2, ppp1cb, ptges3) were classified in both glucan and glycogen biosynthetic processes. Similarly, transcript overlap was found in axon development and axonogenesis, or also in axon guidance and neuron projection guidance pathways. Thus, it is important to examine the data using additional methodologies. We supplemented the GO analyses with GSEA and additional cell type-specific comparisons, which confirmed the significant downregulation of axon-guidance and semaphorin-plexin pathways in the PL cortex and the regulation of translation pathways in the dHC.
We also found that a greater proportion of transcripts were downregulated 1 h after training in the PL (55% downregulated/45% upregulated) compared with the dHC (21% downregulated/79% upregulated). The reason for these decreases remains to be understood in future studies. This downregulation may result from decreased transcription because of repression, or decreased mRNA stability. As indicated by previous studies, transcriptional repressors, such as MeCP2, Sin3a, HDAC2, DNA methyltransferases, and protein kinase GCN2, are induced by learning and play an active role in forming long-term memory (Costa-Mattioli et al., 2005; Miller and Sweatt, 2007; Guan et al., 2009; Bambah-Mukku et al., 2014). Transcriptional downregulation following contextual learning in the HC has also been reported using microarray or RNAseq at multiple early time points after training (Levenson et al., 2004; Cho et al., 2016). For example, Levenson et al. (2004) showed that hippocampal transcripts were mostly upregulated at 1 h after contextual fear conditioning but were mostly downregulated starting at 2 h and extending to 4 h after training in mice. Nevertheless, we should keep in mind that the differential expression patterns we found in the PL cortex and dHC are snapshots of dynamic processes of transcription regulation induced by learning. A complete understanding of the biological changes underlying memory consolidation should include omic analyses conducted on numerous time points following training.
Of note was the finding that oligodendrocyte-enriched transcripts, as well as transcripts associated with axonogenesis and axon guidance, were significantly downregulated 1 h after training, especially in the PL cortex. These data are in disagreement with recent studies reporting that spatial or schema learning leads to an increase in oligodendrogenesis in the adult PFC in the context of experience-dependent myelin remodeling in cortical regions underlying long-term memory (Hasan et al., 2019; Steadman et al., 2020). A possible explanation for this discrepancy is that these studies assessed myelin proliferation on the scale of days to weeks after training and did not explore the transcriptional dynamics that occur immediately after training. Perhaps myelin-based mechanisms first need to decrease to allow a subsequent increase in oligodendrocyte changes that would be needed to stabilize newly established connectivity. In line with this idea, a decrease in myelin proteins was reported by proteomic analysis of the dentate gyrus 24 h after contextual fear conditioning in mice (Houyoux et al., 2017). If our hypothesis is correct, it is possible that a delayed wave of oligodendrocyte transcript upregulation in the PL cortex may take place to support systems consolidation. It is also possible that myelination and axonal mechanisms, both in terms of temporal progression and types of molecules involved, may vary greatly among different brain areas, including distinct cortical subregions. For example, Bero et al. (2014) found an increase in axon guidance GO biological processes in the ACC (a subregion that is adjacent to the PL cortex and involved in long-lasting memory storage) 1 h after contextual fear conditioning in mice following RNAseq analysis.
One other intriguing observation that emerged from our data was that the PL cortex showed a significant downregulation of transcripts specific for interneuron markers at 1 h after IA training. One of these was reln, which encodes a secreted extracellular matrix glycoprotein that is expressed predominantly in inhibitory neurons in the adult brain and plays a prominent role in synaptic function by modulating calcium currents in NMDARs (Y. Chen et al., 2005; Stranahan et al., 2013; Pohlkamp et al., 2014). Another transcript that was significantly decreased in the PL cortex after training was rab3b, which encodes a synaptic vesicle protein enriched at hippocampal inhibitory synapses and required for LTD (Tsetsenis et al., 2011). These data led us to speculate that learning may trigger a rapid decrease in inhibition, hence boosting excitability in this region.
We found that ∼50% of DEGs regulated in the dHC and PL cortex 1 h after training are localized in the dendritic compartment (Middleton et al., 2019). Additionally, ∼21% of all DEGs in the PL cortex at 1 h after training encoded proteins related to synaptic function, compared with only ∼11% of all DEGs in the dHC. These data indicate that the dHC, and to a greater extent the PL cortex, rapidly regulates compartment-specific and synaptic plasticity transcripts on learning. This is in line with findings showing that both regions are rapidly activated by training (Holloway and McIntyre, 2011; Lesburgueres et al., 2011; Tse et al., 2011; Bero et al., 2014; Halder et al., 2016; Rizzo et al., 2017).
All transcriptional changes returned to baseline by 6 d after IA training in both regions, a time point at which IA memory retention is still highly significant, as it persists for several weeks after training (D. Y. Chen et al., 2011; Inda et al., 2011; Ye et al., 2017). These data indicated that, despite the fact that cortical regions, such as the PFC and ACC, undergo systems consolidation for several weeks in rodents (Frankland and Bontempi, 2005; Squire et al., 2015; Kitamura et al., 2017), the transcriptomic profiles induced by learning do not remain significantly changed 6 d after training. However, it cannot be excluded that post-transcriptional changes are still engaged. This was, for example, documented in previous studies from our laboratory in which a sustained upregulation in BDNF and synaptic adhesion molecules neuroligin 1 and neuroligin 2 proteins were found at 1 week after IA training in the PL cortex, despite their mRNA levels remaining unchanged (Ye et al., 2017).
Future studies should assess whether our RNAseq data are generalizable to other learning tasks and species by using different aversive or nonaversive paradigms. These studies should also determine whether the distinct transcriptional profiles of the dHC and PL cortex regulate distinct or convergent biological functions.
In conclusion, our systems-level biology approach indicates that the biological composition of the different brain regions is significantly distinct. Furthermore, the mRNA changes taking place in different brain regions to support memory consolidation are mostly divergent. Our data also identified numerous novel pathways as potentially implicated in long-term memory formation. We conclude that a comprehensive systems-level biology approach is required for understanding the biology of long-term memory.
Footnotes
This work was supported by National Institutes of Health Grant MH065635 to C.M.A. A.K. was supported by National Institutes of Mental Health Grant F31MH116585. The RNAseq analysis was performed using High Performance Computing resources at New York University Langone Health and library preparation and sequencing through the Genome Technology Center. This shared resource is supported in part by Laura and Isaac Perlmutter Cancer Center Cancer Center Support Grant P30CA016087.
The authors declare no competing financial interests.
- Correspondence should be addressed to Cristina M. Alberini at ca60{at}nyu.edu
