Abstract
Experience-dependent gene expression reshapes neural circuits, permitting the learning of knowledge and skills. Most learning involves repetitive experiences during which neurons undergo multiple stages of functional and structural plasticity. Currently, the diversity of transcriptional responses underlying dynamic plasticity during repetition-based learning is poorly understood. To close this gap, we analyzed single-nucleus transcriptomes of L2/3 glutamatergic neurons of the primary motor cortex after 3 d motor skill training or home cage control in water-restricted male mice. “Train” and “control” neurons could be discriminated with high accuracy based on expression patterns of many genes, indicating that recent experience leaves a widespread transcriptional signature across L2/3 neurons. These discriminating genes exhibited divergent modes of coregulation, differentiating neurons into discrete clusters of transcriptional states. Several states showed gene expressions associated with activity-dependent plasticity. Some of these states were also prominent in the previously published reference, suggesting that they represent both spontaneous and task-related plasticity events. Markedly, however, two states were unique to our dataset. The first state, further enriched by motor training, showed gene expression suggestive of late-stage plasticity with repeated activation, which is suitable for expected emergent neuronal ensembles that stably retain motor learning. The second state, equally found in both train and control mice, showed elevated levels of metabolic pathways and norepinephrine sensitivity, suggesting a response to common experiences specific to our experimental conditions, such as water restriction or circadian rhythm. Together, we uncovered divergent transcriptional responses across L2/3 neurons, each potentially linked with distinct features of repetition-based motor learning such as plasticity, memory, and motivation.
- L2/3 glutamatergic neuron
- motor engram
- mouse primary motor cortex
- plasticity
- single-nucleus RNA sequencing
- water deprivation
Significance Statement
Learning is ultimately supported by coordinated gene expression changes in response to experience. Using single-nucleus RNA sequencing, we characterized divergent transcriptional states that L2/3 motor cortical neurons undergo during long-term motor learning. In particular, we discovered two salient transcriptional states that present signatures of motor engram-like reactivation and arousal-related metabolism, likely reflecting the unique behavioral states of animals in our study such as motor learning, arousal, and water deprivation. Our rich gene expression dataset will provide valuable information for elucidating electrophysiological and molecular mechanisms regulating multistage neuronal plasticity and developing interventions for behavioral disorders related to abnormal plasticity and learning.
Introduction
Salient experiences induce structural and functional changes in neurons, shaping how they process and transmit experience-related information, thereby supporting learning. For example, motor learning enhances the plasticity of dendritic spines (excitatory synapses) on neurons in primary motor cortex (M1; Xu et al., 2009; Yu and Zuo, 2011; Fu et al., 2012; Peters et al., 2014; Hedrick et al., 2022), corresponding to new, emergent patterns of neuronal output (Peters et al., 2014; Hwang et al., 2019). Suppressing such plasticity in M1 impairs learning (Hasan et al., 2013; Chen et al., 2015), supporting a deep connection between these phenomena.
Long-term maintenance of plasticity ultimately requires changes in gene expression, so neurons must encode the occurrence of plasticity as one or many transcriptional states. This can be observed in part at the level of activity-regulated genes (ARGs), the expression of which reflects recent neuronal activity (West et al., 2002; Cao et al., 2015; Yap and Greenberg, 2018). The primary ARGs, known as immediate early genes (IEGs), are rapidly transcribed in response to neuronal activity and facilitate both direct effects and secondary transcriptional responses that support diverse forms of plasticity. Since plasticity can be gated or shaped by numerous signals, including neuromodulators, the ARG network could potentially incorporate information about internal states such as attention, motivation, and reward prediction for coordinating plasticity (O’Donnell et al., 2012; Brzosko et al., 2017; Speranza et al., 2021). Whether learning-related transcriptional responses are sufficiently diverse to reflect such information remains poorly understood.
Previous studies using microarray-based transcriptomics have uncovered broad gene expression changes in M1 during motor learning (Cheung et al., 2013; Hertler et al., 2016). However, technical limitations prevented such efforts from differentiating between cell types. This is critical, as different types respond (or not) to learning-inducing experiences in various ways (Hrvatin et al., 2018). Additionally, plasticity phenotypes are divergent even within neuron types: both excitatory and inhibitory neurons display divergent forms of structural and functional plasticity over learning according to their subtype (Chen et al., 2015; Peters et al., 2017; Tjia et al., 2017; Ren et al., 2022). Therefore, different neuron types and those within the same type undergo heterogeneous activity changes and plasticity phenotypes, which is likely reflected in their transcriptional responses.
Among the various plasticity-related phenotypes, those related to forming engrams are of particular interest. Engram cells are repeatedly activated during continued learning and are thought to support memory recall (Semon, 1923; Tonegawa et al., 2015). Cells with such characteristics have been observed in many brain areas, including motor cortices. In the secondary motor cortex, for example, the initial expression level of the IEG Arc predicts a neuron's probability for reactivation on subsequent motor training days (Cao et al., 2015), implicating a transcriptional response reflective of engram formation. A recent study of novel environment exposure discovered a more comprehensive set of genes associated with reactivation in hippocampal dentate gyrus neurons (Jaeger et al., 2018), suggesting that engram-like reactivation involves a coordinated transcriptional program of many genes. In these studies, only a subset of cells exhibited the characteristics of engram cells, indicating that learning-induced transcriptomic changes are heterogeneous across neurons.
The recent advent of single-cell and single-nucleus RNA sequencing (scRNA-seq or snRNA-seq) has permitted the analysis of transcriptional responses in a way that intrinsically captures the heterogeneity across neurons. Using this approach, recent studies have begun to reveal transcriptional programs that differ both across (Hrvatin et al., 2018; Chen et al., 2020) and within (Chen et al., 2022) cell types during novel experiences and learning. However, most such studies have examined transcriptional responses to single-trial learning; thus, it remains to identify the various plasticity states and their underlying transcriptomic profiles that evolve during repetition-based motor learning in M1.
Therefore, our collaborative project MEMOry from NETwork (MEMONET) set out to investigate how the transcriptomes of individual neurons are altered in mouse M1 during repeated motor learning using snRNA-seq.
Materials and Methods
Animals
All animal procedures were performed in accordance with guidelines and protocols approved by the University of California San Diego (UCSD) Institutional Animal Care and Use Committee and the National Institutes of Health (NIH). Mice (C57Bl/6) were grouped and housed in disposable cages with standard bedding in a temperature-controlled and humidity-controlled room (∼21°C and 42% humidity) with a reversed light cycle (10.00–22.00: dark). All experiments were performed during the dark cycle.
Surgery
The headbar implant surgery was performed on six adult mice (9 weeks old, male, both train and control groups) under general anesthesia (1.5%–2% isoflurane). After the skull was exposed by an incision of the scalp, the periosteum was removed, and a custom headbar was glued to the skull near the lambdoid sutures. The exposed skull area was covered with dental cement, followed by subcutaneous injections of buprenorphine (0.1 mg/kg) and Baytril (10 mg/kg).
Water restriction
In 9–10 d after the surgery, water restriction started and lasted for 13 d. All mice received a controlled amount of water (1 ml/d) on a Petri dish in their home cages on the first 10 d. On days 11–13, the control group mice (n = 3) continued to receive 1 ml/d in their home cage, whereas the train group mice (n = 3) received 1 ml/d, in part as a reward during a lever-press task (0.3–0.8 ml) and the rest in their home cage immediately after completing the task.
Lever-press task
Mice in the train group performed a lever-press task as previously described. Briefly, a lever was positioned for head-fixed mice to press with their left paw (Fig. 1a). A trial began with a 6 kHz tone marking a cue period (10 s) during which a successful lever press was rewarded with water (∼15 μl/trial), followed by an intertrial interval (8–12 s). Success was defined as lever press crossing two thresholds (∼1.5 mm and ∼2.6 mm below the resting position) within 200 ms. With these criteria, a distinct and brisk downward press was considered a successful lever press. An unsuccessful lever press during the cue period triggered a loud white noise and an intertrial interval. Lever press during intertrial intervals was neither rewarded nor punished. Each daily session ended after mice completed 100 trials.
Experimental design and cell type classification. a, After 10 d of water restriction, mice (N = 3 for each group) are either subjected to 3 d of motor skill training or 3 d of home cage control. During motor training, mice must press a lever past the threshold after the presentation of an auditory cue, which results in a water reward (i.e., success trial). Failure to press the lever past the threshold results in no reward. After training on day 3, the primary motor cortex (M1) is extracted in both train and control groups for single-nucleus RNA sequencing (snRNA-seq). b, Lever-press trial success rate. Left: the average fraction of trials where the train mice gained reward from correctly pressing the lever past the threshold (mean ± SEM). By session 3 (1 session/day), the fraction rewarded is significantly increased (r = 0.69, p = 0.039, Pearson's correlation). Middle: Reward fractions are consistent with previously published data. The red lines indicate the performance of each mouse in the current study. The gray lines indicate the performance of mice from a previously published study (Hedrick et al., 2022). Right: change in reward rate (session X–session 1). The red lines indicate the change in reward rate per mouse in this study. The black line is the mean change in reward rate in the previously published data; the gray shaded area represents SEM. Change in reward rate over the three training sessions is comparable to that seen in Hedrick et al. (2022). c, Good consensus of the AIBS (left) and MEMONET datasets (right) in common UMAP space. The AIBS dataset is visualized in UMAP space computed by Seurat. During single-cell reference mapping, the MEMONET dataset was projected onto the AIBS UMAP space. Cells are colored based on cell type. Predicted cell types and prediction scores from mapping are shown in Extended Figure 1-1. d, Comparable neuronal cell type fractions are observed between the AIBS and MEMONET datasets.
Figure 1-1
Single-cell reference mapping results. https://figshare.com/s/cb98b2d50ad98ba6d158. Download Figure 1-1, XLSX file.
Tissue extraction and cryopreservation
All sacrifices were performed within 30 min from the completion of water consumption on water restriction day 13 (i.e., training day 3 for the train group). Mice were anesthetized via combined intraperitoneal injection of 150 mg/kg ketamine and 15 mg/kg xylazine. Once unconscious, mice were transcardially perfused with ice-cold, carbogen-bubbled (95% O2, 5% CO2), nuclease-free, 0.22 µm sterile-filtered artificial cerebrospinal fluid (ACSF) with a composition of 93 mM N-methyl-D-glucamine, 2.5 mM KCl, 1.2 mM NaH2PO4, 30 mM NaHCO3, 20 mM HEPES, 25 mM glucose, 5 mM sodium ascorbate, 2 mM thiourea, 3 mM sodium pyruvate, 13.2 mM trehalose, 12 mM N-acetyl-cysteine, 0.5 mM CaCl2, 10 mM MgSO4, and 93 mM HCl, at pH 7.3–7.4 (Ting et al., 2014; Tasic et al., 2018). Following transcardial perfusion, brains were immediately extracted and submerged into ice-cold carbogen-bubbled ACSF, with <5 min between the beginning of perfusion and final submersion after extraction.
All brains were serially sectioned in ice-cold, carbogen-bubbled ACSF on a VT1000S vibratome (Leica) with polytetrafluoroethylene-coated razor blades (Ted Pella) at 0.15 mm/s and 100 Hz, dividing the whole cerebrum into 400 µm coronal slices. Each slice between 0 and 1.2 mm anterior to bregma (2–4 slices per animal) was dissected under a stereomicroscope using a microsurgical stab knife to cut out the caudal forelimb area in the motor cortex of the right hemisphere (∼1–2 mm from the midline). All extracted tissue was recovered in ice-cold, nuclease-free, 0.22 µm sterile-filtered cryoprotective nuclear storage buffer, composed of 0.32 M sucrose, 5 mM CaCl2, 3 mM magnesium acetate, 10 mM Trizma hydrochloride buffer (pH 8.0), 1 mM dithiothreitol, 0.02 U/µl SUPERase•In RNAse Inhibitor (Invitrogen), and 1× cOmplete Protease Inhibitor Cocktail with EDTA (Roche). Tissue was then flash-frozen using a metal CoolRack M90 (BioCision) prechilled to −80°C and stored at −80°C until nuclear isolation.
Single-nucleus RNA sequencing
Single-nucleus RNA-seq was provided by the Center for Epigenomics–UCSD using the Droplet-based Chromium Single Cell 3′ solution (10x Genomics, v3.1 chemistry). Nucleus preparation was adapted from Lacar et al. (2016). The flash-frozen mouse motor cortex tissues were thawed on ice for ∼30–50 min before centrifugation for 3 min at 100 rcf in a 4°C, swinging bucket centrifuge with setting run 3/3 (5920, Eppendorf). The supernatant was removed, and the tissue was resuspended in 1 ml of douncing buffer consisting of 0.25 M sucrose (S1888, Sigma), 25 mM KCl (AM9610G, Invitrogen), 5 mM MgCl2 (194698, MP Biomedicals), 10 mM Tris-HCl pH 7.5 (15567027, Thermo Fisher Scientific), 1 mM DTT (D9779, Sigma), 1× protease inhibitor (05056489001, Roche), 0.1% Triton X-100 (T8787-100ML, Sigma), and 0.2 U/μl RNasin RNase inhibitor (PAN21110, Promega) in molecular biology grade water (46000-CM, Corning). The tissue was then transferred to a dounce homogenizer on ice and dounced 25 times with a loose plunger and 25 times with a tight plunger. The suspension was then passed through a Celltrix 30 μM filter (04-004-2326, Sysmex) and washed with 300 μl of douncing buffer before being centrifuged for 10 min at 1,000 rcf. The pellet was then washed with an additional 1 ml of douncing buffer without Triton X-100 and centrifuged again for 10 min at 1,000 rcf. The supernatant was then discarded, and the pellet was resuspended in 500 μl of sort buffer consisting of 1 mM EDTA (15575020, Invitrogen), 0.2 U/μl RNasin, and 1% BSA in PBS before staining with 3 μM DRAQ7 (#7406S, Cell Signaling Technology). Nuclei were then incubated on ice for 10 min, and ∼20–30,000 nuclei were sorted using a 100 μm chip in an SH800 sorter (Sony) into 50 μl of collection buffer consisting of 1 U/μl RNasin and 5% BSA in PBS. Samples were then centrifuged for 15 min at 1,000 rcf; the supernatant was removed, leaving behind ∼20 μl, and an additional 25 μl of reaction buffer consisting of 0.2 U/μl RNasin and 1% BSA in PBS was added for a total volume of ∼45–50 μl.
The prepared nuclei were visually inspected and manually counted using a hemocytometer before loading 15,000 onto a chromium controller for 10x GEM generation in the Single Cell 3′ Kit v3.1 (1000268, 10x Genomics). Libraries were generated using the Chromium Single Cell 3′ Library Construction Kit v3.1 (1000190, 10x Genomics) according to manufacturer specifications. cDNA was amplified for 12 PCR cycles. SPRISelect reagent (Beckman Coulter) was used for size selection and clean-up steps. Final library concentration was assessed by Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific), and fragment size was checked using Tapestation High Sensitivity D1000 (Agilent) to ensure that fragment sizes were distributed normally ∼500 bp. Libraries (one for each animal) were barcoded independently so that each cell is linked to the mouse it was derived from, before being pooled on the same Illumina lane for sequencing. Libraries were first shallow sequenced using a NextSeq 500 (Illumina) using the following read lengths: Read 1, 28 bp; Index 1, 10 bp; Index 2, 10 bp; Read 2, 90–91 bp. Libraries were then deep sequenced using a NovaSeq 6000 (Illumina) using the following read lengths: Read 1, 101 bp; Read 2, 101 bp.
Preprocessing of snRNA-seq data
Raw reads were processed using the 10X Genomics Cell Ranger v4.0.0 pipeline. Research associates at the Center for Epigenomics at UCSD used the mkfastq command to demultiplex the binary base call files and create fastq files; the count command was used, with parameter “–expect-cells” set to 7,000, to perform alignment, barcode, and unique molecular identifier (UMI) counting to generate feature-barcode matrices. The aggr command was used to aggregate the data for all mice. The parameter “–none” was used to turn off depth normalization, due to the requirement of the DESC clustering package for unnormalized counts as input.
To determine suitable filtering criteria, violin plots in Seurat v4 (Hao et al., 2021) were used to visualize the total number of unique genes, UMI counts, and mitochondrial content. The original dataset comprised 31,053 genes and 51,926 nuclei. Low-quality nuclei were removed based on the following criteria, with the remaining nuclei indicated in parentheses: total unique detected genes <200 (51,926 nuclei), total number of UMI counts <800 (51,329 nuclei), UMI counts >30,000 (50,688 nuclei), and mitochondrial content >1% (50,576 nuclei). The UMI counts >30,000 thresholds is a method suggested by Seurat to eliminate doublets (Satija Lab, 2023), which are pairs of cells that receive the same barcode in error. Genes that occur in <3 cells were also removed (24,901 genes remain).
AIBS-based cell type identification
Seurat v4 in R offers functions for single-cell reference mapping (Hao et al., 2021). It compares the gene expression profiles of cells in both the reference (Allen Institute for Brain Science, or AIBS, dataset; Yao et al., 2021) and query (our MEMONET dataset) and uses anchors to project the query cells into the reduced principal component analysis space of the reference dataset. Anchors are pairs of cells that represent the same cell state across both datasets (Stuart et al., 2019). They are found by jointly reducing the dimensions of both datasets and then searching for pairs of cells that are mutually closest to each other. These anchors help align the rest of the cells from the query to the reference. Every query cell receives a prediction score for every reference label, and the label with the largest score (prediction.score.max) becomes the predicted label for a given cell. Query cells are also projected onto the reference uniform manifold approximation and projection (UMAP) space computed in Seurat.
When mapping AIBS L2/3 cells onto MEMONET L2/3 cells, UMAP coordinates computed from DESC clustering were used in place of the Seurat-generated UMAP for visualization of the projection (Fig. 6a).
To assess the accuracy of Seurat's single-cell reference mapping on the AIBS dataset, each cell type (“cluster_label” from the AIBS dataset metadata file) was downsampled to 25% and used as a query against the remainder of the dataset. This was done 100 times, and the percentage of false classifications was calculated (Fig. 2).
Single-cell reference mapping accuracy. The accuracy of single-cell reference mapping was assessed by mapping test cells (25% of each cell type in the AIBS dataset) onto a reference formed with the remaining 75% of the AIBS dataset. The x-axis represents the original cell type of test cells (i.e., their original cell type label in the AIBS data set), and the y-axis represents their inferred cell type from the single-cell mapping to the reference. Off-target predictions can be seen as deviations from the diagonal. Counts were normalized per column to add up to a frequency of 1.
Logistic discriminant analysis
Logistic regression of the weighted linear sum of gene expression values was performed so as to minimize classification error of a given cell as originating from the “train” or “control” populations. Fitting was performed using the MATLAB function fitclinear with a logistic “learner” setting. Input gene expression values were z-scored in order to standardize the distributions. Regularization of logistic models was performed using the lasso method, with the regularization parameter λ selected by minimizing classification loss across a logarithmic range of λ inputs. The generalizability of the model was examined by performing a 10-fold cross-validation of the prediction accuracy of each k-fold logistic model. To this end, cells were indexed evenly into k = 10 groups using the MATLAB function “crossvalind,” with each of the k-fold groups being iteratively held out as test data for the kth logistic model trained on the kept 90% of the data. Regularization was performed independently for each k-fold model. Prediction accuracy for each k-fold model was calculated as the fraction of correct classifications (i.e., the sum of true positives for both labels divided by the total number of cells). To rank the contribution of each gene to the classification, the median value of each gene's calculated coefficient across the 10 models was sorted in descending order of their absolute values. The number of “target” genes was selected by optimizing the prediction accuracy of the cross-validated (CV) model against the number of top-weighted genes. “Top 3,000-only” and “top 3,000-excluded” models were constructed by using either the optimized number of target genes (“top 3,000-only”) or the remaining bottom-weighted genes (“top 3,000-excluded”), respectively.
“Shuffled” variants of the logistic model were generated by randomizing the labels of train versus control prior to model fitting and then repeating this process 100 times. All iterations converged on the expected chance value without any deviation, corresponding to the proportion of “control” cells.
Clustering
Cells identified as L2/3 were clustered using a deep learning Python package, DESC v2.4, to identify the subtypes (Li et al., 2020). DESC uses a deep neural network to reduce the multidimensional data into a low-dimensional space and is effective at removing batch effects, which arise when samples from different individuals or samples sequenced on different sequencing lanes have RNA counts that are not directly comparable due to technical noise.
Cell-level normalization is calculated as the UMI count for each gene in each cell divided by the total number of UMIs in the cell, multiplied by 10,000, and then transformed to a natural log scale. The dataset subjected to DESC consisted of the top 3,000 genes identified by the logistic classifier, termed experience-discriminating genes (EDGs). The expression level of each gene was standardized by subtracting the mean across all cells and dividing by the standard deviation across all cells for the given gene.
After normalization, a deep neural network is trained with the data, in which the network iteratively clusters the cells to gradually remove batch effects and refine clusters. The values for the “encoder” and “tol” parameters were chosen based on the information provided in Supplementary Table 2 by Li et al. (2020). Based on the clustering dataset size of 6,433 cells, encoder values were set to 6,432 and tol = 0.001. The n_neighbors parameter was set to 25 and louvain_resolution was set to 0.65. All other parameters were kept at default values.
Cell proportion difference between the train and control groups per cluster
To assess significance, the chance proportions of train and control cells assigned to each cluster in the EDG space were calculated using a permutation test without replacement, by shuffling the train/control labels 1,000 times. For each shuffle, the proportion difference between train and control was calculated for each subtype/cluster and used as the distribution for the null hypothesis (i.e., there is no difference between train and control groups).
p-values were calculated using the following formula detailed by Phipson and Smyth (2010):
Differential expression analysis
Differential gene expression (DE) analysis was run using DESeq2 v1.26.0 in R (Love et al., 2014). DESeq2 default parameters were optimized for bulk sequencing data, requiring some changes to be made based on recommendations tested by Ahlmann-Eltze and Huber (2021) and Zhu et al. (2019). The following DESeq arguments were set: useT = TRUE, minmu = 1 × 10−6, and minReplicatesForReplace = Inf. Set sizeFactors using scran::calculateSumFactors (this is the equivalent of computeSumFactors but works for non-sce objects) instead of the default function estimateSizeFactors.
Statistical significance is determined using the Wald test, which compares log-fold changes between two groups. A significance threshold of <0.05 was used for the adjusted p-values, which were based on the Benjamini–Hochberg method for multiple test correction. Log-fold change values were shrunken using the “apeglm” method.
Enrichment analysis (statistical comparison of cluster compositions)
Gene ontology (GO) enrichment analysis was performed using clusterProfiler v4.2.2 in R (Yu et al., 2012; Wu et al., 2021), and visualization was performed using functions in enrichplot v1.14.2 in R (Yu, 2022).
The function enrichGO() was used for running the overrepresentation test on differentially expressed gene (DEG) lists of individual clusters (or the EDG list). The following parameters were set: ontology = “BP” (biological process), pAdjustMethod = “fdr,” pvalueCutoff = 0.05, and qvalueCutoff = 0.05. The universe was set to all genes identified in the L2/3 cells after DESeq2 QC filtering. The results were visualized using the enrichplot function treeplot().
Reactivation score
The DEG list in Supplementary Table 3, tab “DG_React_NotReact” by Jaeger et al. (2018) was used as the reactivation-related gene list. Four numbers were calculated for each cluster: the number of upregulated genes in the cluster that are also upregulated in the reactivation-related gene list (up/up); the number of downregulated genes in the cluster that are also downregulated in the gene list (down/down); the number of upregulated genes that are downregulated (up/down); and the number of downregulated genes that are upregulated (down/up). The up- or downregulated genes in each cluster were identified from DE analysis of the given cluster in reference to all others. Then, the reactivation score was calculated using the ratio between the number of genes regulated in the same direction (proreactivation) as the reference and the number in the opposite direction (antireactivation) as follows:
Single-cell trajectory analysis
Single-cell trajectory analysis was performed using the Python package scanpy (Wolf et al., 2018), applied to DESC latent embedding of the cells from clusters C0–C3 (omitting C4 and C5). Specifically, we employed the method called diffusion pseudotime (DPT; Haghverdi et al., 2016) within scanpy to evaluate the cell pseudotimes. Here, the cell neighborhood graph was first constructed using UMAP (McInnes et al., 2020). The root cell—the starting cell with a pseudotime of 0—was identified from a force-directed graph representation of the DESC latent embedding, specifically using the Reingold–Tilford tree layout (Reingold and Tilford, 1981) and by choosing the cell from cluster C0 that was located at the root of this graph. Finally, to visualize the cell trajectory as informed by the pseudotime, we used a force-directed graph representation using the Large Graph Layout (LGL; Adai et al., 2004).
Experimental design and statistical analysis
Details of experimental procedures are described above (see “Water restriction” and “Lever-press task”). Statistical design for the logistic classifier can be found in Materials and Methods “Logistic discriminant analysis.” Statistical tests for cluster proportion are described in Materials and Methods “Cell proportion difference between the train and control groups per cluster.”
Data/code availability
The data and analysis code that support the findings of this study are available on request. RNA-seq data is available at GEO under accession number GSE224686.
Results
We performed snRNA-seq of M1 in two groups of water-restricted mice: train and control (Fig. 1a; N = 3 in each group). The train group was exposed to a well-established motor skill learning paradigm (a lever-press task) for 3 d (∼30 min/d), while the control group continued to receive a restricted daily amount of water in their home cage (Methods). In the lever-press task, mice needed to press the lever as a distinct, brisk movement within 10 s from the auditory cue to receive a water reward (Fig. 1a and Materials and Methods). The success rate for the trained group slightly but significantly increased over 3 d (r = 0.69, p = 0.039, Pearson's correlation), consistent with the trend observed in mice that performed the same task for 14 d in previous studies (Fig. 1b). While behavioral improvement becomes more apparent with further training, learning-related synaptogenesis peaks during the first 3–5 d in the lever-press task (Peters et al., 2014), suggesting that gene expressions driving learning-related plasticity may be most discernible before a rapid behavioral change. In support of this idea, previous bulk RNA assay in M1 found that the number of differentially expressed genes by motor training peaks just before the rapid phase of behavioral improvement began (Cheung et al., 2013). Thus, we extracted the M1 forelimb region on training day 3 for snRNA-seq (Materials and Methods).
While many cell types in M1 display various forms of plasticity that are likely critical to motor learning (Fu et al., 2012; Chen et al., 2015; Peters et al., 2017; Ren et al., 2022), we focus on L2/3 glutamatergic neurons as they are considered a critical locus where learning-related changes take place (Chen et al., 2015), and the time course of their structural and functional plasticity has been well established during the same learning paradigm (Peters et al., 2014; Hedrick et al., 2022).
Identification of L2/3 glutamatergic neurons based on single-cell reference mapping
Prior to analysis, quality control steps were taken to remove low-quality nuclei based on total RNA counts, unique gene counts, and mitochondrial gene percentage. The number of nuclei was reduced from 51,926 to 50,576. We then assigned cell type on a per-cell basis by using single-cell reference mapping to the previously published AIBS dataset generated from M1 snRNA-seq (Yao et al., 2021; Materials and Methods). This allowed us to levy the detailed cell type annotations in the AIBS dataset. For example, there are three identified L2/3 subtypes: L2/3 IT_1, L2/3 IT_2, and L2/3 IT_3. A side-by-side comparison of the AIBS dataset (Fig. 1c, left) and our dataset (hereafter referred to as MEMONET data) projected onto the same UMAP space (Fig. 1c, right) shows good consensus of the datasets with regards to cell type distributions. Expectedly, glial clusters are enriched in the MEMONET dataset compared with AIBS due to the AIBS dataset selecting for neurons before sequencing (Yao et al., 2021). Neuronal cell type compositions are comparable between the datasets (Fig. 1d). Inhibitory neurons account for ∼19% of neurons based on the mapping, similar to ∼16% in the AIBS dataset (Yao et al., 2021) and ∼17% according to Bakken et al. (2021).
The median prediction scores for the mapped cell type (ranging from 0 to 1) in our MEMONET dataset are 0.98 ± 6.9 × 10−4 across all cells and 0.97 ± 2.8 × 10−3 across the 6,507 L2/3 glutamatergic neurons (referred to as L2/3 neurons hereafter; Extended Fig. 1-1 and Table 1). The accuracy of the mapping pipeline was checked by mapping random subsets that contain 25% of the AIBS dataset onto the remainder of the dataset 100 times. On average, 93.3% of cells were consistently classified (Fig. 2 and Extended Fig. 1-1). The median prediction scores are 0.98 ± 9.3 × 10−5 across all cell types and 1.0 ± 3.2 × 10−17 for L2/3 neurons (Table 1). The comparable cell type compositions and prediction scores between the two datasets validate the integrity of our mapping pipeline and the MEMONET dataset.
Single-cell reference mapping prediction score summary
To further ensure that we described the heterogeneity among true L2/3 neurons, we excluded those cells with low mapping reliability (defined as the sum of scores for the three L2/3 subtypes defined by the AIBS dataset: L2/3 IT_1, L2/3 IT_2, and L2/3 IT_3 less than or equal to 0.3) in our subsequent analyses. We confirmed that the thresholding removes cells with low average expression of L2/3 marker genes and high average expression of the glial marker gene Mertk (Table 2). This thresholding filters out 74 cells from the predicted L2/3 neurons, leaving a total of 6,433 cells for downstream analysis.
Cutoff mapping scores excluded from L2/3
Identification of EDGs in L2/3 neurons
We first compared the ARG expression in L2/3 neurons between the train and control group, which revealed that Arc, BDNF, Nptx2, Ntrk2, and Scg2 show significantly higher expression levels in the train group (Table 3). Notably, significant training-induced activation of Arc, an important IEG for activity-dependent plasticity, is detected only in L2/3 neurons. Nevertheless, the inevitable variability in gene expression, reflecting both biological and technical noise, poses a significant challenge when attempting to accurately distinguish the experimental conditions based solely on the expression of a single gene. To address this, we employed a population-code approach commonly used in neural data analysis, which deciphers information encoded collectively by a group of neurons rather than by each neuron in isolation. Population-code methods can leverage the correlation structure among different neurons to enhance decoding accuracy (Averbeck et al., 2006). Similarly, to account for both individual gene expression levels and the correlations between them, we trained a logistic classifier to discriminate “train” versus “control” neurons based on normalized expression values of all genes (Materials and Methods). The logistic classifier identified the weight of each gene such that the weighted sum of all gene expressions and thus its transformed probability to be “train” yielded good separation between train and control neurons (Fig. 3a). Consequently, the classifier showed a high CV group discrimination accuracy of ∼87% (Fig. 3b, left). The discrimination accuracy was consistently high across all mice [97 ± 1.7% (mean ± SD); 94–99% (range)]. To understand the predictive power of this model over chance, we randomized the labels of train and control neurons 100 times and trained a logistic classifier to discriminate these labels as before. In this context, prediction accuracy consistently converged to ∼60% for each shuffle (Fig. 3b, right). This percentage is identical to the proportion of control neurons in the population, as the classifier defaulted to always predicting “control” to maximize prediction accuracy. Therefore, the high prediction accuracy of the logistic classifier in our unshuffled data is not due to potentially spurious effects such as a random chance to have a certain number of genes that are differentially expressed between any two randomly divided groups of neurons. Instead, this result indicates the existence of genes whose expression levels reflect the true difference between the two experiences (motor training vs home cage control), and the logistic regression classifier could discover those genes that discriminate the experiences. Thus, logistic regression of gene expression values can accurately predict whether an L2/3 neuron is derived from a recently trained animal.
Logistic classifier details. a, Left: placement of sub-sampled cells on the logistic curve fit to the full model. The x-axis represents the weighted sum of all genes with the weights identified by the full model (i.e., evidence for category “train” vs “control”). The y-axis represents the logistic transformation of the evidence (i.e., the probability to be labeled “train”). The colored labels correspond to the true identities (i.e., “train” vs “control”) of individual neurons. Right: histogram of the probability to be labeled “train” for train (green) versus control (gray) cells. Control and train cells are well separated along the y-axis, where 0.5 corresponds to the decision boundary to be labeled as “train.” b, Left: confusion matrix of the full logistic model (i.e., with all genes included) showing the accordance between true cell labels (train vs control) and labels predicted from the logistic model. The colors represent the fraction of cells for a given true label (rows) occupying each quadrant. Right: the same as left but for randomly shuffled cell labels. c, Left: sorted absolute values of gene weights (i.e., regressor coefficients) for all genes from the full model plotted against the number of genes. “Top genes” were defined based on sorted weights. The red dashed line indicates the ultimate selection of the cutoff for inclusion in future analyses (defined on the right). Right: optimization of the number of genes for consideration based on logistic model prediction accuracy. Logistic regression was repeated using a varying number of top genes as inputs, with the resulting prediction accuracy for each iteration shown. Prediction accuracies reached a plateau when including ∼3,000 top genes. d, Confusion matrices for logistic models using only the top 3,000 genes (left) or excluding the top 3,000 genes (right). e, The probability to be labeled “train” in the “top 3,000 genes only” model per cell, sorted/colored by mouse, showing model performance across animals. Prediction probabilities were calculated using the mean coefficients across all cross-validations. The three shades of green represent the three train mice. The three shades of gray represent the three control mice. f, The prediction accuracy of the model with various subsets of genes included. The full model (gray) includes all genes; top 3,000-excluded model (black) has the top 3,000 genes removed; top 3,000-only model (red) includes only the top 3,000 genes. The dots represent values for each cross-validation. Prediction accuracy for the top 3,000-excluded model is chance, indicated by the cyan line (0.6 is the proportion of control cells in the population). The value of the cyan line was calculated by randomizing the labels of train and control cells 100 times and training the logistic classifier to discriminate the labels. CV discrimination accuracy for the top 3,000-only model is 96.92%. The results of DE analysis for the top 3,000 genes are shown in Extended Figure 3-1, and GO enrichment analysis results are shown in Extended Figure 3-2.
Figure 3-1
Differential expression analysis results. https://figshare.com/s/4b52bbb88a5a664cb9e4. Download Figure 3-1, XLSX file.
Figure 3-2
Gene Ontology enrichment results. https://figshare.com/s/44963f27d654243a2759. Download Figure 3-2, XLSX file.
IEG expression changes induced by training across neuron types
We next sought to identify the genes that contributed the most to discrimination between “train” and “control” neurons. Since both up- or downregulation of gene expression could contribute to successful discrimination, we ranked the median absolute values of the coefficients (weights) in the logistic regressor for each gene across all CV iterations. As shown in Figure 3c (left), the magnitude of gene coefficients decreased sharply in an approximately double-exponential fashion. While many thousands of genes showed nonzero coefficients, we sought to select the top-weighted genes that maximized prediction accuracy. We therefore repeated logistic regression while varying the number of top genes as inputs. Prediction accuracy began to plateau with the inclusion of 3,000 top genes (Fig. 3c, right), so we proceeded to use this value as a cutoff (Materials and Methods). To test the validity of this cutoff, we retrained logistic models using either the top 3,000 genes only (“top 3,000-only model”) or instead using all genes except the top 3,000 (“top 3,000-excluded model”; Fig. 3d). The top 3,000-excluded model produced CV prediction accuracies at chance level, while the top 3,000-only model showed a high prediction accuracy (Fig. 3d). Consistently, the probability to be “train” in the top 3,000-only model was well separated between train and control neurons across all mice (Fig. 3e). Furthermore, we found that the top 3,000-only model—despite having fewer regressors—showed increased prediction accuracy compared with the full model (Fig. 3f). Thus, while trace discriminating information is likely to present in the remaining gene list, the top 3,000 genes returned by the full logistic regression model contain the majority of discriminating information between control and train neurons. Of these 3,000 genes, 1,147 showed significantly different expression levels in L2/3 neurons between the train and control group according to the single-gene basis DE analysis (Extended Fig. 3-1). Therefore, the classifier analysis could identify gene sets that are coregulated by experience, such that their correlated response patterns can discriminate experience even though some individual genes, by themselves, may not reach a significant difference between different experiences. For example, Kalrn, Acp4, Syp, Star, and Egr1 are identified by the classifier analysis but not DE analysis, and are related to regulation of neuronal synaptic plasticity (Extended Fig. 3-2).
To infer the biological processes associated with these 3,000 genes, we performed a GO enrichment analysis (Yu et al., 2012; Wu et al., 2021). This analysis yielded enrichment for terms such as learning (padj < 4.7 × 10−2), regulation of synaptic plasticity (padj < 3.5 × 10−3), and synapse organization (padj < 3.2 × 10−7) (Fig. 4 and Extended Fig. 3-2) and highlighted coregulation between commonly studied plasticity-related genes such as kinases (e.g., Camk2b, Ntrk2) and glutamate receptors (e.g., Gria1, Grin2a; Fig. 4, red text). These results indicate that the subset of genes discriminating between the train and control neurons in our L2/3 population are related primarily to cellular plasticity and learning. We therefore hereafter refer to these 3,000 genes as EDGs (Extended Fig. 4-1). We performed GO analysis separately for positively versus negatively weighted EDGs and found many GO pathways that are significantly associated with the positively weighted EDGs, but none with the negatively weighted EDGs (Extended Fig. 3-2). The pathways associated with the positively weighted EDGs include many plasticity-related terms, consistent with the idea that training results in the upregulation of genes that support learning-related plasticity.
GO enrichment analysis of the 3,000 EDGs. Eight pathways of interest are shown with associated genes. Well-studied genes are highlighted in red. The expression level of train and control cells is given as the mean z-score per gene, where green and magenta indicate positive and negative scores, respectively. Each pathway is broken down by the positively weighted and negatively weighted genes produced by the logistic model, and genes are sorted based on the magnitude of their regression weights. Full GO results are shown in Extended Figure 3-2. The 3,000 EDGs are shown in Extended Figure 4-1.
Figure 4-1
EDG list. https://figshare.com/s/b851d113dc339c6ad57c. Download Figure 4-1, XLSX file.
Distinct modes of coregulation among EDGs in L2/3 neurons
During the lever-press task, activity patterns differ greatly among L2/3 neurons (Peters et al., 2014); thus, it is unlikely that training affected the EDGs uniformly across all L2/3 neurons. Instead, within each neuron, a subset of EDGs might be coregulated in a manner dependent on their activity, such that neurons with similar activity patterns show a similar mode of EDG coregulation. To test this prediction, we applied unsupervised clustering to the 6,433 L2/3 neurons based only on the expression of the EDGs defined using the above logistic classification (Materials and Methods). This approach generated six clusters of neurons (C0–C5) representing distinct coregulation profiles of the EDGs (Fig. 5a), each with varying proportions of train and control neurons (Fig. 5b). By comparing the observed versus expected proportions of train and control neurons in each cluster, we find that C3 has significantly more train neurons, while C0 and C1 have significantly more control neurons (Fig. 5b). The increased proportion of C3 in the train group is consistent across mice, as indicated by the lowest C3 proportion in the train group being still higher than the highest in the control group. Despite the proportional differences between the two groups, nonnegligible representations of both train and control neurons were present in all clusters. These results suggest that motor learning does not generate a unique transcriptional state per se but rather shifts the relative occupancy of states (i.e., the transcriptional landscape) along a conserved transcriptional axis related to learning. This is consistent with previous observations that plasticity events, such as dendritic spine dynamics, are expressed in animals that are not undergoing specific learning tasks, albeit at lower levels compared with task-learning animals (Hedrick et al., 2022).
Distinct clusters of L2/3 glutamatergic neurons. a, L2/3 clusters. Left: Clustering was performed using the Python package DESC, based on the normalized expression of the EDGs. Six distinct clusters (C0–C5), each color-coded, were visualized in UMAP space. Right: The same as the left, but the train and control cells are green and gray, respectively. b, Proportion of train and control cells in each cluster. C0 and C1 are significantly enriched with control cells (p = 0.012 and p = 0.003, respectively). C3 is significantly enriched with train cells (p = 0.001). *p < 0.05, **p < 0.01, ***p < 0.001 (permutation test). The dots indicate the proportion of each mouse, where train mice are shades of green and control mice are shades of gray. c, Average IEG expression across clusters. Z-scoring was applied on a per-gene basis to counts normalized by DESeq2 before averaging across cells in each cluster. Significance is based on DE analysis of one cluster versus the others, padj < 0.05. The dull results of the DE analysis are shown in Extended Figure 3-1.
What is the biological basis of the distinct transcriptional states represented by these clusters? One attractive possibility is that such delineations reflect sequential stages of plasticity through which neurons transition as animals learn. Such a model meshes well with the observation that L2/3 neurons advance from a stage of enhanced activity and synaptogenesis to synapse elimination and sparsened activity as motor learning progresses (Peters et al., 2014). We therefore hypothesized that the transcriptional states captured by our clusters would mirror the ongoing dynamics of activity and plasticity in these neurons. To investigate this possibility, we first compared the expression levels of IEGs that are known to be transcribed by neuronal activity. We found that IEGs were indeed differentially expressed across clusters (Fig. 5c and Extended Fig. 3-1), supporting the possibility that different clusters might represent neurons with different levels of activity and thus engaging different activity-dependent genes underlying plasticity. Specifically, we find that C3, which was enriched by training, shows the highest expression in 12 out of the 16 profiled IEGs, indicating that the neurons of this cluster had likely undergone strong recent activity. Indeed, among the upregulated IEGs was Npas4, the induction of which is highly specific to neuronal activity as opposed to other activities (e.g., trophic factor stimulation; Lin et al., 2008) and primarily reflects very recent activity. Thus, C3 likely corresponds to neurons that were active on the day of tissue collection. Additionally, C2 and C5 also show some IEG upregulation, although to a lesser extent than C3. By contrast, C0, C1, and C4 show lower expression levels in IEGs, suggesting that these clusters correspond to a relatively inactive population. Thus, different groups of L2/3 neurons identified from the EDG space also differ in the levels of recent activity.
Having observed the diversity of IEG expression, we next sought to define a putative “baseline” cluster that best represents the neurons with relatively low activity levels and plasticity. Given the relatively sparse encoding of movements by L2/3 neurons in this lever-press task (typically 20–30% of neurons are significantly correlated with movements; Peters et al., 2014), we suspected that any such baseline population would comprise a large fraction of all neurons and would show the lowest IEG expression. Consistent with these criteria, C0 represented among the largest clusters (∼23% of the total population; Table 4), and overall showed the lowest expression levels of IEGs (Fig. 5c). Thus, we hypothesized that C0 may represent the baseline state with low neural activity and low activity-dependent gene expression. If so, we predict that AIBS L2/3 neurons sequenced after no explicit novel experience would occupy C0 more frequently. To further test this idea, single-cell reference mapping was used to map AIBS L2/3 neurons onto our MEMONET L2/3 clusters, annotating the AIBS cells by MEMONET clusters (Fig. 6a and Extended Fig. 1-1). The majority of AIBS L2/3 neurons are assigned to C0, C1, and C2 (Fig. 6b), with almost half of AIBS neurons (49.6%) mapped to C0, further supporting C0 as a putative baseline state.
Single-cell reference mapping of AIBS L2/3 neurons to MEMONET L2/3 clusters. a, AIBS L2/3 neurons projected onto MEMONET UMAP space, colored based on the mapped cluster. Predicted cluster assignment and prediction scores from mapping are shown in Extended Figure 1-1. b, Proportion of AIBS L2/3 neurons mapped to MEMONET L2/3 subclusters. The majority of AIBS neurons map to C0–C2.
Distribution of L2/3 neurons across subclusters C0–C5
We, therefore, performed a DE analysis comparing each cluster to C0 as a reference for the baseline state and performed GO enrichment analysis for the upregulated DEGs of each cluster (Extended Figs. 3-1 and 3-2). In this context, the top enriched GO terms for C4 include characteristics of many non-L2/3 cell types, for example, osteoblast differentiation and kidney epithelium development (Fig. 7), leading us to speculate that C4 may contain remaining contaminants that were not filtered out with our L2/3 prediction score cutoff described earlier. To investigate this, we compared the average expression per cluster of a handful of L2/3 marker genes (Cux2, Otof, Rtn4rl1, Slc30a3, Cacna2d3) and the glial marker Mertk. C4 has the lowest expression of L2/3 marker genes and the highest expression of Mertk, suggesting that C4 may represent contaminating cell types (Table 5). Therefore, C4, which constitutes <5% of MEMONET L2/3 neurons (294 cells) and only nine AIBS L2/3 neurons, was removed from subsequent analyses.
GO enrichment analysis of C4. The top 30 most significant terms (based on adjusted p-value) are shown. The upregulated differentially expressed genes from DE analysis of C4 versus C0 were used as input for GO analysis. The majority of the top terms include characteristics of other cell types. The full GO results are shown in Extended Figure 3-2.
Glia contamination in the L2/3 subcluster
Notably, C3 showed significant upregulation for 15 out of 16 IEGs (Fig. 8a). The other clusters also present stronger IEG signals when compared with C0. C1 showed significant upregulation of multiple IEGs, sharing Arc, Nefm, Nr4a2, and Nr4a3, while C2 and C5 have 7–10 additional IEGs significantly upregulated compared with C0. Specifically, C2, C3, and C5 show upregulated Npas4 and Fos transcription, indicating strong recent activity. Taken together, these analyses support the notion that C0 corresponds to a relative floor of IEG expression and thus recent activity, substantiating its consideration as a baseline state.
Plasticity-related phenotypes of distinct L2/3 clusters. a, IEG expression is shown as Log2 fold change, a metric used by DESeq2 to assess how much a gene's expression has changed between the two comparisons. A positive value indicates upregulation compared with C0 while a negative indicates downregulation. In the DE analysis of each cluster versus C0, most clusters show upregulation of a majority of IEGs. Significance is based on padj < 0.05. The full results of the DE analysis are shown in Extended Figure 3-1. b, GO enrichment analysis of selected plasticity-related terms across L2/3 clusters. Top: mean z-scored expression of all genes contributing to each term (mean ± SEM). Bottom: heat maps of z-scored expression of individual genes in each term. For visualization purposes, a maximum of 40 genes are shown in heat maps. Data for line plots include all genes for a given term. The full GO results are shown in Extended Figure 3-2.
To better understand the biological processes in which each cluster is involved, we examined the enriched terms from GO enrichment analysis for the upregulated DEGs of each cluster with respect to C0 (Extended Fig. 3-2). Learning- and plasticity-related terms dominated the most enriched results for C1–C3, implicating these clusters as specific transcriptional states relevant to motor learning. To identify potential differences in the gene expression profiles associated with such terms across clusters, we next compared the mean expression of associated genes for each term across clusters. C1–C3 showed varying levels of genes associated with different aspects of plasticity, including long-term potentiation (LTP), long-term depression (LTD), and actin filament-based process(es) implicating dendritic spine morphogenesis (Fig. 8b). More specifically, LTP- and LTD-related gene expression peaks at C2, indicating that C2 might represent the neurons that were strongly activated recently and undergoing early LTP and LTD (Fig. 8b). Figure 9 shows an expanded list of behavioral processes enriched in C2 and C3. C2 is notable for its strong expression of genes associated with changes in the intrinsic excitability of neurons (Fig. 9), consistent with reports that such changes can co-occur with synaptic plasticity (Xu et al., 2005). At a higher level, C2 shows the highest expression for genes related to locomotory behavior and walking behavior, implicating the transcriptional state represented by C2 in volitional limb movements related to task performance and/or enhanced exploratory locomotion observed in water-restricted mice (Tsunematsu et al., 2008). Most genes in C1 show intermediate levels of expressions between C0 and C2, suggesting that C1 phenotypes might be an intermediate stage between C0 and C2.
Expanded list of plasticity-related GO terms across clusters. The line plots and heat maps of mean z-scored expression levels for genes associated with plasticity-related GO terms, as in Figure 8b. Terms were loosely categorized based on their functional relevance (vertical labels at left), as follows: (1) “functional (synaptic)” refers to processes related to the efficacy of chemical synaptic transmission at excitatory (glutamatergic) synapses; (2) “structural (synaptic)” denotes processes related to the structural basis of dendritic spines or their filopodial precursors; (3) “molecular” broadly refers to molecular pathways, events, or signal transduction; (4) “intrinsic” refers to nonsynaptic processes that govern the intrinsic excitability of neurons relevant to action potential firing; (5) “behavioral” denotes higher-order phenomena relevant to the intact, behaving animal. As in Figure 8b, the line plots correspond to the mean z-scored expression for all genes from a given GO term, while the heat maps depict a maximum of 40 genes for visualization purposes. For the heat maps, genes were sorted in ascending order according to the clusters that showed the highest expression. For GO terms with >40 associated genes, the gene list was first sorted according to the regression weights from the “top 3,000 genes only” logistic model in order to highlight genes that best discriminated between “train” versus “control” neurons, and the top 40 of these sorted genes were used. The full GO results are shown in Extended Figure 3-2.
Like C2, C3 exhibits high expression of genes associated with numerous aspects of synaptic plasticity, including LTP/LTD and related processes, structural dynamics of dendritic spines, and terms related to learning and memory. Consistently, both clusters also display higher expression of genes related to molecular pathways relevant to plasticity, such as kinase activity and neurotrophin signaling (Fig. 9). In contrast, however, C3 shows particularly strong expression of genes associated with the regulation of actin-related process(es) (Fig. 8b), including the regulation of filopodium assembly (Fig. 9). As actin is the primary cytoskeletal constituent of dendritic spines and filopodia represent the structural predecessors of dendritic spines, these results might reflect the increased requirement for actin polymerization to support the ongoing spinogenesis expected at this stage of motor learning (Peters et al., 2014; Hedrick et al., 2022). Further delineating C3, we also found that this cluster presents the highest expression levels of genes associated with mRNA processing (Fig. 8b), posttranscriptional regulation of gene expression, and regulation of translation (Fig. 9). These processes collectively suggest that C3 represents a more advanced stage of plasticity than C2, such as late LTP, which requires de novo protein synthesis (Kelleher et al., 2004; Ho et al., 2011; Baltaci et al., 2019). Compatible with the idea that this cluster corresponds to later stages of plasticity, C3 also showed the highest expression levels of genes pertaining to proteasome-mediated ubiquitin-dependent protein catabolic process (Fig. 8b), a form of posttranslational control of protein levels that has been implicated in the maintenance of late stages of LTP (Dong et al., 2008) and new spine formation (Hamilton et al., 2012). Further, C3 also presents with high expression of genes related to histone modification (Fig. 9), raising the possibility that the plastic changes observed in this cluster extend to more enduring changes through epigenetic modification for the control of future gene expression. C3 thus represents a functional convergence of plasticity-related phenotypes with upregulated posttranscriptional and posttranslational mechanisms associated with the maintenance of learning-related structural and functional changes to neurons.
Similar to C3, C5 shows high-level IEG expression (Fig. 8a) but low-level expressions of plasticity-related genes (Fig. 8b). Instead, C5—equally represented in both train and control groups—has the highest expression of genes involved in the ATP biosynthetic process and cellular respiration. In fact, many of the top enriched GO terms for C5 are related to metabolic processes such as mitochondrion organization and oxidative phosphorylation (Fig. 10a). Brain energy metabolism supporting neural activity and circuit remodeling is largely regulated by brain states (DiNuzzo and Nedergaard, 2017). In the active awake state with high arousal, for example, intracellular ATP levels in cortical excitatory neurons and their firing rates are increased (Natsubori et al., 2020). The active awake state has also been associated with increased glycolysis in the brain (DiNuzzo and Nedergaard, 2017). We found upregulated expressions of glycolysis-related genes in C5 (Fig. 10b), suggesting that C5 may mediate brain state-related metabolic changes in M1. The control of arousal and neuronal metabolism highly relies on norepinephrine (NE), a neuromodulator broadcast mainly from locus ceruleus (LC; O’Donnell et al., 2012; DiNuzzo and Nedergaard, 2017). Linking NE, arousal, metabolism, and neural activity, we found the highest level of NE receptor gene expressions in C5 (Fig. 10b). These results suggest that experimental conditions common to the train and control groups—but absent in the AIBS experiments—might enhance the release of arousal-related NE, revealing a subpopulation of neurons that responds with increased metabolic activity and neural activity (viz., C5). Among the notable experimental differences in our experimental approach that might explain C5 are water restriction and circadian rhythm timing (as our experiments were conducted in the dark cycle, in contrast to the AIBS data, which was acquired in the light cycle; Hongkui Zeng, personal communication; Yao et al., 2021, 2023).
GO enrichment analysis of C5. a, Top 30 most significant terms (based on adjusted p-value) are shown for C5. The upregulated differentially expressed genes from DE analysis of C5 versus C0 were used as input for GO analysis. The full GO results are shown in Extended Figure 3-2. b, Upregulated processes in C5. Top: mean z-score of associated genes (mean ± SEM). Bottom: the z-scored expression of individual genes. The gene set for glycolysis was derived from the PANTHER database (Thomas et al., 2022).
Sequential relationship among transcriptionally distinct L2/3 clusters
As described earlier, some clusters were significantly enriched or de-enriched by training (Fig. 5b), which raises the possibility that individual neurons can dynamically switch between the states represented by our clusters and thus that the fraction of neurons in each cluster can change depending on recent experiences. In particular, C0–C3 is associated with baseline low neural activity, early LTP or LTD stages, and late-stage LTP, respectively, suggesting that they represent different stages of activity-dependent plasticity that gradually and sequentially evolve. To test this idea, we performed a single-cell trajectory analysis that derives the sequential relationship among loosely connected clusters, relying on the transcriptional similarity between L2/3 neurons in the EDG space (Wolf et al., 2018; Fig. 11a; Materials and Methods). We found a sequential order among the 4 clusters in terms of pseudotime as follows: C0→C1→C2→C3. These results suggest that C3 neurons might have been in a C2 state on previous days, while some of the C2 neurons would transition to C3 on the following days. In other words, L2/3 neurons may transition between these different clusters depending on the recent experience (i.e., plasticity of transcriptional responses).
Trajectory analysis and reactivation score of clusters. a, Single-cell trajectory analysis (DPT) reveals sequential transition between clusters. The cell trajectory informed by the pseudotime is visualized via a force-directed graph representation using the LGL. b, The number of significantly regulated genes in the direction indicating proreactivation versus antireactivation in each cluster. Genes that the reference (Jaeger et al., 2018) previously found to be significantly regulated in reactivated neurons were considered. Up/up: upregulated in both the MEMONET and reference dataset. Down/down: downregulated in both datasets. Up/down: upregulated in MEMONET but downregulated in reference. Down/up: downregulated in MEMONET but upregulated in reference. The significance of regulation in the MEMONET dataset was based on the DE analysis of each cluster versus the others. Numerical data are shown in Extended Figure 11-1. c, Reactivation score per cluster, calculated as (up/up + down/down)/(down/up + up/down) − 1. A score greater than 0 indicates that a larger number of genes are regulated in the proreactivation direction.
Figure 11-1
Reactivation score. https://figshare.com/s/dc4fb5de81625aaa96e5. Reactivation score was calculated as follows: (up/up + down/down) / (down/up + up/down) −1. A score greater than 0 indicates that a larger number of genes are regulated in the direction consistent with reactivation. MEMONET/REFERENCE columns indicate the number of differentially expressed genes (DEGs) in each cluster corresponding to pro-reactivation and anti-reactivation. Up/up: upregulated in the cluster, i.e., MEMONET, and upregulated in the reference. Down/down: downregulated in MEMONET and downregulated in the reference. Up/down: upregulated in MEMONET but downregulated in reference. Down/up: downregulated in MEMONET but upregulated in reference. DEGs for each cluster were chosen from the differential expression (DE) analysis of each cluster vs. the others. Reference gene list is taken from (Jaeger et al., 2018) (see Methods). Download Figure 11-1, XLSX file.
In the AIBS dataset, L2/3 cells were further differentiated into three subtypes IT_1–3 (Yao et al., 2021). To examine how these subtypes correspond to our functional clusters, we performed reference mapping of these subtypes onto our clusters C0–C5 (Fig. 6). We found that IT subtypes were mapped almost exclusively to C0–C2, consistent with our earlier interpretations that (1) C3 and C5 are uniquely emerged functional states in our study, and (2) C4 primarily comprises contaminated nonneuronal cells. Furthermore, IT1 mapped predominantly to C0, and IT3 predominantly mapped to C2, while IT2 mapped to both C1 and C2 (Table 6). This observation suggests that IT subtypes in the AIBS dataset may be differentiated, at least in part, by activity-dependent functional states underlying the spontaneous plasticity of L2/3 cells.
Distribution of AIBS L2/3 subtypes across C0–C5
Learning-induced enrichment of reactivation
Strong synaptic activity can induce the potentiation of dendritic spines via NMDAR-Ca2+-dependent, actin-based morphogenesis (Matsuzaki et al., 2004; Okamoto et al., 2004), facilitating their reactivation in the future (Chen et al., 2013). A previous study has shown that reversing the plasticity of previously potentiated spines can erase motor memories, supporting the notion that the potentiation of dendritic spines facilitates the generation of synaptic motor engrams (Hayashi-Takagi et al., 2015). Consistent with this notion, previously synaptically activated neurons are more likely to be reactivated during motor learning (Cao et al., 2015). We hypothesized that C3, showing a signature of late LTP [i.e., high actin filament-based process(es) and posttranscriptional/translational processes; Fig. 8b] and elevated IEG expression (Fig. 8a), might represent repeatedly active (i.e., reactivated) neurons over the course of the 3 d training. To test this hypothesis, we examined the expression profile of previously reported reactivation genes (Jaeger et al., 2018) and assessed the likelihood of reactivation for each cluster (i.e., “reactivation score”; Materials and Methods). The reactivation score was the largest in C3, while other clusters show no evidence of reactivation, agreeing with our hypothesis (Extended Fig. 11b–c). It has also been shown that the fraction of reactivated neurons is higher in motor skill-trained animals than in control (Cao et al., 2015). Thus, the significant enrichment of C3 in our MEMONET dataset compared with the AIBS dataset, coupled with the enrichment of C3 by motor training versus control as described earlier, is compatible with the reactivation hypothesis (Figs. 5b, 6b). These results collectively suggest multiple interrelated processes in C3 neurons: (1) the upregulated expression of genes associated with mRNA processing and posttranslational processes in C3 neurons, likely reflecting the protein synthesis and maintenance during late LTP; (2) the de novo synthesis of proteins involved in forming, enlarging, and maintaining potentiated dendritic spines; and (3) the reactivation of stably potentiated dendritic spines and consequently the soma of C3 neurons. Additionally, C3 shows upregulated expression of Baz1a, which a previous study of L2/3 neurons in S1 found to be a marker for neurons that can maintain plasticity in response to stimuli (Condylis et al., 2022). Consistent with the single-cell trajectory analysis above, we propose that L2/3 neurons progressively enter more advanced plasticity states from C1 to C2 to C3, depending on their levels of reactivation. Together, these findings suggest that C3 might contribute to the consolidation of motor learning across the multiday training, representing possible engram cells.
Discussion
L2/3 neurons have been shown to exhibit high plasticity during the early stage of motor skill learning, showing an increase in dendritic spines by day 3 of learning (Peters et al., 2014). However, the patterns of gene expression that allow for plasticity states to arise are not fully understood. Here, we used snRNA-seq of mouse M1 after 3 d of learning a lever-press task to investigate the transcriptomic profiles of L2/3 neurons during their initial increase in plasticity. Gene expression levels inherently vary across cells, even within the same cell type and under identical experimental conditions, due to biological and measurement noise. To overcome the inherent variability of a single gene and accurately classify cells from different experimental groups (i.e., train or control), we employed a logistic regression-based approach. This method identifies the optimal linear combination of gene expression levels that maximizes separation between these groups by utilizing both the DE of individual genes as well as correlations between them. The advantages of logistic classifiers over other commonly used methods have been reported (Ntranos et al., 2019). The logistic regression-based classifier revealed EDGs that together can reliably predict the group identity of each L2/3 neuron. These EDGs were coregulated in multiple distinct modes, forming separate clusters among L2/3 neurons. We found that these clusters represent different levels of neural activity, metabolic activity, and phenotypes consistent with different phases of plasticity across L2/3 neurons. Learning-inducing experiences produce transcriptional responses of EDGs, changing the fraction of neurons in these clusters.
The putative baseline state, C0, that showed the lowest IEG expression in reference to the other clusters was enriched in the control group, suggestive of quiescent neurons that are not stimulated or activated by recent experience. Consistent with this idea, approximately 50% of AIBS L2/3 neurons mapped to this cluster, indicating that C0 is the most dominant state in the absence of salient experience. In our dataset, fewer neurons occupied C0 (∼23%), suggesting that many neurons under the water restriction and motor learning during the dark cycle were active, transitioning to states enriched in our dataset such as C3 or C5. In addition, 23% is also lower than the estimate from the fraction of movement-related neurons observed in the same lever-pressing task (Peters et al., 2014). This discrepancy is accountable, at least in part, as movement-related neurons are a subset of active neurons. In other words, the fraction of active neurons would be larger than movement-related neurons. By using trajectory analysis with C0 as the starting point, we found the progression between clusters to be C0→C1→C2→C3. Trajectory analysis results imply that C1 might be a transitory state that bridges the baseline (C0) and activity (C2 and C3) states. A confirmation of this idea will require future longitudinal analyses of gene expression over the course of learning.
C2 appears to represent neurons activated during motor training on day 3, but not (consistently) activated on previous days. While C2 had some IEG upregulation in reference to the other clusters and even more upregulation in reference to C0, its reactivation score was near chance, suggesting these neurons have only been activated on 1 d of motor training. The upregulation of multiple IEGs, and in particular Npas4, in reference to baseline suggests that C2 was activated on day 3, rather than day 2, as such an IEG signature would not be expected when sequenced a full day later. That is because, unlike other IEGs that may be upregulated by such signals as trophic factors, Npas4 has a specific relationship with neuronal action potentials: it is only induced via calcium influx (Lin et al., 2008; Ramamoorthi et al., 2011). Additionally, it displays a fast half-life, with mRNA being degraded within hours (Lin et al., 2008). Thus, Npas4 is a dependable signal for neurons actively firing, which in this case is presumed to arise due to repetitive execution of unfamiliar movements. Genes involved in LTP and LTD processes expressed peak expression in C2, suggesting these neurons are undergoing early LTP and LTD. There exist L2/3 neurons in the AIBS dataset that map to C2 (∼26%). Given that C2 is also not significantly enriched in train or control groups, this state can perhaps arise in the absence of a learning task simply by the mouse experiencing and interacting with the environment.
On the other hand, C3 was more or less unique to the MEMONET dataset, representing only 1.68% of AIBS L2/3 neurons. This unique state seems to be attributed to reactivation during the learning process. One caveat to using the reactivation genes from Jaeger et al. (2018) is that that study profiled a different brain region (dentate gyrus) with a different learning paradigm (novel environment exposure). Our current understanding of the mechanism of learning is limited to appreciating the translatability of their gene list to our study. It is possible that there is a different, though maybe overlapping, set of genes regulating reactivation potential in M1. Nevertheless, with the gene list available to us, we found a reactivation signal in C3 that appears to be supported by DE and GO analyses. Besides the reactivation score, we found other lines of evidence linking C3 with reactivation, including dendritic spine morphogenesis that is preferably targeted to previously activated spines during motor learning in M1 (Hayashi-Takagi et al., 2015) and the enrichment by training (Cao et al., 2015). Taken together, the upregulation of genes for mRNA processing machinery and ubiquitin-proteasome processes in C3 neurons might support the late LTP-associated proteins that enlarge and maintain the previously potentiated and/or newly formed synapses, resulting in a reliable reactivation of previously activated neuronal ensemble. That is, C3 may represent motor engram cells with transcriptional responses that support the consolidation of learning. In line with this idea, the fraction of C3 neurons in our train group (∼25%) closely matches the previously observed fraction of movement-related neurons in the same lever-press task (20%–30%).
Another cluster that is unique to our dataset is C5. However, this cluster was equally represented in both the train and control groups. An intriguing question is which experimental manipulations in our study elicit transcriptional responses common to both groups, revealing C5. A prominent candidate is water restriction, which motivates mice to explore to find water sources. Water deprivation intensifying thirst promotes locomotor activity via orexin neurons in the hypothalamus that densely project to monoaminergic nuclei including LC, suggesting that it can stimulate NE neurons (Tsunematsu et al., 2008). Water restriction is also known to induce large-scale, global changes in cortical activity (Allen et al., 2019). However, the mechanisms underlying such a massive effect remain largely unknown. Our finding raises the possibility that thirst-driven global changes might be mediated by C5-like neurons across the cortex via the NE neuromodulatory system. Alternatively and/or additionally, C5 could be explained by differences in circadian rhythm; our experiments were conducted during the dark cycle to facilitate wakeful engagement in the task, while AIBS experiments were performed in the light cycle (Hongkui Zeng, personal communication; Yao et al., 2021, 2023). Circadian changes in gene transcription are well established, are highly conserved across species, and can even be regulated by neuronal activity (Xu et al., 2021). The transcriptional state represented by C5 may therefore at least partially reflect circadian differences. Disambiguating the contribution of water restriction and circadian rhythm to this transcriptional state will therefore require additional experiments wherein these features are specifically controlled. For instance, a number of prominent marker genes in C5 may be examined using RNAscope in four controlled groups of samples: (1) water-deprived and collected during the light cycle, (2) water-deprived and collected during the dark cycle, (3) nondeprived and collected during the light cycle, and (4) nondeprived and collected during the dark cycle.
Taken together, we found different clusters in the EDG space that represent various activity levels and activity-dependent plasticity phenotypes, including low-IEG baseline (C0), recent activity and LTP/LTD (C2), and reactivation (C3). The single-cell trajectory analysis indicates that individual neurons move along these different transcriptional states as they undergo different phases of neuronal plasticity. Furthermore, we found a cluster (C5) that may promote plasticity in M1 by increasing exploratory neural activity in response to the NE-mediated arousal signal. State transitions between clusters are unlikely to be an automatic program; instead, each transition must be regulated by experience and internal state, and only selective neurons in a given state must move to the next depending on their functional relation to the experience. Certain functional relations may also render some neurons prone to return to the baseline state during or after learning, although such transition is missing in the current trajectory analysis.
Our interpretations thus far rely on GO analysis and existing knowledge, the validation of which requires refined experiments for direct testing in the future. Nonetheless, our data demonstrate the emergence of diverse yet classifiable gene expression patterns within L2/3 neurons during long-term motor learning and provide a comprehensive set of coregulated genes associated with each cluster. We focused on L2/3 neurons in our dataset, leaving many other cell types awaiting exploration. The extensive time required for analysis and interpretation of results prevented in-depth exploration of the other cell types for the current study; however, sequencing data are publicly available for further studies.
In sum, we conclude that learning-inducing experiences like thirst-driven motor skill training induce a transient change in the transcriptional landscape of the same cell type neurons, as individual neurons visit divergent transcriptional states that underlie distinct phenotypes of neural activity and plasticity. The gene expression profiles we inferred for each plasticity phenotype will help future studies to delineate the molecular mechanisms that drive specific phenotypes and state transitions. In particular, investigating the causal role of the genes that are characteristic of each transcriptional state will illuminate the mechanisms regulating neural plasticity and thus learning.
Footnotes
Author contributions: R.G., T.K., Y.Z., and E.J.H. designed research; K.L.D., N.G.H., S.M., C.M.R., J.R.H., J.W., R.G., T.K., Y.Z., and E.J.H. performed research; C.M.R. contributed unpublished reagents/analytic tools; K.L.D., N.G.H., S.M., J.W., Y.Z., and E.J.H. analyzed data; K.L.D., N.G.H., T.K., Y.Z., and E.J.H. wrote the paper. The authors declare no competing financial interests.
We thank Dr. Hanquing Liu for his technical assistance with single-cell RNA sequencing datasets published by the BRAIN Initiative Cell Census Network. N.G.H. was supported by NIH (K99 NS114175). S.M. and R.G. were supported by NSF (1940162). T.K. was supported by NSF (1940202) and NIH (R01 NS091010, R21 NS112750, R01 NS125298). Y.Z. and K.L.D. were supported by NSF (1939992). E.J.H. was supported by the Alfred P. Sloan Foundation.
↵* K.L.D. and N.G.H. are co-first authors.
The authors declare no competing financial interests.
- Correspondence should be addressed to Eun Jung Hwang at eunjung.hwang{at}rosalindfranklin.edu or Ying Zhang at yingzhang{at}uri.edu or Takaki Komiyama at tkomiyama{at}ucsd.edu.
