Persistent transcriptional programmes are associated with remote memory

Abstract

The role of gene expression during learning and in short-term memories has been studied extensively^1-3, but less is known about remote memories, which can persist for a lifetime⁴. Here we used long-term contextual fear memory as a paradigm to probe the single-cell gene expression landscape that underlies remote memory storage in the medial prefrontal cortex. We found persistent activity-specific transcriptional alterations in diverse populations of neurons that lasted for weeks after fear learning. Out of a vast plasticity-coding space, we identified genes associated with membrane fusion that could have important roles in the maintenance of remote memory. Unexpectedly, astrocytes and microglia also acquired persistent gene expression signatures that were associated with remote memory, suggesting that they actively contribute to memory circuits. The discovery of gene expression programmes associated with remote memory engrams adds an important dimension of activity-dependent cellular states to existing brain taxonomy atlases and sheds light on the elusive mechanisms of remote memory storage.

Extended Data Fig. 1 |. Fidelity of the TRAP2; Ai14 line and sequencing quality metrics.

a, Fidelity with which the tdTomato reporter in TRAP2; Ai14 mice captures endogenous cFos expression during fear memory encoding by in situ hybridization of cFos and iCre in the mPFC directly following fear conditioning. (Left) Percentage of cFos⁺ cells in either fear conditioned (FC) or homecage (HC) mice as seen by in situ staining. Percentage of cFos⁺ cells that are also iCre⁺. (Right) Representative in situ hybridization images of FC mice (mPFC). b, Representative images of regions of the mPFC analysed in this study (anterior cingulate cortex (ACC), the prelimbic (PrL) and the infralimbic cortex (IL)). Scale bars are 1mm (and 0.5 mm in the insets). c, Violin plots of the number of reads and number of genes per biological replicates (m1-m4), per cell. d, Scatterplot depicting the lack of strong relationship between number of genes detected and number of reads obtained per cell. e, Scaled expression of canonical markers in non-neuronal cells (Cldn5-BECs, Pdgfra-OPCs, Cx3cr1-microglia, Aqp4-astrocytes)

Extended Data Fig. 2 |. Distribution of cell numbers and neuronal subtypes across various training conditions.

a, Representative flow cytometry plot of the amount of tdT⁺ events per training condition. In scatter plot, each point represents one mouse. b, Number of cells from each biological replicate that were annotated as one of 7 defined neuron subtypes (C0-C6). c, Number of TRAPed and Inactive cells (as defined by non-zero expression of tdT mRNA) collected per neuron subtype, in either fear-recall (FR) or no-fear (NF) mice. d, Representative images of subtype marker genes (C0 to C6) in the Allen Brain ISH atlas. e, Number of TRAPed (tdT mRNA+) cells collected in each experimental condition that fall in one of 7 neuronal subtype categories.

Extended Data Fig. 3 |. Differential gene expression in distinct neuronal subtypes (FR over NF TRAPed populations).

a, DEGs in each neuron subtype when inactive (tdT^-) neurons are compared between FR and NF mice. b, (Left) Volcano plots of DEGs in FR vs NF mice for each neuron subtype (C0 to C4). DEGs found when all replicates are pooled analysed in a combined manner are shown in red. Recall-dependent DEGs (defined as being differentially expressed in at least ¾ replicates, when analysed replicates are analysed individually) are labelled in black. (Right) Permutations are performed for every recall-dependent DEG for each neuronal population. Upregulated DEGs lying above the y = x line (red) and downregulated DEGs lying below the y = x line (blue) are considered to be above the 99^th percentile of the permuted distribution. c, GO enrichment analysis of all up- and downregulated DEGs (941 DEGs up, 384 DEGs down, all neuron subtypes combined) when all replicates are pooled. Bars show the enrichment scores (GeneAnalytics) for the GO pathway and dot indicates the number of DEGs involved in that pathway. d, GO enrichment analysis of only the upregulated remote-memory-specific DEGs (from Fig. 3c) from all neuron subtypes combined. Bar indicates enrichment score (Gene Analytics) and number indicates number of recall-dependent DEGs involved in that pathway.

Extended Data Fig. 4 |. Analysis of TRAPed ensembles in food salience (S) versus no salience (NS) mice.

