The Genomic Architecture of Pregnancy-Associated Plasticity in the Maternal Mouse Hippocampus

Abstract Pregnancy is associated with extraordinary plasticity in the maternal brain. Studies in humans and other mammals suggest extensive structural and functional remodeling of the female brain during and after pregnancy. However, we understand remarkably little about the molecular underpinnings of this natural phenomenon. To gain insight into pregnancy-associated hippocampal plasticity, we performed single nucleus RNA sequencing (snRNA-seq) and snATAC-seq from the mouse hippocampus before, during, and after pregnancy. We identified cell type-specific transcriptional and epigenetic signatures associated with pregnancy and postpartum adaptation. In addition, we analyzed receptor-ligand interactions and transcription factor (TF) motifs that inform hippocampal cell type identity and provide evidence of pregnancy-associated adaption. In total, these data provide a unique resource of coupled transcriptional and epigenetic data across a dynamic time period in the mouse hippocampus and suggest opportunities for functional interrogation of hormone-mediated plasticity.


Introduction
Pregnancy is a period of remarkable physiological adaptation in the maternal brain. Pregnancy-associated neuroplasticity is observed across mammals, with variations in temporal and regional specificity (Leuner and Sabihi, 2016). The profound molecular and cellular changes in the maternal brain are presumed to promote offspring survival and alter maternal behavior (Brunton and Russell, 2008), but these changes also have implications for maternal mental health. Critical periods for plasticity, like pregnancy, are also periods of vulnerability for neurologic disease. Pregnant and postpartum women are at increased risk for psychiatric disorders, including depression, anxiety, and psychosis (Hillerer et al., 2014). Therefore, insight into physiologic pregnancy plasticity may also shed light on potential mechanisms for pregnancy-associated mental illness. Maternal brain plasticity is predominantly driven by dynamic hormone levels in the peripartum period, as pregnancy-associated neurologic changes have been reproduced in virgin female rodents through exogenous hormone treatment (Keyser-Marcus et al., 2001;Kinsley et al., 2006). Steroid hormones, which are in flux during and after pregnancy, markedly affect neurologic function via intracellular and membrane-bound steroid receptors expressed throughout the brain (Moraga-Amaro et al., 2018). Steroids such as progesterone regulate cell survival, synapse formation and dendritogenesis, whereas estrogens modulate synaptic plasticity and learning (Rossetti et al., 2016). Similarly, female sex steroids have notable effects in non-neuronal cells, such as promoting myelination, glutamate uptake in astrocytes, and the inflammatory action of microglia (Kipp et al., 2016). While the precise trajectory of hormone concentrations vary between rodents and humans, both progesterone and estradiol increase during gestation in mammals, ultimately peaking around parturition (Duarte-Guterman et al., 2019). Similarly, glucocorticoids follow an upward trend during pregnancy, and in mice, glucocorticoid levels increase dramatically from mid to late gestation (Solano and Arck, 2020).
There is a major gap in knowledge in the molecular underpinnings of adaptation in the brain during pregnancy. Given that pregnancy-associated plasticity is highly brain region-specific, we focus on the hippocampus because of its importance in postpartum maternal care behavior, as well as memory, mood, and stress resiliency (Duarte-Guterman et al., 2019). Neither the cell type-specific transcriptional nor cis-regulatory mechanisms underlying maternal brain plasticity have been previously explored. To begin to address this gap, we leveraged complementary single nuclei RNA-sequencing (snRNA-seq) and snATAC-seq of the maternal hippocampus before, during, and after pregnancy. We probe cell network ligand-receptor interactions during and after pregnancy and characterize enriched transcription factor (TF) motifs that define cell type identity and function in the murine hippocampus. This joint transcriptional and epigenetic atlas provides molecular support of pregnancy-associated structural remodeling of hippocampal architecture, as well as changes in neurogenesis and neurotransmission.

Materials and Methods
Tissue collection and generation of single-nuclei suspensions All animal procedures were performed in accordance with the Harvard Medical School animal care committee's regulations. Virgin female mice were euthanized at eight weeks of age. Pregnant mice were euthanized at embryonic day (E)18, and postpartum female mice were euthanized three weeks after parturition. The hippocampus was freshly dissected and flash frozen at À80°C. Frozen tissue was thawed in 500 ml buffer HB (0.25 M sucrose, 25 mM KCl, 5 mM MgCl2, 20 mM Tricine-KOH pH 7.8, 0.15 mM spermine tetrahydrochloride, 0.5 mM spermidine trihydrochloride, and 1 mM DTT). The tissue was transferred to a 2 ml dounce; 500 ml 5% IGEPAL CA-630 (Sigma) and 1 ml HB were added to the tissue, and the tissue was homogenized with a tight pestle 10-15 times. The sample was transferred to a 15 ml tube and total solution brought to 9 ml with HB. In a Corex tube (Fisher Scientific), 1 ml 30% iodixanol layered on top of 1 ml 40% iodixanol. The 9-ml sample was layered on top of the iodixanol cushion. The sample was spun at 10,000 Â g for 18 min; 1 ml of sample at the 30-40% iodixanol interface was collected. After counting nuclei with a hemocytometer, the sample was diluted to 100,000 nuclei/ml with 30% iodixanol (with Rnasin) and subjected to single nuclear droplet encapsulation with inDrop (Klein et al., 2015).