a, Schematic of experimental paradigm for generating TRAPed neuronal ensembles as a result of food salience (food deprivation followed by food return (salience) or no food return (no salience)). b, Percentage of events in flow cytometry that were tdT⁺, by experience. c, tSNE of the merge of data from fear-recall experiments and food-salience experiments colored by neuron subtype (left) and experimental paradigm (right). d, Subtype composition differences between TRAPed fear-recalled ensembles and TRAPed food salience ensembles, as compared to background Inactive ensembles. e, Heatmap of the average log2 fold change of DEGs in each neuron subtype when comparing fear-recall vs no-fear, and salience vs no salience. Only DEGs with FDR <0.01 are shown.

Extended Data Fig. 5 |. DEGs when comparing ensembles from food salience (S) to no salience (NS) mice.

Volcano plots showing the log2 fold change and adjusted p-values (in log10 scale) of genes when comparing food salience over no salience groups within each neuron subtype (C0 through C4). Top DEGs per neuron subtype are labelled in red. Positive log2 fold change indicates upregulation in food salience (S) group.

Extended Data Fig. 6 |. Neuron subtype-specific activation programs, hypothesized protein–protein interactions and upstream regulatory motifs.

a, Fraction of cells in each neuron subtype that are induced with the transcriptional program (that is, DEGs) from a neuron subtype. Overall, the activation program of each TRAPed neuron subtype is found to be more specific to it than the inactive population, or other subtypes. b, (Left) De novo regulator motif discovery: analysis was performed using HOMER on the subset of 99 remote-memory-associated DEGs by looking at the sequences −400 to +100 bp from the TSS. 12 de novo and 2 known motifs were found (only motifs with an enrichment p-value <1e-2) were kept). Heatmap depicts the “motif score” of each DEG for each motif, and genes and motifs were clustered via the ward.D method. (Right) Bar graph depicting the % of the DEGs (target sequences) that possess a match for the motif within −400 to +100 bp from the TSS, vs the % of background sequences. For de novo motifs, the best match gene is listed on the right. HIF1b and HIF1a are matches to known motifs. c) (Left) Hypothesized protein–protein interactions of a subset of recall-dependent DEGs (TRAPed FR/NF) using the STRING database (https://string-db.org/). Only genes that are connected at a confidence level of 0.4 (medium) are shown. Connections indicate a possible existence of an interaction between two proteins. Genes are colored by up of downregulation in FR/NF. (Right) Same network plot, with nodes colored by the neuron subtype which differentially regulates the DEG.

Extended Data Fig. 7 |. In situ validation of tdT levels, neuronal subtype compositions and remote-memory-specific DEGs in the mPFC.

a, Ratio of Nuclei that are tdT⁺ (mRNA level) per training condition. Each data point represents one region of interest. b, Ratio of TRAPed cells that are positive for a neuronal subtype marker obtained either via the RNA-scope method, or by scRNA-seq (see Fig. 2). TRAPed cells are defined as DAPI⁺/tdT⁺ in RNAscope quantification, and as tdT mRNA count >1 in scRNA-seq (post-QC). No significant differences are found between FR and NF within either RNA-scope or scRNA-methods, indicating no major changes in neuronal subtype composition of active populations in different training conditions. c, in situ hybridization of key remote-memory specific DEGs including Stx1b in Rprm⁺/tdT⁺ cells, Syt13 in Tnfaip8l3⁺/tdT⁺ cells, Vamp2 in Tesc⁺/tdT⁺ cells. Scale bars = 100 μm.

Extended Data Fig. 8 |. DEGs and potential cell–cell interactions in non-neuronal cells during memory (re)consolidation.

a, Volcano plots of non-neuronal cell types when comparing cells in FR over NF nice. DEGs (FDR >0.01, log2FC >1) are labelled in red, and exemplary DEGs (high log2FC and log10FDR) are labelled in black. b, Number of non-neuronal cells collected in this study, for each cell type and experimental condition. c, Heatmap of a subset of neuronal ligands and glial receptors that are found to be differentially perturbed upon memory consolidation. Only receptors and ligands which were found to be (differentially) expressed are shown. d, (Left) Heatmap of the log2FC of DEGs (FR over NF) in neurons that are classified as ligands. (Middle and right) Sankey plot of known ligand-receptor pairs and heatmap of the average scaled expression level of the corresponding receptors in each cell type. e, (Left) Heatmap of the log2FC of DEGs (FR over NF) in neurons that are classified as receptors. (Middle and right) Sankey plot of known ligand-receptor pairs and heatmap of the average scaled expression level of the corresponding ligands in each cell type.

Extended Data Fig. 9 |. Comparison of remote-memory DEGs with previously published data sets of experience-dependent transcriptional activity.

a, (Left) Heatmap of the log2FC of all 1292 DEGs (FDR <0.05, FR over NF, all cells pooled) in this manuscript, and their log2FC values in previously published data sets of experience-dependent DEGs in activated neurons during: recent fear memory retrieval (Rao-Ruiz, 2019), associative fear-learning (Cho, 2016), post-visual stimulus (Hrvatin, 2018), or novel environment exposure (Lacar, 2016). A value of zero log2FC indicates the gene was not differentially expressed in a data set. (Right) Same as (left), but now DEGs are filtered down to the “Recall-dependent DEG ” set derived from this manuscript. Only genes differentially expressed in ¾ replicates are remaining. b, (Left) Log2 fold change heatmap of the recall-dependent DEGs between tdT⁺ vs tdT^- neurons in FR mice (genes are differentially expressed in >3/4 replicates) undergoing remote fear memory consolidation. (Right) The log2FC values of these genes if they are found in previously published data sets of experience-dependent DEGs (see (a)). A value of zero log2FC indicates the gene was not differentially expressed in that data set.

Extended Data Fig. 10 |. Comparison of remote-memory-specific DEGs and fear-experience-related DEGs.

(Left) Log2 fold change heatmap of the union of DEGs from comparing FR vs NF (remote-memory specific) and NR vs NF (fear-related). Sustained transcriptional changes from the fear-experience itself is shown in the yellow highlighted columns. (Right) Zoomed in view of the portion of the heatmap within the green boxes where most fear experience-related DEGs are located.

Fig. 1 |. Labelling and collection of single memory engram cells via the TRAP2; Ai14 line.

a, The experimental paradigm includes generating remote fear-memory traces via contextual fear conditioning, isolation of TRAP⁺-activated neurons via flow cytometry and plate-based single-cell RNA sequencing (scRNA-seq). b, Representative image of tdT⁺ TRAPed (red) cells in the anterior cingulate cortex (ACC) region of the mPFC 9 days after an injection of 4-OHT (at the time of remote memory recall). Scale bars, 1 mm (0.5 mm for zoomed in images). c, Degree of freezing upon return to the novel context 16 days after fear conditioning (FR) or no conditioning (NF) (n = 4 mice per condition, ***p= 4.89 × 10⁻⁹). d, Representative flow gating for tdT⁺ TRAPed cells (~1.5% of events in FR) in one mouse from the FR condition (left). Post-sequencing analysis of 722 example cells from a representative FR brain shows enrichment of tdT mRNA (pos) in the positive sort gate (orange violin) (middle). Enrichment of neuronal cells in the positive gate and neuronal and non-neuronal cell types in the negative gate in all brains (right). e, All sorted cells, with neuronal cells identified via the expression of Snap25 mRNA. All training and control conditions are represented in all cell clusters.

Fig. 2 |. Molecular identification of active neurons during (re)consolidation of a remote memory trace.

a, Dimensional reduction of all Snap25⁺ neurons (n = 3,691 cells) reveals seven distinct neuronal subtypes (C0–C6). b, Neuronal subtypes fall into two distinct categories: excitatory (glutamatergic) and inhibitory (GABAergic). The Glut–GABA signature is calculated based on the difference of the scaled expression level of Gad1 and Slc17a7. c, Expression levels of the top marker gene for each neuron subtype (C0–C6). d, Differences in the composition of neuron subtypes of TRAPed populations in FR and NF conditions, as well as inactive populations in FR mice (1 mouse per condition per replicate, n = 5 replicates, two-sided t-test) in the composition of TRAPed populations between FR and NF conditions are found. FR TRAPed populations are composed of significantly more C2 (GABAergic) neurons and less C3 (glutamatergic) neurons than the inactive population (error bars indicate SEM).

Fig. 3 |. Transcriptional programmes activated by (re)consolidation of remote memories are distinct per neuron subtype.

a, DEGs in FR versus NF in the C0 neuron subtype (n = 126 cells (NF), n = 289 cells (FR)) with all replicated pooled. Differential expression is defined by FDR < 0.01 and abs(log₂FC) > 0.3 (red points) via a two-sided Mann–Whitney test. Remote-memory-associated DEGs (that is, DEGs that remain differentially expressed in three-quarters or more of replicates and when FR is compared to the NR and HC groups (see Extended Data Fig. 3)) are labelled in black. b, The number of DEGs per neuron subtype (left) and the number of shared DEGs between each neuron subtype (using pooled cells) (right). c, log₂FC of remote-memory-associated DEGs (FR versus NF) per neuron subtype. Each gene is further annotated with potential functions. DE, differentially expressed; EXH, excitatory; GO, Gene Ontology; INH, inhibitory; NS, not significant. d, The percentile in which a TRAPed neuron (from C0) lies in the distribution of expression of a C0 DEG for all TRAPed C0 cells. The box plots (median ± s.d.) show the log₂cpm distribution for each C0 DEG. Hierarchical clustering reveals one common C0 transcriptional programme that is concertedly upregulated. e, The fraction of cells in each neuron subtype that is activated with the transcriptional programme (that is, DEGs) from C0.

Fig. 4 |. Remote memory (re)consolidation is associated with specific markers for vesicle exocytosis.

a, A subset of remote-memory-associated DEGs (FDR < 0.05 and log₂FC > 0.2) that are known to regulate vesicle membrane fusion and exocytosis at the presynaptic terminal. The violin plots are overlayed with a drumstick plot indicating the average expression per mouse. The red points represent the median. The bubble plots depict the log₂FC and the degree of significance (FDR) for each replicate and are coloured by such. In addition to the FR/NF log₂FC (first column), DEGs were also confirmed to be differentially expressed when compared to NR (second column), and their activation is specific only to the active (inactive (In)) (last column). The red dots indicate which neuronal subtype these particular genes are upregulated in. b, Representative in situ images of Serinc3 expression (purple) in Dkkl1⁺/tdT⁺/DAPI⁺ cells in the mPFC. Scale bars, 100 μm. c, In situ validation of key genes involved in vesicle exocytosis in various neuron subtypes, including Serinc3 (in the Dkkl1⁺ subtype, n = 268 (FR), n = 62 (NF) cells), Stx1b (in the Rprm⁺ subtype, n = 342 (FR), n = 144 (NF) cells), Syt13 (in the Tnfai8lp⁺ subtype, n = 244 (FR), n = 44 (NF) cells) and Vamp2 (in the Tesc⁺ subtype, n = 326 (FR), n = 292 (NF) cells). Each point represents the normalized integrated intensity of the probe per cell analysed.

Fig. 5 |. Transcriptomic changes in non-neuronal cells associated with remote memory (re)consolidation.

a, tSNE of all non-neuronal cells reveals five non-neuronal cell types (astrocyte, endothelial (EC), microglia, OPC and oligodendrocyte (oligo)) that were collected in an unbiased manner through tdT-negative flow cytometry gates. b, DEGs in non-neuronal cells (FR versus NF) (left). DEGs are defined as log₂FC > 1 and FDR < 0.01. The number of DEGs that satisfy these criteria in each non-neuronal cell type is also shown (right). The top DEGs (FR versus NF) for glial cells (astrocytes and microglia) that are also differentially expressed in FR versus NR are labelled. c, DEGs (determined by two-sided Mann–Whitney test) that are upregulated and downregulated in astrocytes (left) and microglia (right) in FR versus NF mice. DEGs (FDR < 0.01, log₂FC > 1) are indicated by red dots, and the top DEGs are labelled in black. d, Pathway analysis of the DEGs (FR versus NF) in microglia and astrocytes. The score is defined as the –log₂(P value) using the GeneAnalytics software.