snRNA-seq and analysis
Individual nuclei were captured and barcoded using the inDrop platform as previously described (Klein et al., 2015). Briefly, single-cell suspensions were fed into a microfluidic device that packaged the cells with barcoded hydrogel microspheres and reverse transcriptase/lysis reagents. After cell encapsulation, primers were photo-released by UV exposure. Two technical replicate libraries were collected for each sample. Indexed libraries were pooled and sequenced on a Nextseq 500 (Illumina) to an average depth of ;30,000 reads/nucleus. Sequencing data were aligned to the genome and processed according to a previously published pipeline (https://github.com/indrops/ indrops). Unique molecular identifiers (UMIs) were used to link sequence reads back to individual captured molecules. All steps of the pipeline were run using default parameters.
All data processing was performed using Seurat (v. 4.0). Raw data were log normalized and scaled using default values. Variable features and principal components were then calculated using default values. Uniform manifold approximation and projection (UMAP) and t-distributed stochastic neighbor embedding (tSNE) dimensionality reductions were performed with default values. Batch effects were not evident in the dimensionality reductions, and therefore, the data were then analyzed as-is without further corrections.

Identification of cell types in snRNA-seq
We used scPred (Alquicira-Hernandez et al., 2019) to train radial basis kernel support vector machines (SVMs) in a one versus all fashion using data deposited to Allen Brain Atlas (Yao et al., 2021). For consistency, we only used cells from the hippocampus (212,939 cells) to train the predictor models. Before training, the training data were independently preprocessed in the same fashion as above. For each cell type annotated in the Allen Brain Atlas, we trained a separate one versus all SVM model with a Radial Basis Function kernel. The training data were split 80/20; 80% of the training data were centered, scaled, and normalized in the same fashion as our dataset. Models were trained using a 5-fold cross validation approach and only models with area under the receiver operating characteristic (AUROC) . 0.95 were used. Using the remaining 20% of the training dataset, we then predicted the cell types of each cell in our dataset. Only cells with in-class prediction (e.g., CA1 vs all others) probability .0.95 were retained. We manually confirmed cell types based on known marker gene expression. Identified cell types (as some cell types corresponded to multiple clusters) with ,100 cells per condition were not used for further analysis.

Differential expression and ligand-receptor analysis
For each identified cell type, we performed differential expression analysis using Seurat with MAST method between different conditions (e.g., virgin CA3 vs E18 CA3). We used Benjamini-Hochberg corrected p-values for significance. Genes with ,5% false discovery rate (FDR) met statistical criteria for significant differential expression. Ligand-receptor analysis was performed using CellChat (version 1.1.3) with default parameters (Jin et al., 2021).

snATAC-seq
Tissue was processed for single nucleus suspensions as described above. Nuclei were processed using the 10Â Chromium scATAC-seq kit (v1.0) per manufacturer's instructions. Libraries were sequenced to an average depth of ;50,000 reads/nucleus. Sequencing reads were demultiplexed and aligned to the mm10 genome using CellRanger with default parameters. Resulting outputs were then processed using Signac (v 1.5; Stuart et al., 2021). Low quality cells were removed using filters peak region fragments ,100,000, percent reads in peaks .40, blacklist site ratio ,0.025, nucleosome signal ,4, and transcriptional start site (TSS) enrichment score .2. Remaining preprocessing and dimensionality reduction was performed according to Signac recommendations using default parameters.
Cell types were determined using the snRNA-seq data labels based on Signac recommendations (Stuart et al., 2021). We used FindTransferAnchors and TransferData functions to annotate cell types. We then again checked for consistency between unsupervised clustering and label transfer. Similar to snRNA-seq data, clusters with high cell type heterogeneity were discarded and cell types were assigned to clusters based on the majority of the labels. Cell types with ,100 cells per condition were not analyzed further.

Differential accessibility and TF motif analysis
Using cell type labels described above, we performed differential accessibility analysis using Signac's FindMarkers function with log fold change threshold 0.1 and n_count_peaks as a latent variable. FDR ,5% was used as the statistical significance threshold. Motif analysis was done using RunChromVAR function in Signac with default parameters. For differential motif activity, we used log2 fold change threshold of 0.1 with a 0 pseudocount. For statistical significance, we used an FDR threshold of 5%.

Integrative analysis with ArchR
In order to perform co-accessibility and peak-to-gene linkage analysis, we processed snATAC-seq data with ArchR (Granja et al., 2021). Data were processed using default parameters. We used the addCoAccessibility and getCoAccessibility functions to generate co-accessibility data frames and figures. We performed peak-to-gene linkage analysis using the addPeak2GeneLinks function, and the plotPeak2GeneHeatmap function was used to generate a heatmap of peak-to-gene results.

Data accessibility
Raw and processed sequencing data have been deposited in the Gene Expression Omnibus (GEO) under accession number GSE198447.

Results
Cell type-specific transcriptomic landscape of the maternal hippocampus To identify cell type-specific programs governing plasticity during pregnancy, we performed droplet-based snRNA-seq from the female hippocampus under three conditions: eight-week-old virgin female mice, E18 pregnant female mice, and P21 (postpartum) female mice (Fig.  1a). These time points were selected to span the continuum of hormonal shifts that occur before, during, and after pregnancy. We analyzed three biological replicates at each time point, for a total of nine independent samples. Virgin mice were synchronized in their estrous cycle before tissue collection.
Nuclei were isolated from frozen tissue (see Materials and Methods) and processed using the inDrop platform (Klein et al., 2015). After computational filtering of low quality and doublet nuclei, we profiled a total of 40,652 nuclei by snRNA-seq. Quality control metrics are displayed in Figure 1b. We used the R package Seurat to perform unsupervised graph-based clustering of all nuclei, followed by cluster visualization with UMAP, a nonlinear dimensional reduction algorithm. To aid in the confident identification of cell types, we leveraged previously published high quality datasets from Allen Brain Atlas (Alquicira-Hernandez et al., 2019). With these expert labeled datasets, we used scPred (Alquicira-Hernandez et al., 2019) to train one versus all SVM models with radial kernels to identify cell types. Models were trained using a 5fold cross-validation approach. Using each of the models, we then predicted the cell types of each cell in our dataset. Only cells with in-class prediction (e.g., CA1 pyramidal neurons versus all others) probability .0.95 were retained. Complete AUROC values are provided in Extended Data Figure 1-1.
In total, we identified clusters corresponding to major hippocampal cell types: glutamatergic (granule cells, pyramidal cells), GABAergic, and non-neuronal (microglia, astrocytes, oligodendrocytes; Fig. 1c). Expression of canonical cell type marker genes affirmed cell identifies (Fig. 1d). Cell cluster assignment did not globally differ between conditions, suggesting that no new cell types were detected between the three conditions (Fig. 1e). Hippocampal cell types were subclustered into glutamatergic, GABAergic, and non-neuronal groups ( Fig. 2a-d).
We implemented differential gene expression analysis to investigate transcriptomic changes across pregnancy Figure 1. Single-cell investigation of the female hippocampus before, during, and after pregnancy. a, Experimental schematic of snRNA-seq and snATAC-seq of the female hippocampus in virgin, pregnant (E18), and postpartum mice. b, Quality control (QC) plots of UMIs and genes per nucleus in nuclei included in the final snRNA-seq dataset. c, UMAP of snRNA-seq data with major cell types identified. Extended Data Figure 1-1 provides AUROC, sensitivity, and specificity data for SVM model for cell type identification in snRNA-seq. d, Violin plots of canonical marker gene expression across cell types, as measured by snRNA-seq. e, UMAPs of nuclei included in snRNA-seq dataset, separated by condition.
Research Article: New Research within each cell type. Statistical significance was assigned to genes with a FDR , 5%. All differential gene expression output is summarized by cell type in Extended Data Figure 3-1. We identified cell type-specific changes in gene expression across the time course of pregnancy and postpartum period (Fig. 3a). We performed gene ontology (GO) analysis of statistically significant differentially expressed genes between timepoints to identify functional pathways that exhibit peripartum cell type-specific changes . Regulation of the actin cytoskeleton (including the profilin complex), endocytosis, trans-synaptic signaling, and the trafficking of neurotransmitter receptors were commonly identified pathways, consistent with hippocampal plasticity during this critical window.
We then focused our interpretation of the snRNA-seq data on three specific themes, neurogenesis, neurotransmission, and structural remodeling, that are known to be perturbed in the female brain during and after pregnancy (Galea et al., 2014;Hillerer et al., 2014;Duarte-Guterman et al., 2019).

Neurogenesis
The effect of pregnancy on the birth of new neurons demonstrates temporal, region, and species specificity (Leuner and Sabihi, 2016). In the maternal mouse brain, hippocampal neurogenesis is decreased at E18 relative to age-matched nonpregnant controls (Kim et al., 2010). This reduction in hippocampal neurogenesis is sustained following parturition, with the number of neural precursor cells (NPCs) returning to prepregnancy levels three weeks postpartum (Rolls et al., 2008). Thus, fluctuations in neurogenesis might contribute to the reductions in hippocampus volume that occur in pregnancy in both rats and humans (Galea et al., 2000;Hoekzema et al., 2017).
The dentate gyrus (DG) is the site of hippocampal neurogenesis. We probed our snRNA-seq dataset of the hippocampus and identified gene expression patterns in DG Figure 2. Subclustering analysis of snRNA-seq of the female hippocampus before, during, and after pregnancy. a, UMAP of glutamatergic cells, with DG granule cells, CA1 pyramidal cells, and CA3 pyramidal cells labeled by color. b, Violin plots of canonical interneuron marker genes, as measured by snRNA-seq. c, UMAP of GABAergic cells, with interneuron subtypes labeled by color. d, UMAP of non-neuronal cells, with subtypes labeled by color.
Research Article: New Research granule cells consistent with decreased neurogenesis during late pregnancy. Specifically, we found reduced expression of the transcriptional activator Sox11 at E18 compared with virgin mice (Fig. 4a). Sox11 is highly expressed in regions of adult neurogenesis and is sufficient to induce the differentiation of NPCs to immature neurons (Haslinger et al., 2009). Notably, the expression of Neurite outgrowth inhibitor (Nogo), also known as Reticulon 4, is increased at E18 relative to virgin mice in DG granule cells (Fig. 4a). This is consistent with the role of Nogo and its receptor NgR1 as negative regulators of NPC proliferation (Rolando et al., 2012). We found that the expression of other reticulon family members, including Reticulon (Rtn) 1 and 3, also exhibit higher expression in DG granule cells during pregnancy when compared with virgin levels (Fig.  4a). Reticulon family members are known to promote apoptosis (Yang and Strittmatter, 2007) which is notable given the importance of apoptosis in shaping cell fate during neurogenesis (Tashiro et al., 2006;Sierra et al., 2010).
We also identified coordinate changes in insulin-like growth factor signaling in the hippocampus during pregnancy. Insulin-like growth factor 2 (Igf2) expression in CA1 and CA3 pyramidal neurons, as well as in GABAergic interneurons, is lower during pregnancy relative to virgin controls ( Fig. 4b-d). In parallel, insulin-like growth factor 1 receptor (Igf1r) expression is also lower during pregnancy in the DG (Fig. 4a). IGF2 primarily binds to IGF1R and regulates cell proliferation during adult hippocampal neurogenesis (Bracko et al., 2012;Ferrón et al., 2015). Specifically, IGF2 in NPCs stimulates self-renewal and increases mRNA expression of the canonical stem cell TFs Oct4 and Sox1 (Ziegler et al., 2012).
We further identified lower expression of cell cycle regulators in DG granule cells of pregnant mice relative to . Differential gene expression analysis to identify candidate regulators of peripartum plasticity. a, Volcano plot of differential gene expression between virgin and E18 condition in DG granule cells. Statistically significant genes (FDR , 5%) are colored, with blue indicating lower expression and pink indicating higher expression at E18. b, Volcano plot of differential gene expression between virgin and E18 condition in CA1. Statistically significant genes (FDR , 5%) are colored, with blue indicating lower expression and pink indicating higher expression at E18. c, Volcano plot of differential gene expression between virgin and E18 condition in CA3. Statistically significant genes (FDR , 5%) are colored, with blue indicating lower expression and pink indicating higher expression at E18. d, Volcano plot of differential gene expression between virgin and E18 condition in GABAergic cells. Statistically significant genes (FDR , 5%) are colored, with blue indicating lower expression and pink indicating higher expression at E18. e, Volcano plot of differential gene expression between virgin and E18 condition in astrocytes. Statistically significant genes (FDR , 5%) are colored, with blue indicating lower expression and pink indicating higher expression at E18. virgin controls, including three subunits of the anaphasepromoting complex (Anapc1, Anapc5, and Anapc7) and ataxia telangiectasia and Rad3 related (Atr; Fig. 4a). In mouse embryos, ANAPC facilitates neurogenesis, and ATR governs the transition from G 2 to M in the mitotic cycle of NPCs (Delgado-Esteban et al., 2013;Enriquez-Rios et al., 2017).

Neurotransmission
In both humans and rodents, memory function is profoundly altered in the peripartum and postpartum periods (Pawluski et al., 2016). Given that memories are stored through enhanced connectivity between groups of neurons, changes in neurotransmission may underlie modifications in memory by affecting synaptic strength (Josselyn et al., 2015). For example, in rats, motherhood induces long-term potentiation (LTP) with higher NMDA and non-NMDA-mediated hippocampal neurotransmission relative to virgin controls (Lemaire et al., 2006). Accordingly, we analyzed our snRNA-seq dataset for molecular mediators of neurotransmission. We found that several conventional protein kinase C isozymes, including a, b , and g , have increased expression in DG granule cells and CA1 pyramidal neurons during pregnancy compared with virgin mice (Fig. 4a,b). PKCa is the requisite isozyme for structural LTP, and PKCg is likely connected to LTP through its phosphorylation of the NMDA receptor (NMDAR) subunit NR1 (Gomis-González et al., 2021). We identified further changes in gene expression related to LTP in the hippocampus during pregnancy. Specifically, at E18 relative to virgin, we noted increased expression of the A and b isoforms of calmodulin-dependent protein kinase II (CaMKII) in DG granule cells and CA1 pyramidal neurons (Fig. 4a,b). This serine/threonine-protein kinase is a vital component of LTP by modulation of AMPAR insertion and conductance (Lisman et al., 2012). Consistent with increased CaMKII, we identified the increased expression of Synapsin 1 (Syn1) and Synapsin 2 (Syn2) in DG granule cells at E18 relative to virgin (Fig. 4a). These phosphoproteins are subject to CAMKII activation and have many neuronal functions, including linking vesicles to the cytoskeleton at the presynaptic terminal (Cesca et al., 2010). Syn1 deficient mice have impaired glutamate release (L Li et al., 1995), and likewise, GABA neurotransmission is distorted in Syn2 deficient mice (Medrihan et al., 2013). We also found increased expression of Ras-related protein Rab-3A (RAB3A) and Rab-3C (RAB3C) in DG granule cells at E18 relative to virgin (Fig. 4a), consistent with changes in the presynaptic terminal. These small GTP binding proteins control calcium-induced exocytosis (Schlüter et al., 2004) and BDNF-dependent glutamate release (Simsek-Duran and Lonart, 2008).
Non-neuronal cells, particularly astrocytes, are essential regulators of neuronal communication because of their ability to sense and respond to neurotransmitter release (Hwang et al., 2021). In astrocytes, we found increased expression of glutamate receptor ionotropic AMPA type subunit 1 (Gria1) in pregnant mice relative to virgin controls (Fig. 4e). Astrocytes regulate glutamate homeostasis through glutamate uptake, enabling them to modulate neighboring neuronal activity (Droste et al., 2017;Mahmoud et al., 2019).

Structural remodeling
The maternal neural circuit architecture exhibits substantial structural plasticity during pregnancy (Hillerer et al., 2014). Notably, membrane-bound estrogen receptors are prevalent in the mammalian hippocampus and act through a defined mechanism to mediate dendritic spine as well as synapse formation (Hara et al., 2015). In our data, we identified differential expression of genes implicated in structural remodeling during and after pregnancy.
Strikingly, we found numerous candidates pertinent to actin restructuring to be more highly expressed during pregnancy. This is relevant because dendritic spines have a complex cytoskeleton of branched and linear actin filaments (Basu and Lamprecht, 2018). Consequently, dendritic size and morphology regulation, such as in response to glutamate stimulation, requires tightly controlled actin polymerization and depolymerization (Matsuzaki et al., 2004;Borovac et al., 2018). The expression of Actin Depolymerizing Factor (ADF) is higher during pregnancy in DG granule cells compared with virgin controls (Fig. 4a). The ADF/Cofilin family governs the structural plasticity of dendritic spines, and active ADF/Cofilin enables increases in spine size during LTP and the diminishing of spines associated with LTD (Rust, 2015). Similar to ADF/Cofilin, the Rho GTPase family modifies dendrite morphology through actin skeleton dynamics (Murakoshi et al., 2011). In DG granule cells, we noted increased expression of RhoB and Rac1 at E18 relative to virgin mice (Fig. 4a). Consistent with this RhoGTPase mechanism, we identified increased expression of p21 activated kinase (Pak1) at E18 relative to virgin in DG granule cells (Fig. 4a). Pak1 is a notable downstream serine/threonine kinase effector of the Rho GTPase family, and its phosphorylation is associated with neurite organization and outgrowth (Rashid et al., 2001;Nikoli c, 2008).
Moreover, in CA1 pyramidal neurons and DG granule cells during pregnancy, we found increased expression of the postsynaptic density protein 95 (PSD-95) and synapse-associated protein 102 (SAP-102) compared with postpartum mice (Fig. 4a,b). Together, these membraneassociated guanylate kinases (MAGUKs) are critical structural components of the postsynaptic density (PSD) and serve to regulate its size (Zheng et al., 2011). Interestingly, the translation of PSD-95 mRNA is stimulated by estrogen signaling, and this process is thought to partly underlie estrogen's positive effect on synaptogenesis (Akama and McEwen, 2003).

Ligand-Receptor analysis of the maternal hippocampus
Cells in the brain communicate through soluble and membrane-bound components. To map ligand-receptor networks in the maternal hippocampus at single cell resolution, we used the open-source R toolkit CellChat (Jin et al., 2021). We identified several signaling networks in the hippocampus that are perturbed during or after pregnancy. The CX3C signaling pathway is upregulated during pregnancy relative to virgin conditions (Fig. 5a). CX3C signaling regulates neuronal network maintenance and is dependent on interactions between the transmembrane chemokine CX3CL1 and its G-protein coupled receptor CX3CR1, found only in microglia (Cardona et al., 2006;Paolicelli et al., 2011;Pawelec et al., 2020). The dominant interaction occurs between the CX3CR1 on microglia and fractalkine secreted by somatostatin interneurons. Thus, we posit that a method of synapse refinement during pregnancy involves microglial phagocytosis driven by the CX3C pathway. Alternatively, CX3C signaling might underlie the increased oligodendrocyte differentiation during pregnancy (Gregg, 2009), as interneurons mediate developmental oligodendrogenesis through secretion of CX3CL1 (Voronova et al., 2017).
We also found that the pleiotrophin (PTN) network is upregulated during pregnancy compared with both virgin and postpartum conditions (Fig. 5a). Interactions between PTN with PTPRZ receptors on oligodendrocyte precursor cells (OPCs) stimulate differentiation to oligodendrocytes, and therefore the capacity for myelination might be increased during pregnancy (Tanga et al., 2019). This is consistent with murine lineage-tracing experiments that indicate a 74% increase in oligodendrocyte formation in the corpus callosum at E18 relative to virgin controls (Gregg et al., 2007).
Three signaling networks (CXCL12/Ackr3, IGF1/IGF1r, Wnt5a/Fzd6) were significantly downregulated in pregnancy (E18) compared with virgin (Fig. 5b). Signaling between chemokine CXCL12 and its receptor ACKR3 on endothelial cells declined in pregnancy but increased postpartum. CXCL12/ACKR3 signaling is important in directing cell migration and mediating angiogenesis (Guyon, 2014;Quinn et al., 2018). Similarly, IGF1/IGF1R signaling declined during pregnancy, which is notable given that IGF1 deficiency in the hippocampus is linked to depression (Torres Aleman, 2005; Mitschelen et al., 2011) Wnt5a signaling in endothelial cells also declined during pregnancy compared with virgin (Fig. 5b). Wnt ligands are important for angiogenesis, tight junction Figure 5. Ligand-receptor analysis of snRNA-seq reveals dynamic signaling networks during pregnancy. a, Chord diagram of cell type-specific signaling pathways that are upregulated at E18 compared with virgin. Cell type identity of the "sending" cell is indicated in the outermost edge of the diagram, while cell identity of the "receiving" cell is indicated by the color second from the outside. Colored arrows in the diagram indicate specific ligand/receptor pathways. b, Chord diagram of cell type-specific signaling pathways that are downregulated at E18 compared with virgin. Cell type identity of the "sending" cell is indicated in the outermost edge of the diagram, while cell identity of the "receiving" cell is indicated by the color second from the outside. Colored arrows in the diagram indicate specific ligand/receptor pathways. c, Chord diagram of cell type-specific signaling pathways that are upregulated postpartum compared with pregnancy. Cell type identity of the "sending" cell is indicated in the outermost edge of the diagram, while cell identity of the "receiving" cell is indicated by the color second from the outside. Colored arrows in the diagram indicate specific ligand/receptor pathways. formation, and blood brain barrier integrity, which is notable given changes in permeability of the blood brain barrier during pregnancy (Daneman et al., 2009;Cipolla et al., 2011).
The Prosaposin (PSAP)/GPR37 signaling network increases postpartum compared with pregnancy (E18; Fig. 5c). PSAP and its receptors (GPR37 and GPR37l1) are expressed in astrocytes, OPCs, and oligodendrocytes. GPR37l1 signaling inhibits glutamate uptake in astrocytes and reduces neuronal NMDAR activity (Jolly et al., 2018). Loss of GPR37 in oligodendrocytes accelerates differentiation and hypermyelination, and therefore, lower levels of PSAP/GPR37 signaling may support increased myelination during pregnancy (Yang et al., 2016;Smith et al., 2017) Cis-regulatory landscape of the maternal hippocampus Chromatin organization governs transcription factor (TF) binding and transcriptional regulation, which are necessary for the dynamic adaptation of neural circuits. To complement the cell type-specific transcriptional atlas of the maternal hippocampus, we performed snATACseq at identical timepoints to those described above. Our snATACseq dataset includes 38,123 nuclei that pass quality filtering (Fig. 6a-d), with an average of 216,990 unique nuclear genome fragments. The proportion of accessible regions within genes, promoters, or distal regulatory sites was 69.6%, 0.1%, 30.3%, respectively. We assigned cell type identity to snATAC-seq clusters using anchoring and label transfer from our snRNA-seq data via Signac (Stuart et al., 2021). Label transfer and integration of snRNA-seq and snATAC-seq demonstrated strong agreement between computational and supervised (marker gene-driven) linkage of the datasets (Fig. 6e,f). To further correlate transcriptional and chromatin dynamics within cell types, we analyzed the combined datasets using ArchR Peak2GeneLinkage (Granja et al., 2021). We identified strong cell type-specific modules of gene-chromatin correlations, as evidenced by the heatmap in Figure 7a, suggesting distinct gene regulatory interactions. Figure 6. Clustering and quality control metrics of snATAC-seq of the murine hippocampus before, during, and after pregnancy. a, UMAP of snATAC-seq data with major cell types identified. b, TSS enrichment score, defined as the ratio of fragments centered at the TSS to fragments in TSS-flanking regions, for all cells included in snATAC-seq analysis, grouped by condition. c, Percentage of reads in peaks for all cells included in snATAC-seq analysis, grouped by condition. d, Ratio of reads in genomic blacklist regions as defined by the ENCODE project, representing reads that are often associated with artefactual signal. e, Agreement between predicted snATACseq cell types, as identified by Seurat label transfer, and supervised (manual) cell type annotation based on the gene activity matrix of canonical cell type-specific marker genes. The color of the squares indicates the fraction of cells in which there is agreement in cell type identity. f, Prediction score of label transfer between snRNA-seq and snATAC-seq, comparing nuclei with "true" identity assignment (concordance between Seurat-based label transfer analysis and manual annotation based on marker genes) and "false" identity assignment (disagreement between Seurat-based label transfer analysis and manual annotation based on marker genes).
First, given limited existing data on the chromatin landscape of the murine hippocampus, we identified chromatin accessibility signatures for distinct cell types in the mouse hippocampus across all timepoints in a pooled analysis (virgin, E18, postpartum). In the GABAergic neuronal population, we noted the marked increase in the accessibility of the Nuclear Respiratory factor 1 (NRF-1) motif relative to glutamatergic neurons (Fig. 7b). This TF serves as an activity-dependent activator for the transcription of the b 1 subunit of the GABA-A receptor (Z Li et al., 2018). Conversely, we identified the enhanced accessibility of the Neurogenin 1 (Neurog1) motif in the glutamatergic population relative to GABAergic neurons (Fig. 7c). Neurog1 is a critical TF for glutamatergic fate specification (Reyes et al., 2008;Shaker et al., 2012). In astrocytes relative to oligodendrocytes, there was increased accessibility of the Interferon regulatory factor 1 (IRF1) motif (Fig. 7d); IRF1 restricts the differentiation of OPCs to oligodendrocytes (Tanner et al., 2011). We also found the accessibility of the Sex determining region Y-box 2 (SOX2) motif to be greater in oligodendrocytes (Fig. 7e). In murine brain development, the proliferation and differentiation of oligodendrocytes are upregulated by SOX2 expression, and in adulthood, SOX2 facilitates remyelination following white matter damage (S Zhang et al., 2018). These results, while largely expected based on known mechanisms of TF activity in brain cell types, support the fidelity of chromatin accessibility signatures in our dataset.
We then constructed chromatin co-accessibility networks for canonical cell type identity genes for the murine hippocampus, again across all timepoints in a pooled analysis (virgin, E18, postpartum). Co-accessibility is a correlation in chromatin accessibility between two peaks across many single cells. We used ArchR to map proximal and distal chromatin regulatory elements associated with target genes. We identified chromatin co-accessibility in regulatory elements linking Olig1 and Olig2 loci, which are genes predominantly expression in oligodendrocyte lineage cells (Fig. 7f). We also mapped the complex cis-regulatory landscape surrounding the Cx3cr1 locus, which is associated with microglial identity (Fig. 7g).
We next asked whether there were cell type-specific changes in chromatin accessibility consistent with pregnancy-associated neuroplasticity. We performed pairwise differential accessibility analysis between conditions (virgin, E18, postpartum) within cell types, and these data are provided in Extended Data Figure 8-1. Broadly, using the peaks with significant changes in accessibility during pregnancy (virgin vs E18), we used GO analysis to identify functional pathways that exhibit dynamic cell type-specific changes (DG granule cells, CA1 pyramidal neurons, and astrocytes; Fig. 8a-c; Extended Data Fig. 8-2). The Wnt/ b -catenin signaling, which is critical in regulating neurogenesis (Arredondo et al., 2020), synapse strength (Okuda et al., 2007), and dendritic structure (Heppt et al., 2020), was a common pathway identified across several cell types and is known to be enhanced by estrogen (Varea et al., 2009). Next, we focused on three themes, similar to the snRNAseq analysis: neurogenesis, neurotransmission, and structural remodeling.

Neurogenesis
Cell type-specific chromatin accessibility analysis corroborated our transcriptional findings associated with impaired neurogenesis during pregnancy. Specifically, the accessibility of the Sox4 locus in DG granule cells is decreased during pregnancy (E18) compared with virgin (Fig. 8d). Sox4 is a member of the SoxC group within the SOX TF family and facilitates neurogenesis both in vitro and in vivo (Mu et al., 2012). We further found the accessibility of Erb-B2 receptor tyrosine kinase 4 (Erbb4) locus to be lower during pregnancy (E18) compared with virgin (Fig. 8d). In the subventricular zone, another neurogenic niche, Erbb4 is necessary for the survival, migration and organisation of NPCs (Anton et al., 2004;H Zhang et al., 2018).
We inferred TF activity through analysis of differential motif accessibility and found evidence to support a decline in neurogenesis during pregnancy. In DG granule cells, we found the decreased accessibility of several FOS and JUN motifs in pregnant mice (E18) relative to virgin controls (Fig. 8f). These proteins are members of the Activator protein-1 (AP1) TF family, broadly regulating cell proliferation and differentiation (Shaulian and Karin, 2002). Specific to neurogenesis, the AP1 complex likely regulates NPC proliferation, and Fos is sufficient to rescue NPC self-renewal in SOX2 knock-out (KO) mice (Pagin et al., 2021). Moreover, there is increased accessibility of the Etv1 motif in DG granule cells of pregnant mice (E18) when compared with virgin controls (Fig. 8f). In vitro, ETV1 represses cell proliferation genes in murine cerebellar granule cells and promotes mature neural circuit formation (Okazawa et al., 2016).

Structural remodeling
We identified changes in cell type-specific chromatin accessibility to support the structural remodeling of the neural circuit associated with pregnancy. In particular, in CA3 pyramidal neurons, we discovered increased accessibility of the Nectin Cell Adhesion Molecule 3 (Nectin3) continued Figure 7. Analysis of TF motifs and differential accessibility from snATAC-seq of the murine hippocampus. a, Integrated analysis of snRNA-seq and snATAC-seq with ArchR peak-to-gene linkage. The heatmap displays k-means clustered genes or chromatin loci as rows. Colored columns indicate cell types. b, Motif feature plot and position weight matrix (PWM) for the TF Nrg1, enriched in GABAergic cells. c, Motif feature plot and PWM for the TF Neurog1, enriched in glutamatergic cells. d, Motif feature plot and PWM for the TF IRF1, enriched in microglia. e, Motif feature plot and PWM for the TF Sox2, enriched in mature oligodendrocytes. f, ArchR-based chromatin co-accessibility analysis of the Olig1 and Olig2 loci. Top panels display genomic tracks of chromatin accessibility, organized by cell type. Bottom panels display location and linkage of co-accessible peaks near the Olig1 and Olig2 loci. g, ArchR-based chromatin co-accessibility analysis of the Cx3cr1 locus. Top panels display genomic tracks of chromatin accessibility, organized by cell type. Bottom panels display location and linkage of co-accessible peaks near the Cx3cr1 locus. locus during pregnancy (E18) relative to virgin conditions. This synaptic cell adhesion molecule is located at the postsynaptic terminal and coordinates the actin cytoskeleton and cadherin-catenin complex during synaptogenesis and synaptic remodeling (R Liu et al., 2019). Likewise, Shank2 is a postsynaptic terminal component and binding partner of PSD-95 (Brandstätter et al., 2004). By snRNA-seq, we found the increased expression of PSD-95 in CA1 pyramidal neurons during pregnancy (E18) compared with virgin. Interestingly, we also noted that the accessibility of the Shank2 locus is greater in CA1 pyramidal neurons during pregnancy (E18) compared with virgin (Fig. 7b). Akin to PSD-95, Shank2 is a postsynaptic scaffold protein localized to excitatory synapses and, along with several other components, organizes the location of metabotropic and ionotropic glutamate receptors (Monteiro and Feng, 2017). . Differential accessibility analysis of snATAC-seq of the murine hippocampus before, during, and after pregnancy. a, GO analysis of differentially accessible peaks (FDR , 5%) between virgin and E18 in DG granule cells. Extended Data Figure 8-1 provides complete results of differential chromatin accessibility analysis, and Extended Data Figure 8-2 provides complete results of GO analysis from snATAC-seq. b, GO analysis of differentially accessible peaks (FDR , 5%) between virgin and E18 in CA1 pyramidal cells. c, GO analysis of differentially accessible peaks (FDR , 5%) between virgin and E18 in astrocytes. d, MA plot (log2 fold change vs log2 mean expression) of differentially accessible peaks in DG comparing virgin versus E18, with select candidates (Erbb4, Sox4, and Grin2a) highlighted. e, MA plot (log2 fold change vs log2 mean expression) of differentially accessible peaks in CA1 comparing virgin versus E18, with select candidates (Grid1, Grid2, Shank2, Grin2a) highlighted. f, Position weight matrices (PWM) for TF motifs with significant changes at E18 compared with virgin. FDR values reflect comparisons in the DG for FOS, JUN, ETV1, ETV4, SP4, and ESR1. FDR values reflect comparisons in CA3 for CTCF and CA1 for ELK1.
With regard to TF motifs, we identified increased accessibility of the ETV4 motif in DG granule cells, as well as CA1 and CA3 pyramidal neurons, during pregnancy (E18) relative to virgin conditions (Fig. 8f). This TF likely mediates BDNF signaling and is highly expressed during hippocampal dendrite development, with ETV4 deletion resulting in decreased dendritic complexity and size (D Liu et al., 2016;Fontanet et al., 2018). We further identified increased accessibility of the SP4 TF motif in CA1 and CA3 pyramidal neurons during pregnancy (E18) compared with virgin (Fig. 7c). In rodents, SP4 has been found to modulate dendritic arborization and pruning as well as NMDAR trafficking (Ramos et al., 2007;Sun et al., 2015). Notably, the accessibility of the Estrogen receptor a (ESR1) motif is increased during pregnancy (E18) in DG granule cells and CA1 and CA3 pyramidal neurons compared with virgin (Fig. 8f). In the hippocampus, this nuclear receptor modulates a broad spectrum of plasticity, including dendritic spine density and the abundance of NMDARs (Mukai et al., 2010). This finding is consistent with a hormone-mediated model of pregnancy-associated neuroplasticity.

Neurotransmission
Changes in the chromatin accessibility landscape during pregnancy also corroborate pregnancy-associated adaptations in neurotransmission. In CA1 pyramidal neurons, there is increased accessibility of the loci for glutamate receptor d -1 and d -2 subunits (Grid1 and Grid2) during pregnancy (E18) compared with virgin (Fig. 8e). Both genes play a role in synaptic plasticity, albeit in different neuronal populations (Fossati et al., 2019) and Grid1, in particular, mediates hippocampal mGlu5 signaling (Suryavanshi et al., 2016). Furthermore, in CA1 pyramidal neurons and DG granule cells, we identified increased accessibility of the NMDAR subunit «À1 (Grin2a) locus ( Fig. 7b) in E18 compared with control. In CA3 pyramidal neurons, we identified the NMDAR subunit « À2 (Grin2b) locus as more accessible at E18 compared with virgin. These subunits are integral to synaptic plasticity, with Grin2a mutant mice experiencing reduced LTP in CA1 (Sakimura et al., 1995), and in Grin2b mutants, NMDARmediated LTD is eliminated (Brigman et al., 2010).
Analysis of differential TF motifs between pregnancy conditions revealed further evidence of pregnancy-associated modifications in neurotransmission. The motif for CCCTC-binding factor (CTCF) is more accessible in CA3 pyramidal neurons of the pregnant hippocampus (E18) when compared with virgin mice (Fig. 8f). CTCF-induced gene expression is thought to underlie adult hippocampal memory via chromatin reorganization, and CTCF KO mice show disrupted LTP (Sams et al., 2016). We also noted the enhanced accessibility of the ETS TF ELK1 in CA3 and CA1 pyramidal neurons during pregnancy (E18) when compared with virgin (Fig. 8f). In the hippocampus, ELK1 phosphorylation via the MAPK/ERK pathway has been associated with induction of LTD and is also thought to mediate LTP-dependent gene regulation (Davis et al., 2000;Besnard et al., 2011).
In summary, our multiomics characterization provides a detailed landscape of complementary layers of regulatory control governing peripartum plasticity and identify candidate regulators of hormone-mediated adaptation.

Discussion
Despite long-standing recognition of the extraordinary plasticity of the female brain during and after pregnancy, there is a shortage of molecular insights into this physiological adaptation. Hormone-driven plasticity during pregnancy has likely evolved to prepare females for childrearing and improve offspring survival (Brunton and Russell, 2008). Human neuroimaging studies have demonstrated pregnancy-associated decreases in gray matter volume that persists for at least two years postpartum (Hoekzema et al., 2017). These volume reductions are particularly prominent in regions implicated in maternal caregiving (Barha and Galea, 2017;Michaelian et al., 2019). Although the existing literature on humans is limited, pregnancy-associated cognitive changes are supported by self-reported symptoms of working memory deficits and general "mental fogginess" (Brown and Schaffir, 2019). Globally, .200 million women are pregnant each year, and it is therefore critical that we understand the mechanisms by which pregnancy shapes the female brain.
In this study, we leverage snRNA-seq and snATAC-seq to identify candidate transcriptional and epigenetic regulators of pregnancy-associated neuroplasticity. We focus on the hippocampus, as this region demonstrates substantial plasticity across the mammalian lifespan, and the DG is one of only two known areas in which neurogenesis occurs in adult mammals (Cameron and Mckay, 2001;Gould, 2007). We identify numerous potential candidate regulators of neurogenesis, neurotransmission and structural remodeling during and after pregnancy, with the most dramatic changes occurring between virgin (prepregnancy) and E18. These results are broadly consistent with previous in vitro and in vivo experiments demonstrating a marked impact of sex hormones on neural plasticity. Specifically, in mice, late pregnancy is associated with reduced neurogenesis in the DG that continues into the early postpartum period (Rolls et al., 2008;Kim et al., 2010). Pregnancy also correlates with increased dendritic spine density in the CA1 neurons of female rats (Kinsley et al., 2006;Pawluski and Galea, 2006), as well as changes in GABAergic neurotransmission in mice during pregnancy (Maguire et al., 2009).
Periods of intense neuroplasticity are often associated with vulnerability for mental illness, as perturbations during so-called "critical periods" can have profound longterm consequences for cognition, memory, and behavior. Therefore, it is not surprising that the peripartum period is also marked by increased incidence of psychiatric illness, including depression, anxiety, and psychosis (Galea et al., 2014;Barba-Müller et al., 2019). Future work should characterize how stress or other perturbations during pregnancy modulate peripartum neuroplasticity, as well as cell type-specific regulators of maternal behavior. For example, in pregnant rats, stress during gestation results in increased neurogenesis and also impacts hippocampal glucocorticoid receptor density (Pawluski et al., 2015).
There are likely other drivers of neuroplasticity during pregnancy that require deeper investigation beyond this study. In particular, immune adaptation across gestation is critical to safeguard the mother from infection and prevent an immune reaction directed at the fetus (Peterson et al., 2020). While classically implicated in maternal-fetal tolerance, T cells also play a role in adult murine neurogenesis, and pregnant mice lacking a T cell population only display a marginal decline in neurogenesis during pregnancy (Ziv et al., 2006;Rolls et al., 2008). Interestingly, the pregnancy-associated decrease in neurogenesis is reversed following T cell restoration in nude mice (Rolls et al., 2008).
The findings of the present study should be interpreted in the context of several technical and analytical limitations. Importantly, changes in the maternal brain are highly species-dependent, and therefore, the results of this study must be interpreted in the context of the murine brain. snRNA-seq suffers from low capture efficiency and gene dropout, which may bias the analysis toward more highly expressed genes. Similarly, chromatin accessibility from snATAC-seq is sparse, and it is challenging to infer precise biological function from changes in accessible regions and putative gene regulatory networks. The comparisons between conditions are pairwise, which limits any conclusion about the precise trajectory across all three peripartum timepoints. While we identified changes between groups that were statistically significant (FDR , 5%), the fold changes are often modest, at best, which may reflect a limitation of the single-cell sequencing approach. Despite these limitations, this dataset is an important resource for mapping cell type-specific adaptation to the unique context of pregnancy.