Dopamine drives persistent remodelling of the maternal brain

Dopamine is identified as a central regulator of long-lasting brain changes in mothers, driving persistent neural adaptations through a mechanism involving H3 dopaminylation. The study used transcriptomic profiling in mice and humans to show that reproductive experience remodels the dorsal hippocampal formation, and that chronic postpartum stress disrupts these adaptations by altering dopamine dynamics.

Abstract Pregnancy and postpartum experiences represent transformative physiological states that impose lasting demands on the maternal body and brain, resulting in lifelong neural adaptations 1,2,3,4,5,6. However, the precise molecular mechanisms that drive these persistent alterations remain poorly understood. Here we used brain-wide transcriptomic profiling to define the molecular landscape of neuroplasticity induced by reproductive experience, identifying the dorsal hippocampal formation dHF as a key site of transcriptional remodelling. Combining single-cell RNA sequencing with a maternal–pup separation paradigm, we additionally found that chronic postpartum stress significantly disrupts dHF adaptations by altering dopamine dynamics, leading to changes in the dopamine-dependent histone post-translational modification, H3 dopaminylation, which causally mediates downstream alterations in gene expression and behaviour. In human dorsal subiculum, a brain structure within the dHF, we uncovered conserved patterns of parity-dependent alterations in H3 dopaminylation and transcription. We further established the sufficiency of dopamine modulation in regulating these adaptations via chemogenetic suppression of dopamine release into the dHF, which recapitulated key epigenomic and behavioural features of reproductive experience in virgin female mice. In sum, our findings establish dopamine as a central regulator of parity-induced neuroadaptations in humans and mice, revealing a fundamental transcriptional mechanism by which female reproductive experience remodels the brain to sustain long-term behavioural adaptations. Similar content being viewed by others Main Matrescence—the physical, emotional, hormonal and social transition to motherhood—is a period of profound transformation that reshapes the maternal body and brain to support pregnancy, birth and offspring care. Although extensive research has focused on acute maternal brain processes that are essential for the onset and maintenance of parenting, precisely how brain adaptations are sustained beyond the postpartum period to persistently influence behaviour remains unclear. Recent human neuroimaging studies have revealed that parity—encompassing prior pregnancies carried to term across diverse modes of conception, delivery and infant feeding—induces long-lasting alterations in brain connectivity and structure, persisting for years or even decades following birth 1,2,3,4,5,6. Similarly, animal models of parity, which capture the composite exposures of pregnancy, parturition, lactation and maternal care hereafter referred to as reproductive experience RE , demonstrate persistent alterations in synaptic remodelling, cell proliferation and behavioural outcomes . However, the specific mechanisms within the extensive repertoire of maternal physiological changes that drive these long-lasting brain adaptations remain poorly understood. 7 ref-CR7 , 8 ref-CR8 , 9 ref-CR9 , 10 ref-CR10 , 11 ref-CR11 , 12 ref-CR12 , 13 /articles/s41586-026-10509-4 ref-CR13 Sustained transcriptional and behavioural adaptations To investigate the persistent effects of RE in maternal brain, we compared dams that experienced a composite of reproductive and maternal experiences breeding, pregnancy, parturition, lactation and pup interactions to age-matched virgin nulliparous NP mice Fig. 1a /articles/s41586-026-10509-4 Fig1 . As hormones return to baseline levels around 7 weeks post-conception 14, we conducted brain-wide transcriptomic profiling 4 weeks following pup weaning 49 days postpartum dpp . Following differential expression analyses of 11 brain regions selected on the basis of prior evidence supporting their involvement in maternal behaviours Extended Data Fig. 1a /articles/s41586-026-10509-4 Fig6 , our data revealed a wide range in the number of differentially expressed genes DEGs; adjusted P < 0.05 across brain regions Fig. 1b /articles/s41586-026-10509-4 Fig1 and Supplementary Table 1 /articles/s41586-026-10509-4 MOESM3 . To assess whether these differential expression patterns could be attributed to common changes in cell types, we performed cell-type deconvolution analyses. We observed changes in cell marker expression across brain regions, regardless of the number of DEGs Fig. 1c /articles/s41586-026-10509-4 Fig1 . Notably, the expression of neuronal markers was downregulated in the brain regions displaying the most DEGs Fig. 1c,d /articles/s41586-026-10509-4 Fig1 and Extended Data Fig. 1b–f /articles/s41586-026-10509-4 Fig6 , implicating changes in neuronal populations that may underlie observed transcriptional responses. As DEG differences across brain regions are influenced by cell-type composition, we repeated our analyses with surrogate cell-type proportion variables included in the model. Incorporating these covariates altered the relative ordering of brain regions by DEG number Extended Data Fig. 1g /articles/s41586-026-10509-4 Fig6 and Supplementary Table 2 /articles/s41586-026-10509-4 MOESM3 . Remarkably, however, the dHF remained the most transcriptionally sensitive region, demonstrating that RE induces robust cell-intrinsic transcriptional remodelling in this region. Consistent with this, pathway enrichment analyses comparing DEGs from the original and composition-adjusted dHF analyses revealed strong concordance in the biological processes identified, indicating that both approaches robustly capture core transcriptional programs engaged by RE Extended Data Fig. 1h,i /articles/s41586-026-10509-4 Fig6 . Notably, variation in litter size greater than zero did not meaningfully alter gene expression patterns within RE Extended Data Fig. 1j /articles/s41586-026-10509-4 Fig6 . We next explored whether the transcriptomic alterations in RE-sensitive regions might be driven by common upstream regulators. We returned to our unadjusted model, which preserves the full extent of RE-dependent transcriptional changes, to compare DEGs across brain regions, enabling the identification of common gene sets altered by RE. This revealed that the greatest overlap occurs in regions displaying the highest levels of transcriptional sensitivity Fig. 1e /articles/s41586-026-10509-4 Fig1 . Next, we grouped the 11 brain regions according to the number of DEGs, the degree of overlap with other regions, and significant fold changes in neuronal marker expression, based on the premise that these regions may be co-regulated to elicit convergent transcriptional responses referred to as high-sensitivity regions: dHF, nucleus accumbens NAc , medial preoptic area mPOA , locus coeruleus LC , ventral hippocampus vHF . Brain regions lacking these criteria were classified as low-sensitivity regions dorsal raphe nucleus DRN , paraventricular nucleus of the hypothalamus PVN , basolateral amygdala BLA , ventral tegmental area VTA , olfactory bulb OB and the medial prefrontal cortex mPFC . Next, we explored the predicted upstream regulators of overlapping DEGs in high- versus low-sensitivity brain regions, which revealed greater enrichment for oestrogen, progesterone, testosterone, dopamine and other ligands with similar receptor affinity or structural homology in high-sensitivity regions Fig. 1f /articles/s41586-026-10509-4 Fig1 . We next performed weighted gene co-expression network analysis WGCNA across all 11 brain regions, resulting in 9 co-expression modules Extended Data Fig. 1k–n /articles/s41586-026-10509-4 Fig6 . Module–trait correlations identified significant correlations between high-sensitivity regions in RE dams and the brown, green, yellow and magenta modules hereafter termed RE-sensitive modules; Fig. 1g /articles/s41586-026-10509-4 Fig1 . Genes within RE-sensitive modules displayed significant enrichment for DEGs from dHF, NAc, mPOA, LC and vHF Fig. 1h /articles/s41586-026-10509-4 Fig1 . Functional annotation analysis of the genes from these modules highlighted signalling pathways related to neuromodulators, such as oestrogen yellow and brown and dopamine green and magenta , in agreement with our upstream regulator analyses Fig. 1i /articles/s41586-026-10509-4 Fig1 and Extended Data Fig. 1o,p /articles/s41586-026-10509-4 Fig6 . We next investigated whether such transcriptomic findings correspond to sustained functional alterations in RE females. We performed behavioural assessments about 16 weeks after pup weaning based on functions associated with brain regions displaying the most robust transcriptional changes, including maternal behaviour, learning and memory. To assess maternal behaviour, we tested pup retrieval in the home cage by placing two donor pups in opposite corners away from the nest. RE dams retrieved pups with significantly reduced latency compared with NP females Fig. 1j /articles/s41586-026-10509-4 Fig1 , consistent with previous findings, although at a more extended timepoint in our studies 15. Next, in a contextual fear conditioning paradigm, RE females exhibited both enhanced acquisition and context recall compared with NP Fig. 1k,l /articles/s41586-026-10509-4 Fig1 , supporting previous reports of cognitive adaptations . No significant differences were observed in the open field, light–dark box or forced swim tests Extended Data Fig. 7 /articles/s41586-026-10509-4 ref-CR7 , 10 /articles/s41586-026-10509-4 ref-CR10 , 12 /articles/s41586-026-10509-4 ref-CR12 , 13 /articles/s41586-026-10509-4 ref-CR13 , 16 ref-CR16 , 17 ref-CR17 , 18 ref-CR18 , 19 ref-CR19 , 20 ref-CR20 , 21 ref-CR21 , 22 /articles/s41586-026-10509-4 ref-CR22 2a–h /articles/s41586-026-10509-4 Fig7 . Additionally, the distribution of oestrous stages was comparable between groups Extended Data Fig. 2i–n /articles/s41586-026-10509-4 Fig7 . Contributions of pregnancy and postpartum events We next sought to resolve the reproductive event s that contributed most substantially to the sustained transcriptional alterations that we observed. To do so, we compared females that were successfully bred to males but did not become pregnant mating + no pregnancy , females that experienced pregnancy and parturition but did not transition to postpartum owing to pup removal on the day of birth mating + pregnancy + birth , and virgin females exposed to 21 days of donor pup interactions pup sensitized; Extended Data Fig. 3a,b /articles/s41586-026-10509-4 Fig8 . We focused our investigations on the dHF owing to its pronounced transcriptional sensitivity to RE. When comparing gene expression profiles in NP versus RE, we observed that the mating + pregnancy + birth group most closely resembled changes observed with RE, suggesting that pregnancy, and potentially its combination with birth, is a primary driver of dHF changes Extended Data Fig. 3c,d /articles/s41586-026-10509-4 Fig8 . However, the magnitude of fold changes did not match that of the RE group, indicating that additional experiences are necessary. Although visualization of differential expression profiles in the mating + no pregnancy and pup sensitized groups did not mimic changes seen with RE, significant DEGs across all groups compared with NP overlapped with those in RE Extended Data Fig. 3e–h /articles/s41586-026-10509-4 Fig8 and Supplementary Table 3 /articles/s41586-026-10509-4 MOESM3 . Furthermore, DEGs from each comparison significantly overlapped with genes from RE-sensitive modules, with the strongest enrichment in mating + pregnancy + birth and pup sensitized groups Extended Data Fig. 3i /articles/s41586-026-10509-4 Fig8 . Thus, despite distinct expression profiles between pup sensitized versus RE, pup interactions are likely to converge on similar regulatory networks, albeit through different DEGs that are critical for RE-driven transcriptional alterations. To further explore these relationships, we performed threshold-free rank–rank hypergeometric overlap RRHO analyses, which demonstrated transcriptional concordance across experiences, with the strongest overlap observed in the pregnancy-exposed group Extended Data Fig. 3j /articles/s41586-026-10509-4 Fig8 . Furthermore, comparing gene ontology GO term analyses of DEGs for each group revealed significant enrichment in overlapping biological processes Extended Data Fig. 3k /articles/s41586-026-10509-4 Fig8 . Notably, certain pathways were unique to the combination of reproductive exposures, including changes in CREB activity, catecholamine secretion and myelination. We next conducted a controlled time-course study to examine the transcriptional dynamics across early and late gestational and postpartum stages in age-matched females Extended Data Fig. 4a /articles/s41586-026-10509-4 Fig9 . Based on the hypothesis that RE-induced transcriptional alterations require multiple reproductive events to achieve full induction, we first compared DEGs at each timepoint versus NP to genes from RE-sensitive modules. DEGs from all comparisons significantly overlapped with the brown module, with postpartum exposures also enriching for the green and yellow modules Extended Data Fig. 4b /articles/s41586-026-10509-4 Fig9 , left . To assess the regional specificity of these dynamic changes, we analysed bulk RNA-sequencing RNA-seq data from the vHF, a region that displayed moderate sensitivity to RE Extended Data Fig. 4b /articles/s41586-026-10509-4 Fig9 , middle , along with the mPFC, which exhibited few DEGs Extended Data Fig. 4b /articles/s41586-026-10509-4 Fig9 , right . These regions had limited overlap with RE-sensitive module genes, emphasizing the specificity of pathways in brain regions displaying heightened sensitivity to RE Supplementary Tables 4 /articles/s41586-026-10509-4 MOESM3 – 6 /articles/s41586-026-10509-4 MOESM3 . To characterize cumulative transcriptional changes, we performed temporal analyses to detect genes displaying sustained or transient expression alterations. This revealed a prominent set of genes showing persistent and progressive downregulation, most pronounced during postpartum Extended Data Fig. 4c /articles/s41586-026-10509-4 Fig9 , that enriched for processes related to glutamatergic synapses, dopaminergic signalling and endocrine-related pathways Extended Data Fig. 4d /articles/s41586-026-10509-4 Fig9 . Conversely, genes with persistent upregulation enriched for pathways associated with cell junction dynamics Extended Data Fig. 4e /articles/s41586-026-10509-4 Fig9 . Additionally, we identified gene sets displaying transient alterations, including downregulation of neurogenesis-related pathways and upregulation of copper ion metabolism Extended Data Fig. 4f,g /articles/s41586-026-10509-4 Fig9 , consistent with other studies reviewed in refs. 23,24 . In total, these findings highlight postpartum as a key window for reinforcing RE-induced alterations. Chronic postpartum stress disrupts dHF adaptations We next examined whether postpartum perturbations may disrupt pathways that reinforce RE-induced alterations in dHF. We implemented a maternal stress paradigm that increases stress hormones and disrupts key postpartum experiences, including nursing and pup interactions, from 10–20 dpp Fig. 2a /articles/s41586-026-10509-4 Fig2 and Extended Data Fig. 5a /articles/s41586-026-10509-4 Fig10 . During this period, dams were separated from their pups for 3 h daily and received limited nesting until pup weaning stress RE . No differences in litter sizes or postpartum weights were observed Extended Data Fig. 5b,c /articles/s41586-026-10509-4 Fig10 . Comparison of stress RE dHF transcriptional profiles revealed intermediate clustering between NP and control RE groups Fig. 2b /articles/s41586-026-10509-4 Fig2 , Extended Data Fig. 5d /articles/s41586-026-10509-4 Fig10 and Supplementary Table 7 /articles/s41586-026-10509-4 MOESM3 , suggesting that postpartum stress disrupts the extent of RE alterations in dHF. Approximately 85% of control RE versus stress RE DEGs showed opposing directionality from NP versus control RE, enriching across pathways associated with long-term potentiation, dopaminergic synapses, oxytocin signalling and other processes Fig. 2c,d /articles/s41586-026-10509-4 Fig2 . Since our pathway analyses indicated changes in long-term potentiation, we focused our behavioural assessments on contextual fear conditioning. Control RE females exhibited enhanced acquisition and context recall compared to NP, whereas stress RE females showed no significant differences versus NP Fig. 2e,f /articles/s41586-026-10509-4 Fig2 . To further evaluate dHF-dependent function, we used the object location task. Given prior findings that adult females require 10 min of training for reliable discrimination 25, we first used this duration to confirm that all groups could successfully identify the moved object Extended Data Fig. 5e /articles/s41586-026-10509-4 Fig10 . We then reduced the training time to 5 min. Under these conditions, control RE mice discriminated the moved object, whereas NP and stress RE females did not Extended Data Fig. 5f /articles/s41586-026-10509-4 Fig10 , an effect that was not attributed to oestrous stage or locomotion Extended Data Fig. 5g,h /articles/s41586-026-10509-4 Fig10 . We next tested whether the same maternal stress paradigm from 2–12 dpp may also disrupt dHF adaptations Extended Data Fig. 6a /articles/s41586-026-10509-4 Fig11 . Similar to the ‘late’ stress RE group, dHF transcriptome profiles from the early stress RE group clustered between those of NP and control RE groups Extended Data Fig. 6b,c /articles/s41586-026-10509-4 Fig11 and Supplementary Table 8 /articles/s41586-026-10509-4 MOESM3 . DEGs with opposing directionality between NP versus control RE and control RE versus early stress RE had about 41% overlap Extended Data Fig. 6d,e /articles/s41586-026-10509-4 Fig11 , and associated with synaptic plasticity pathways Extended Data Fig. 6f /articles/s41586-026-10509-4 Fig11 . Thus, early and late postpartum stress similarly perturb dHF RE transcriptomic adaptations, with early stress exerting a less extensive effect. Consistent with this finding, there were no differences between early stress RE versus NP in contextual fear conditioning or the object location task Extended Data Fig. 6g–i /articles/s41586-026-10509-4 Fig11 . Dopaminergic modulation underlying dHF plasticity To identify the mechanisms underlying RE-induced dHF plasticity, we next leveraged the transcriptomic and behavioural signatures shared between NP and stress RE groups, focusing all subsequent analyses on the late postpartum stress paradigm given stronger disruption to RE adaptations. To identify the dHF cell-types that integrate hormone and neurotransmitter signalling to drive RE-related plasticity, we next performed single-nuclei RNA sequencing snRNA-seq . Following rigorous quality control assessments to remove confounding sources of variation, such as mitochondrial mapping percentage, doublets, and cell-type contamination Extended Data Fig. 7a–f /articles/s41586-026-10509-4 Fig12 , we retained 109,334 nuclei for downstream analysis NP n = 37,070 ; control RE n = 35,631 ; stress RE n = 36,633 . Cell-type annotation resulted in 16 distinct clusters Fig. 2g /articles/s41586-026-10509-4 Fig2 and Extended Data Fig. 7g–i /articles/s41586-026-10509-4 Fig12 , including excitatory neuron subtypes dentate gyrus DG , CA1, CA3, subiculum and mossy cells; 43.3% of total , γ-aminobutyric acid-producing GABAergic neurons GABA.1–GABA.4; 9.1% of total , neural progenitors neural stem cells and radial glia-like cells; 9.7% of total and non-neuronal or glial populations astrocytes, oligodendrocytes, OPCs, microglia and immune cells; 37.9% of total . Since our cell-type deconvolution analyses indicated downregulation of neuronal markers, we first confirmed shifts in proportional differences across major cell types, with GABA.2 and GABA.4 clusters reduced in control RE dHF compared to both NP and stress RE Extended Data Fig. 7j,k /articles/s41586-026-10509-4 Fig12 . Following pseudobulk aggregation, we identified DEGs across both neuronal and non-neuronal subtypes, with excitatory neuron alterations predominantly observed between NP versus control RE, particularly within the subiculum Extended Data Fig. 7l /articles/s41586-026-10509-4 Fig12 . By contrast, the most pronounced neuronal changes between control and stress RE were observed in the GABA.2 and mossy cell populations Extended Data Fig. 7l /articles/s41586-026-10509-4 Fig12 . GO term analysis of DEGs from cell-types displaying proportional differences enriched for pathways associated with responses to endogenous stimuli, cell proliferation, synaptic function and metabolic processes Extended Data Fig. 7m,n /articles/s41586-026-10509-4 Fig12 . Notably, DEGs shared between single-cell and bulk analyses enriched for similar pathways, indicating cross-platform convergence Extended Data Fig. 7o–q /articles/s41586-026-10509-4 Fig12 . We next focused on the GABA.2 cluster owing to shared proportional and transcriptional changes observed in comparison to both NP and stress RE groups, which suggested a common mechanism limiting dHF plasticity. Marker gene analysis revealed significant enrichment of dopamine receptor D1 Drd1 and D2 Drd2 mRNA expression Extended Data Fig. 7r /articles/s41586-026-10509-4 Fig12 , aligning with our prior bioinformatic analyses implicating dopaminergic regulation. To further resolve this population, we subclustered the GABA.2 cells and identified distinct subpopulations expressing Drd1 and Drd2 , both of which appeared to be altered in control RE mice Fig. 2h /articles/s41586-026-10509-4 Fig2 . In a separate cohort, we validated these changes using RNA fluorescence in situ hybridization FISH in dorsal CA1 Fig. 2i /articles/s41586-026-10509-4 Fig2 and Extended Data Fig. 7s,t /articles/s41586-026-10509-4 Fig12 , where dopamine receptor-expressing interneurons have been described 26,27. We also examined whether such changes extend to other dopamine-sensitive neuronal populations. RNA FISH in the DG, which contains excitatory Drd1 and Drd2 -expressing neurons in granule and hilar mossy cells, respectively , similarly revealed significant changes in the distribution of Drd1 and Drd2 puncta across groups, alongside expected regional variations consistent with the sparser distribution of dopamine receptor-expressing neurons in CA1 Fig. 28 /articles/s41586-026-10509-4 ref-CR28 2j /articles/s41586-026-10509-4 Fig2 and Extended Data Fig. 7u /articles/s41586-026-10509-4 Fig12 . Since these findings established dopaminergic signalling as a potential mediator of persistent RE-induced adaptations, this prompted us to examine whether maternal stress disrupts this process within its defined postpartum window. Given prior evidence that pup separation elevates dopamine levels in NAc 29,30, we first investigated whether a similar response occurs in dHF Fig. 2k /articles/s41586-026-10509-4 Fig2 . Dopamine levels increased during pup separation, with expected regional differences in NAc versus dHF reflecting the strength of innervation to these regions Fig. 2l /articles/s41586-026-10509-4 Fig2 . Since acute and chronic stress produce different biological responses, we next assessed the effect of repeated separation stress, revealing elevated baseline dopamine in both dHF and NAc tissues Fig. 2m /articles/s41586-026-10509-4 Fig2 and Extended Data Fig. 7v /articles/s41586-026-10509-4 Fig12 . Thus, decreased Drd1/2 -expressing neurons may reflect a compensatory response to lower basal dopamine, consistent with our findings of reduced dopamine-related transcriptional activity in control RE. Parity alters H3 dopaminylation in mice and humans Together, changes in dopamine modulation and gene expression point to a role for dopamine-dependent epigenetic mechanisms in shaping persistent transcriptional states. Motivated by this idea, we turned to a recently characterized class of histone post-translational modifications that is dependent on intracellular pools of biogenic monoamines, including serotonin, dopamine and histamine, termed monoaminylations 31,32,33. This process involves transamidation of monoamines onto the glutamine 5 residue of histone H3 H3Q5 , a site adjacent to the canonically permissive lysine 4 tri-methylation H3K4me3 , by tissue transglutaminase 2 TG2; Extended Data Fig. 34 /articles/s41586-026-10509-4 ref-CR34 8a–c /articles/s41586-026-10509-4 Fig13 . Monoaminylation at H3Q5 can coexist with H3K4me3 and is responsive to environmental stimuli such as stress or drugs of abuse, with downstream effects on transcriptional output . Using our previously validated H3K4me3Q5dop antibody 31 ref-CR31 , 32 ref-CR32 , 33 /articles/s41586-026-10509-4 ref-CR33 , 35 ref-CR35 , 36 ref-CR36 , 37 ref-CR37 , 38 ref-CR38 , 39 ref-CR39 , 40 /articles/s41586-026-10509-4 ref-CR40 , we performed CUT&RUN-seq to test the involvement of H3 dopaminylation in RE transcriptional adaptations. The majority of H3K4me3Q5dop peaks annotated to genic loci, with a large proportion occurring within 2 kb of transcriptional start sites TSSs; Extended Data Fig. 32 /articles/s41586-026-10509-4 ref-CR32 8d–g /articles/s41586-026-10509-4 Fig13 . Across all groups, higher H3K4me3Q5dop enrichment associated with increased gene expression Extended Data Fig. 8h,i /articles/s41586-026-10509-4 Fig13 . We next performed differential peak analysis and observed that the majority of group-specific changes in H3 dopaminylation were downregulated in control RE dHF compared with both NP and stress RE Extended Data Fig. 8j /articles/s41586-026-10509-4 Fig13 and Supplementary Tables 9 /articles/s41586-026-10509-4 MOESM3 and 10 /articles/s41586-026-10509-4 MOESM3 . Integrating these peaks with DEGs from the same comparisons revealed significant overlap between downregulated peaks and genes decreased in control RE Extended Data Fig. 8k /articles/s41586-026-10509-4 Fig13 . H3K4me3Q5dop enrichment at the TSS of RE-regulated DEGs was also reduced at most loci compared to NP and stress RE Fig. 3a,b /articles/s41586-026-10509-4 Fig3 and Extended Data Fig. 8l /articles/s41586-026-10509-4 Fig13 . Stratifying DEGs by directionality revealed enrichment of canonical repressors among downregulated genes and ligand-responsive nuclear receptors among upregulated genes, pointing to ligand-responsive transcriptional regulators that may engage additional chromatin remodelling pathways Fig. 3c,d /articles/s41586-026-10509-4 Fig3 . To assess the translational relevance of these findings, we next performed transcriptomic and H3K4me3Q5dop profiling in human dorsal subiculum dSub , a subregion of the dHF that serves as its primary output and was identified in our snRNA-seq analysis as the most transcriptionally responsive subregion to RE. Comparing age-matched NP individuals with parous individuals who had 1–2 previous births, we identified robust differential gene expression indicative of long-lasting effects of parity in human brain Fig. 3e /articles/s41586-026-10509-4 Fig3 and Supplementary Tables 11 /articles/s41586-026-10509-4 MOESM3 and 12 /articles/s41586-026-10509-4 MOESM3 . These transcriptional signatures showed strong concordance with mouse pseudobulk expression profiles Fig. 3f /articles/s41586-026-10509-4 Fig3 , and were enriched for pathways related to extracellular matrix remodelling and metabolism Extended Data Fig. 8m /articles/s41586-026-10509-4 Fig13 . Furthermore, both upregulated and downregulated DEGs in human dSub were enriched for upstream regulators similarly identified in mouse Fig. 3g /articles/s41586-026-10509-4 Fig3 . We also performed H3K4me3Q5dop CUT&RUN-seq Extended Data Fig. 8n–r /articles/s41586-026-10509-4 Fig13 in human brain, where differential binding analysis showed widespread reduction in dopaminylation signal in parous dHF, corresponding with DEG analyses Fig. 3h,i /articles/s41586-026-10509-4 Fig3 , Extended Data Fig. 8s /articles/s41586-026-10509-4 Fig13 and Supplementary Table 13 /articles/s41586-026-10509-4 MOESM3 . Pathway analysis of differential peaks revealed significant enrichment for loci involved in maternal behaviour, behavioural fear responses, hormone regulation and dopamine transport Fig. 3j /articles/s41586-026-10509-4 Fig3 , suggesting conserved, RE-associated remodelling of dopamine-sensitive transcriptional programs in dHF. Dopamine suppression recapitulates RE outcomes We next hypothesized that persistent reductions in dopamine tone may engage adaptive chromatin mechanisms to support enduring transcriptional and behavioural plasticity. As the VTA showed similar response patterns to pup separation Extended Data Fig. 7w /articles/s41586-026-10509-4 Fig12 , we chemogenetically suppressed the dopaminergic VTA projection to dHF by bilaterally expressing a floxed inhibitory neuronal hM4Di-DREADD fused to mCherry versus floxed mCherry controls in the VTA, coupled with retrograde adeno-associated viral vectors AAVs expressing a tyrosine hydroxylase TH -dependent Cre recombinase into the dHF Fig. 4a /articles/s41586-026-10509-4 Fig4 . We validated this approach by confirming selective mCherry expression in TH+ neurons of the VTA, as well as reduced FOS expression following administration of the selective DREADD agonist deschloroclozapine DCZ; Fig. 4b /articles/s41586-026-10509-4 Fig4 and Extended Data Fig. 9a,b /articles/s41586-026-10509-4 Fig14 . To minimize injection-related stress, we restricted the window of chronic DCZ administrations. Since maternal stress from 10–20 dpp was sufficient to disrupt RE-induced plasticity, we posited that dHF dopamine suppression during this period may be sufficient to phenocopy RE-induced adaptations. Therefore, we administered DCZ to virgin NP females versus postpartum dams in parallel NP-mCherry, NP-hM4Di, RE-mCherry and RE-hM4Di; Fig. 4c /articles/s41586-026-10509-4 Fig4 . After confirming that our approach suppressed dHF dopamine Fig. 4d /articles/s41586-026-10509-4 Fig4 and did not induce gross metabolic changes Extended Data Fig. 9c /articles/s41586-026-10509-4 Fig14 , we assessed behavioural outcomes one month after the final DCZ injection. In the pup retrieval test, failure to retrieve either pup occurred in 8 out of 13 NP-mCherry females, a result that occurred in only 1 out of 15 NP-hM4Di, 0 out of 12 RE-mCherry, and 1 out of 11 RE-hM4Di females Fig. 4e /articles/s41586-026-10509-4 Fig4 . In contextual fear conditioning, NP-hM4Di females demonstrated enhanced acquisition, with no significant differences in the context recall test Fig. 4f /articles/s41586-026-10509-4 Fig4 and Extended Data Fig. 9d /articles/s41586-026-10509-4 Fig14 . Oestrous stage analyses confirmed that the observed behavioural outcomes were driven by RE, independent of hormonal state or viral expression Extended Data Fig. 9e,f /articles/s41586-026-10509-4 Fig14 . Next, we examined whether dopamine suppression influences cellular receptor expression in the dHF. Similar to control RE, NP-hM4Di females had greater distribution of cells with low expression of Drd1 and Drd2 , compared with NP-mCherry Fig. 4g,h /articles/s41586-026-10509-4 Fig4 and Extended Data Fig. 9g–i /articles/s41586-026-10509-4 Fig14 . Following transcriptional profiling to assess the degree to which dopamine suppression mirrors RE-associated gene expression patterns, we observed that NP-hM4Di exhibited an intermediate expression pattern relative to NP-mCherry versus RE-mCherry expression profiles Extended Data Fig. 9j /articles/s41586-026-10509-4 Fig14 and Supplementary Table 14 /articles/s41586-026-10509-4 MOESM3 . RRHO analysis further supported transcriptional concordance between RE-mCherry and NP-hM4Di groups compared with NP-mCherry Fig. 4i /articles/s41586-026-10509-4 Fig4 . Of note, NP-mCherry versus NP-hM4Di DEGs overlapped with genes from the RE-sensitive brown module Extended Data Fig. 9k /articles/s41586-026-10509-4 Fig14 . Functional annotation of these DEGs revealed enrichment in pathways governing synaptic plasticity and responses to endogenous stimuli Fig. 4j /articles/s41586-026-10509-4 Fig4 , similarly mirroring patterns altered in control RE dHF. Finally, to assess whether the effects of dopamine suppression on dHF transcription are modulated in part through H3 dopaminylation, we conducted H3K4me3Q5dop CUT&RUN-seq in virally infected NP tissues. Examining H3K4me3Q5dop enrichment at the TSS of both RE-regulated DEGs and genes altered between NP-mCherry and NP-hM4Di groups, we observed significantly reduced signal at all loci examined in NP-hM4Di tissues Fig. 4k /articles/s41586-026-10509-4 Fig4 , Extended Data Fig. 9l,m /articles/s41586-026-10509-4 Fig14 and Supplementary Table 15 /articles/s41586-026-10509-4 MOESM3 . We identified significant transcription factors regulating DEGs between NP-mCherry and NP-hM4Di, overlapping with those observed in parous individuals in both human dSub and our mouse model Fig. 4l /articles/s41586-026-10509-4 Fig4 . Reversing H3 dopaminylation rescues stress effects Finally, to determine whether H3 dopaminylation causally regulates RE adaptations, we manipulated H3 dopaminylation in dHF by transducing AAVs expressing histone H3.3 constructs—which incorporate into chromatin in post-mitotic neurons—encoding either the wild-type glutamine at position 5 H3.3 WT, control or a glutamine-to-alanine substitution H3.3 Q5A , which functions as a dominant negative mutant to reduce monoaminylated H3Q5 31,32,33,35,39. We targeted H3 dopaminylation in stress RE dHF, as this group showed higher H3 dopaminylation levels versus control RE, and a more pronounced increase versus NP. Accordingly, stress RE dams received either H3.3 WT or H3.3 Q5A AAVs, and were compared to control RE dams receiving H3.3 WT as a control Fig. 5a,b /articles/s41586-026-10509-4 Fig5 . To resolve the precise genomic loci affected by H3.3 Q5A expression, we performed H3K4me3Q5dop CUT&RUN-seq on infected dHF tissues. Remarkably, differential H3 dopaminylation sites that distinguished control RE H3.3 WT from stress RE H3.3 WT were reversed in stress RE H3.3 Q5A Fig. 5c /articles/s41586-026-10509-4 Fig5 . As expected from prior studies validating this approach 31,32,33,35,39, the majority of dopaminylation peaks altered in stress RE H3.3 Q5A were downregulated compared with stress RE H3.3 WT Fig. 5d /articles/s41586-026-10509-4 Fig5 , Extended Data Fig. 10a /articles/s41586-026-10509-4 Fig15 and Supplementary Tables 16 /articles/s41586-026-10509-4 MOESM3 and 17 /articles/s41586-026-10509-4 MOESM3 . These significantly reversed loci enriched for pathways involved in dopamine transport, glutamatergic signalling and learning Fig. 5e /articles/s41586-026-10509-4 Fig5 , suggesting that H3 dopaminylation contributes to cognition-related behaviours. Similar to our chromatin data, bulk transcriptomic analyses in dHF showed that H3.3 Q5A transduction reversed gene expression profiles between control RE H3.3 WT and stress RE H3.3 WT Fig. 5f /articles/s41586-026-10509-4 Fig5 and Supplementary Table 18 /articles/s41586-026-10509-4 MOESM3 . DEGs from both comparisons overlapped with the RE-sensitive green modules, indicating that H3 dopaminylation modulates gene networks central to RE-associated plasticity Extended Data Fig. 10b /articles/s41586-026-10509-4 Fig15 . Indeed, hierarchical clustering of all DEGs revealed that H3.3 Q5A expression in stress RE dHF partially restores transcriptional patterns observed in control RE Fig. 5g /articles/s41586-026-10509-4 Fig5 . To evaluate whether changes in H3 dopaminylation correspond to transcriptional alterations, we next intersected differential H3K4me3Q5dop binding sites with genes overlapping between control RE H3.3 WT versus stress RE H3.3 WT and stress RE H3.3 WT versus stress RE H3.3 Q5A comparisons. Functional annotation analysis of these overlapping genes showed significant enrichment for learning-related pathways, including synaptic plasticity and glutamate signalling processes Extended Data Fig. 10c /articles/s41586-026-10509-4 Fig15 . Because transcriptomic profiles at 49 dpp reflect both primary and secondary consequences of reducing H3 dopaminylation, we next examined upstream regulators that may contribute to the broader gene expression changes induced by H3.3 Q5A transduction. We identified transcription factors predicted to regulate DEGs between stress RE H3.3 WT and stress RE H3.3 Q5A , including several shared with parous human dSub and our mouse model Fig. 5h /articles/s41586-026-10509-4 Fig5 and Extended Data Fig. 10d /articles/s41586-026-10509-4 Fig15 , suggesting that H3 dopaminylation shapes RE-dependent plasticity by modulating not only target genes, but also transcription factors and chromatin remodellers that coordinate downstream expression. Finally, we assessed whether reversing H3 dopaminylation in stress RE dams could influence contextual fear conditioning outcomes. Both control RE H3.3 WT and stress RE H3.3 Q5A females showed enhanced acquisition compared to stress RE H3.3 WT Fig. 5i /articles/s41586-026-10509-4 Fig5 . During the context test, all groups demonstrated recall, with control RE H3.3 WT and stress RE H3.3 Q5A females showing higher freezing responses than stress RE H3.3 WT females Fig. 5j /articles/s41586-026-10509-4 Fig5 . Together, these data directly link RE-dependent dHF transcriptional plasticity with persistent behavioural adaptations, and establish H3 dopaminylation as a key regulator of transcriptional networks impacted by RE and postpartum stress Extended Data Fig. 10e /articles/s41586-026-10509-4 Fig15 . Discussion Through controlled brain-wide, time-course and cell-type-specific transcriptomic analyses—paired with robust behavioural outputs—we identified gene networks, reproductive events and neuromodulatory-dependent epigenetic processes that dynamically shape regional sensitivity across pregnancy and postpartum experiences, thereby promoting plasticity long after these stages have ended. In particular, our study identified the dHF as a key site of sustained plasticity, consistent with its established roles in spatial cognition, novelty detection and sensory integration, which may bestow parous dams with heightened sensitivity to salient environmental cues relevant for resource foraging, nest navigation and offspring survival 41,42,43,44. Consistent with our findings, lesion studies corroborate hippocampal involvement in pup retrieval behaviour , potentially through downstream regulation of the NAc and ventral pallidum—key regulators of pup retrieval 45 /articles/s41586-026-10509-4 ref-CR45 , 46 /articles/s41586-026-10509-4 ref-CR46 . Notably, human studies show that maternal early life stress alters hippocampal responses to infant cues 44 /articles/s41586-026-10509-4 ref-CR44 , 47 /articles/s41586-026-10509-4 ref-CR47 , supporting chronic stress as a common interferer in the normal trajectory of maternal dHF neuroadaptations. Building on this, our maternal–pup separation stress model revealed that dopamine dynamics are essential for maintaining maternal adaptations in the dHF. Indeed, we established the necessity and sufficiency of maintaining appropriate dopamine tone for dHF plasticity, demonstrating that chronic dopamine elevation through maternal stress and chemogenetic suppression within virgin dHF bidirectionally modulate RE-related changes at transcriptional, epigenetic, behavioural and cellular levels. Moreover, we establish conserved transcriptional and epigenetic changes in the human dSub, a subregion of the dHF. In our postpartum stress mouse model, reversing aberrant changes in H3 dopaminylation was sufficient to restore RE adaptations. These findings implicate dopamine as a conserved driver of maternal brain remodelling, linking its transient roles during pregnancy and postpartum to enduring neuroplasticity in critical brain regions, including the dHF. 48 /articles/s41586-026-10509-4 ref-CR48 , 49 /articles/s41586-026-10509-4 ref-CR49 Although our study focused largely on dopamine, it is notable that our chemogenetic and viral manipulations recapitulated key transcriptional and behavioural outcomes of RE, albeit not to the same extent. This underscores dopamine as a crucial but non-exclusive driver of RE-induced adaptations. Indeed, our data also implicate hormonal contributions, including oestrogen, progesterone and oxytocin signalling, in shaping transcriptional programming. Consistent with this idea, previous studies have linked parity to oestradiol-driven neuronal alterations 12, whereas others have demonstrated that blocking oxytocin signalling disrupts RE-induced spatial learning enhancements , supporting an integration of multiple neuromodulatory pathways. However, given that oestrogen signalling during postpartum regulates oxytocin release, which in turn modulates VTA signalling, we propose that dopamine adaptations serve as a key downstream mechanism within a broader network driving sustained maternal neuroplasticity. Moreover, although we aimed to parse the contributions of individual events to RE-dependent dHF programming, fully disentangling these factors remains experimentally challenging. Mating that results in pregnancy may exert distinct effects from mating that does not result in pregnancy, and isolating pregnancy, birth or lactation would require invasive procedures that introduce confounds for brain outcomes. Nevertheless, our data indicate that pregnancy induces substantial dHF transcriptional changes, although they are not characteristic of the full composite of RE, emphasizing the additional contribution of healthy postpartum exposures to dHF programming. This raises the possibility that restoring elements of healthy postpartum exposure, through pharmacological interventions or pup-associated sensory cues, may mitigate disruptions imposed by postpartum stress. More broadly, our work suggests that RE drives a fundamental shift in brain state that may interact with additional experiences. As parity is implicated in both risk and resilience to brain disorders 21 /articles/s41586-026-10509-4 ref-CR21 , future research should consider parity status as a key variable shaping differential outcomes, particularly in its interactions with other risk factors. In particular, our findings highlight an interplay between RE and stress in shaping maternal brain outcomes, revealing epigenetic mechanisms and neuromodulatory pathways that underlie long-lasting adaptations, and underscore the importance of stress mitigation during pregnancy and postpartum. 50 /articles/s41586-026-10509-4 ref-CR50 Methods Human participants Brain tissues used in this study were provided by the Douglas Brain Bank RRID: SCR 025991 https://scicrunch.org/resolver/SCR 025991/ with ethical approval from the Research Ethics Board REB of the Centre Intégré Universitaire de Santé et de Services Sociaux CIUSSS de l’Ouest-de-l'Île-de-Montréal. Analyses were conducted in accordance with all relevant guidelines and regulations. Informed consent from next of kin was obtained for each individual included in this study. Psychological autopsies, considered the gold standard for obtaining information on deceased individuals 51,52, were conducted. In brief, these consist of a series of proxy-based, structured interviews assessing psychopathology with next of kin and complemented by reviews of medical records, as previously described . Groups were matched for depression diagnosis, and were otherwise neurotypical individuals who died suddenly without prolonged agonal periods and did not have evidence of axis I disorders. Groups were matched for postmortem interval, tissue pH and RNA Integrity number. Frozen histological grade samples of grey and white matter were dissected from the subiculum by expert neuroanatomists and stored at –80 °C. Dissections were performed on 0.5-cm-thick coronal sections with the guidance of a human brain atlas 51 /articles/s41586-026-10509-4 ref-CR51 53 /articles/s41586-026-10509-4 ref-CR53 http://www.thehumanbrain.info/brain/bn brain atlas/brain.html http://www.thehumanbrain.info/brain/bn brain atlas/brain.html . Subiculum samples were obtained from sections equivalent to plate 43 of the atlas level of lateral geniculate nucleus , by dissecting through the hippocampal fissure with a slight upward angle, up to the beginning of the CA1 region. All procedures complied with ethical regulations for research involving human tissue. Demographic information can be found in Supplementary Table 11 /articles/s41586-026-10509-4 MOESM3 . Animals Wild-type C57BL6/J mice were purchased from Jackson Laboratories at 8 weeks old, and maintained on a 12 h:12 h light:dark cycle throughout the entirety of the experiments. Mice were housed in temperature- and humidity-controlled conditions approximately 22 ± 2 °C and 50 ± 10% humidity with ad libitum access to food and water throughout the entirety of the experiments. All behavioural testing occurred during the light cycle. Experimenters were blind to experimental group, and the order of testing was counterbalanced during behavioural experiments. All animal procedures were performed in accordance with NIH guidelines and with the approval of the Institutional Animal Care and Use Committee of the Icahn School of Medicine at Mount Sinai. Breeding Adult virgin female mice were pair bred in house with age-matched males. Males were removed after a maximum of 5 days, and pregnant females were singly housed at least 2 days before parturition. Pups were counted on the day of birth 0 dpp and weaned at 21 dpp. Only dams with litters between 4–10 pups were used for all experiments. On the day of weaning, dams were group-housed into cages of 3–5 mice with animals of the same experimental condition. For timed breeding, copulation plugs were checked every morning within 1 h after lights on, where confirmation of a plug was designated as embryonic day 0.5 E0.5 , signalling the immediate removal of the female to her own cage with a nestlet. Virgin NP females were age-matched for each experimental cohort. Mating + no pregnancy females were confirmed for the presence of a copulation plug. For mating + pregnancy + birth dams, pups were removed at 0 dpp within 3 h after lights-on to minimize maternal-offspring interactions. For comparison of RE dams to pup sensitized virgin females, litters were culled to 4 pups to equate litter size with the number of pups used for each sensitized female. In all other experiments, litters were not culled. Postpartum stress paradigm Stress RE females were subjected to limited nesting and maternal separation from 10–20 dpp late stress RE or 2–12 dpp early stress RE , as previously described 54,55,56. On each day of separation, the entire litter was removed to a clean cage with Sani-Chip bedding for 3–4 h. Separations occurred during the light cycle, and the timing varied each day to minimize predictability and acclimatization. EnviroDri nesting material was depleted to one-third of the amount in control cages during the days of separation. Following pup weaning on 21 dpp, the nesting material was restored to normal levels, and dams were group-housed into cages of 3–5 with mice of the same experimental condition. Pup sensitization Pup sensitization was conducted as previously described 16,57. On the first day of pup sensitization, virgin NP females were presented with 4 pups postnatal day 5 from a cage consisting of a lactating donor dam and her litter. Donor pups were exchanged for satiated pups from the same litter every 8–12 h for 21 days. Donor pups were weighed daily to ensure continual weight-gain throughout the experiment, and only pups that exhibited consistent weight gain were used. Behavioural observations were conducted during the first four days of sensitization. Each dam was observed for 30 min during the light phase for the following maternal behaviours: licking/grooming, crouching and nestbuilding. During these observation periods, females were closely monitored to ensure that they did not display aggressive behaviour toward the donor pups. On the last day of sensitization, donor pups were weaned from the cage, and sensitized females were group-housed into cages of 3–5 mice with animals of the same experimental condition. Brain tissue collection Mice were euthanized by rapid decapitation. Whole brains were flash-frozen with cold 2-methylbutane and stored at −80 °C until further processing. Flash-frozen brains were sectioned at −20 °C using a 1 mm mouse coronal brain matrix Stoelting . Tissues enriched for the brain region of interest were micropunched using a hollow needle Ted Pella according to the Allen Brain Atlas. Brain regions were selected based on previous work supporting their involvement in maternal behaviours 58,59. Behavioural analyses All mice 18–30 weeks, depending on the time from pup weaning were handled for 2 min for two consecutive days prior to initial behavioural testing. Mice were habituated to the testing room for 1 h prior to each behavioural assay. All testing occurred during the light phase. Pup retrieval All mice were individually housed for 24 h prior to testing. Following habituation, 2 pups from a donor litter with a lactating dam aged 4–6 dpp were placed in different corners opposite to the nest in the home cage. Pup-directed interactions were recorded from above for 15 min or until both pups were successfully retrieved into the nest. Retrieval latency was calculated as time pup was first placed into the nest by mouse − time first pup was placed into the cage by experimenter . Contextual fear conditioning On day 1, mice were habituated for 10 min to the testing chamber, which consisted of a square plexiglass box with a metal grid inside a sound-attenuating cabinet wiped with 70% ethanol Med Associates . On day 2, following a 3 min baseline measurement, mice were trained with 5× 2.0 s, 0.7 mA foot shocks delivered with an intertrial interval of 90 s. Testing for conditioned fear responses freezing for a total of 5 min occurred 24 h following training. Freezing was measured using ANY-maze software v7.51 connected to a camera positioned above the testing chamber. Freezing is expressed as a percentage of the total test time or as a percentage of the 60 s prior to each shock during conditioning. Mice with freezing levels exceeding 40% prior to the first shock were excluded to remove potential confounding effects of heightened baseline responsivity that could interfere with accurately assessing conditioned fear. Open field Mice were placed in a 16 × 16 cm square arena under dim lighting for 5 min. A camera positioned overhead recorded the total distance and time spent in the centre versus periphery using Ethovision XT 11 software. Object location task For training, mice were allowed to freely explore two identical objects placed equidistant from adjacent corners of a 16 × 16 cm square arena for 5 or 10 min, before being returned to the home cage. Following a 1 h or 24 h delay, the mice were returned to the arena for testing, during which time one object was moved to the opposing corner. During the test, the mice were allowed to freely explore the objects for 5 min. A camera positioned overhead recorded time spent exploring each object using Ethovision software. The discrimination score was calculated as time spent with moved object – time spent with unmoved object / time spent with moved object + time spent with unmoved object . Light–dark box Anxiety-like behaviour was tested in an apparatus containing two interconnected 20 × 20 cm compartments Omnitech Electronics . One compartment was illuminated during the session ‘light’ side , while the other was covered by an opaque black perspex lid ‘dark’ side . The distance, time spent, and number of crossovers in each compartment were automatically recorded by Fusion software during the 10 min test. Forced swim test Mice were placed in a 4 l glass beaker with 2 l of room temperature water for 8 min. Each session was recorded and scored by a blinded observer. The total number of seconds that mice were immobile during the last 5 min of the test were recorded, as previously described 60. Observation of pup-directed behaviours in the home cage Confirmation of maternal behaviour in pup sensitized females was conducted based on prior studies 61. Mice were observed from the first to fourth day of pup sensitization, based on published work that maternal sensitization in C57BL6/J females occurs following 4 days of pup exposures . Observations occurred in the light period in the first 30 min after pups were placed in the cage. Each mouse was scored every 3 min, with the observed behaviour recorded as one or more of the following categories: nestbuilding, grooming, sniffing, crouching, nursing, eating, drinking, self-grooming or no-contact. RE dams were scored for the first 4 days following birth 0–3 dpp concurrently for comparison. 62 /articles/s41586-026-10509-4 ref-CR62 , 63 /articles/s41586-026-10509-4 ref-CR63 Oestrous cycle testing Vaginal samples were taken on each day of behavioural testing. Fifteen microlitres of sterile PBS was gently pipetted into the vagina and mounted on a glass slide. Vaginal smears were stained with crystal violet dye, washed twice with water, and cover slipped with glycerol. Three 20× images per sample were acquired on a light microscope. Oestrous stage was determined using the pretrained network of EstrousNet 64, and confirmed afterwards by a trained experimenter based on cytology of nucleated, cornified, or leukocytic cells. P values from Fisher’s exact tests confirming that oestrous stage distributions did not differ across groups are provided in Supplementary Table 19 /articles/s41586-026-10509-4 MOESM3 . All staging information is available in the Source Data. Viral transduction Mice were anaesthetized with ketamine 100 mg kg−1 and xylazine 10 mg kg−1 intraperitoneally and positioned in a stereotaxic frame Kopf instruments . Given the widespread changes in dopamine receptor expression across dHF subregions, we did not restrict our viral manipulations to a specific subregion. For chemogenetic studies, 2 µl of retrograde AAV.rTH.PI.Cre.SV40 titre ≥ 7 × 1012 viral genomes vg per ml, Addgene 107788-AAVrg was infused bilaterally into the dHF at 0.2 μl min−1 using the following coordinates: 7° angle; anterior-posterior AP −2.2 mm, medial-lateral ML ±2.0 mm, dorsal-ventral DV −2.0 mm. One microlitre of pAAV-hSyn-DIO-mCherry titre ≥ 4 × 1012 vg ml−1, Addgene 50459-AAV2 or pAAV-hSyn-DIO-hM4D Gi -mCherry titre ≥ 5 × 1012 vg ml−1; Addgene 44362-AAV2 was bilaterally infused into the VTA using the following coordinates: 7° angle; AP −3.3 mm, ML ±0.9 mm; DV −4.6 mm. For rescue of H3 dopaminylation, pAAV2-CMV-H3.3WT-IRES-GFP or pAAV2-CMV-H3.3 Q5A -IRES-GFP viruses were generated and validated as previously described 31,32,33,35,39. Approximately 3–5 days post-weaning, mice were prepared for surgery, as described above, and 1 µl of virus was bilaterally infused into the dHF using the same coordinates at a rate of 0.1 μl min −1. Needles remained in place for 7 min following injection to minimize virus diffusion. Viral validations were conducted at least 21 days post-surgery to allow for optimal viral expression and recovery. Chemogenetic manipulation Following 7 days of recovery, surgerized female mice were randomly assigned to NP or RE groups. RE females were pair bred with naïve male mice for a maximum of 5 days, and individually housed prior to parturition. From 10–20 dpp, RE females were injected subcutaneously with DCZ Tocris 7193 —to minimize off-target effects 65—at 1 μg kg −1in 1% DMSO or vehicle saline . NP females were injected with DCZ concurrently. Following weaning on 21 dpp, RE dams were group-housed with mice of the same experimental condition. Behavioural testing occurred beginning at 28 days post-weaning 49 dpp , and brain tissues were collected following contextual fear conditioning for further analyses. RNAscope in situ hybridization and analysis Fresh frozen brains were cut into 12–16-μm-thick slices in the coronal plane with a cryostat Leica CM3050-S , mounted on charged Superfrost Plus microscope slides Thermofisher P36934 , and stored at −80 °C until processing. Sections were post-fixed with 4% paraformaldehyde PFA for 1.5 h at 4 °C, and permeabolized with hydrogen peroxide 10 min, room temperature and Protease III 30 min at room temperature . The RNAscope Multiplex Fluorescent Reagent Kit v2 ACD Bio was used according to manufacturer’s instructions to sequentially stain sections with the following probes: Drd1a Mm-Drd1a-C1, 461901 and Drd2 Mm-Drd2-C3, 406501-C3 . RNAscope probes were visualized using TSA Vivid Fluorophores Tocris at 1:750 dilution. Sections were counterstained with DAPI ACD Bio and mounted using ProLong Gold Antifade Mountant Thermo . Confocal images 3 images per mouse, 1,024 × 1,024 pixels were acquired on a Zeiss LSM 780 upright microscope using a 40× objective with Zen Black software, with 5 × 1 tiled images. Images were averaged across 8 consecutive acquisitions at a bit depth of 16 bits, with 2 z -stacks acquired per image. Subcellular quantification of individual puncta per 100 μm nucleus—identified by DAPI staining—was performed for each maximum intensity projected image using QuPath software v0.5.1 66. Regions of interest CA1 and DG were annotated and determined by superimposing images onto the Allen Brain Atlas. Linear mixed-effects models were used to assess group-, region-, and probe-dependent effects on puncta density. For each cell, puncta density puncta per 100 μm 2nuclear area was modelled with Group, Region, and Probe as fixed effects, including all interaction terms, and mouse identity was included as a random intercept to account for within-animal clustering of cells. Type III tests with Satterthwaite’s approximation were used for inference. Following the mixed-effects analysis, the distribution of puncta within cells was compared across groups. Cells were classified into low <1 , medium 1–5 , or high 5 puncta expression bins. For each region, contingency tables crossing expression bin with Group were analysed using Pearson’s chi-square test. Standardized Pearson residuals were used to identify bins that were over- or under-represented in each group. For image source data, see Supplementary Fig. 2 /articles/s41586-026-10509-4 MOESM1 . Immunofluorescence and analysis Mice were anaesthetized with isoflurane and perfused with cold 1× PBS followed by 4% PFA. Brains were post-fixed in 4% PFA overnight and then transferred to a 30% sucrose/PBS solution for 2 days. Brains were sectioned at 40 µm thickness using a Leica CM3050-S cryostat, with serial sections collected from the VTA. For each subject, 2–3 brain slices were blocked for 2 h 0.1% Triton X-100, 10% normal donkey serum , followed by overnight incubation at 4 °C with primary antibodies: chicken anti-tyrosine hydroxylase 1:500, Aves Labs TYH , rabbit anti-FOS 1:2,000, Synaptic Systems 226-008 , mouse anti-HA 1:1,000, Abcam ab18181 , and/or rabbit anti-H3K4me3Q5dopaminyl 1:250, Millipore ABE2590 . The next day, slices were incubated for 2 h at room temperature with fluorescent-conjugated secondary antibodies 1:1,000, donkey anti-chicken Alexa Fluor 488, Invitrogen A78948; goat anti-rabbit Alexa Fluor 546, Invitrogen A11010; donkey anti-goat Alexa Fluor 680, Invitrogen A10043; donkey anti-mouse Alexa Fluor 680, Invitrogen A10038 . Slices were counterstained with DAPI 1:10,000, Thermo Scientific 62248 and mounted with ProLong Gold Antifade Mountant Thermo Fisher P36934 . Confocal images 2–3 replicates per mouse, 1,024 × 1,024 pixels were acquired on a Zeiss LSM 780 upright microscope using a 40× objective with Zen Black software. Images were averaged across 8 consecutive acquisitions at a bit depth of 16 bits. Image analysis was conducted using FIJI software NIH . The FOS and TH channels were thresholded using the MaxEntropy method to define regions of interest. Following identification of colocalized FOS+/TH+ signal, FOS intensity was measured and averaged from 4–6 images per mouse. Brightness and contast were adjusted for representative images. For image source data, see Supplementary Fig. 3 /articles/s41586-026-10509-4 MOESM1 . Dopamine ELISA Brain tissue dopamine levels were assessed in response to 3 h of pup separation in the home cage. Independent biological replicates were collected at each time point, with separate mice used for each measurement. Brains were rapidly collected see ‘Brain tissue collection’ at baseline prior to pup removal, 0 min , during pup removal 30 and 180 min , and 30 min after pups were returned to the home cage 210 min . Equal weights of micropunched brain tissues were homogenized in lysis buffer 0.01 N HCl, 1 mM EDTA, 4 mM sodium metabisulfite . Tissue dopamine levels were assessed using the Dopamine Research ELISA Kit ALPCO Diagnostics according to manufacturer’s instruction. HPA axis assessment Plasma corticosterone levels were assessed in response to 3 h of pup separation in the home cage. Testing was initiated within 2 h after lights on. Tail blood was collected prior to pup removal 0 min , during pup removal 30 and 180 min , and 120 min after pups were returned to the cage 300 min from the same mice. Samples for control mice in the home cage, used to account for the effects of handling-induced stress, were collected concurrently. Blood samples were immediately mixed with 50 mM EDTA and centrifuged at 5,000 rpm for 10 min. Plasma was collected and stored at −80 °C until analysis. Corticosterone levels were quantified using a Corticosterone ELISA kit ENZO Life Sciences according to manufacturer’s instruction. Clonal TGM2 knockout in HeLa cells HeLa cells ATCC, CCL-2 were cultured at 37 °C with 5% CO2 in DMEM medium high-glucose, ThermoFisher 11965118 supplemented with 10% FBS Sigma-Aldrich and 500 U penicillin and streptomycin. The CRISPR guide RNA targeting exon 5 of TGM2 was purchase from IDT, containing the sequence ACGCTGGGACAACAACTACG. 120 pmol of guide RNA and 100 pmol of Alt-R Streptococcus pyogenes Cas9 Nuclease V3 IDT, 1081058 were premixed for 15 min in 5 μl total with PBS . HeLa cells 200,000 were washed in PBS, before resuspending in 20 μl of nucleofector solution SE Cell Line 4D-Nucleofector X Kit S, Lonza V4XC-1032 . The 5 μl Cas9/sgRNA mix was added, as well as 1 μl of Alt-R Cas9 Electroporation Enhancer IDT, NC1395977 , and all mixed gently. The entire mixture was transferred to a 16-well Nucleocuvette strip Lonza, PDH-2104 gently, and nucleofected using the Lonza 4D-NucleofectorUnit, using the default settings for HeLa cells. Directly after nucleofection, 80 μl of prewarmed culture media was added to the cuvette. The entire mixture was immediately transferred to a 6-well plate with 2 ml of prewarmed culture media, and incubated cultured at 37 °C with 5% CO2 for 48 h. Nucleofection efficiency was assessed by using a second positive control reaction with a pMax-GFP plasmid Addgene, 177825 . The pool of cells was diluted to a concentration of 1 cell per 200 μl, and 100 μl was aliquoted into each well of 10× 96-well plates, and cultured at 37 °C with 5% CO2. After 14 days, single-clones were identified using a light microscope. Single-clone containing wells were expanded and targeting assessed by extracting genomic DNA and performing PCR and PCR sequencing over the targeted site forward: GGCTCCAGCCCCCACCATCTGCCGCAC; reverse: GCCACATAGCGCATTGAGAGTGTTGGT . PCR sequencing results were assessed using Synthego’s ICE analysis. Clones which were identified as introducing premature stop codons were expanded further. Assessment of TGM2 -knockout cell line Western blot for TG2 In brief, HeLa cells were collected and lysed using high-salt buffer 20 mM HEPES pH 7.9, 500 mM KCl, 10 mM MgCl2, and 1 % NP-40 , followed by brief pulse-sonication. One-hundred micrograms of protein was run on a 4–12% Bis-Tris gel Invitrogen, NW04122 for 45 min at 150 V. Protein was transferred to a 0.2 µm nitrocellulose membrane using a Trans-Blot Turbo Transfer System Bio-Rad following the manufacturers protocol. The membrane was blocked in 5% milk in TBS for 1 h, and primary antibody anti-TG2, Abcam 2386, 1:500 in 1.5% milk in TBS overnight at 4 °C. The blot was washed 3× 15 min in TBS-T, and then incubated with secondary antibody Goat anti-Mouse IgG H + L Cross-Adsorbed Secondary Antibody, Alexa Fluor 647, Invitrogen A-21235 at room temperature for 1 h, followed by washing 3× 15 min in TBS-T. The blot was imaged using a Bio-Rad ChemiDoc MP, and then the blot was stained with amido-black stain to assess total protein loading. Transamidation activity assay One-hundred micrograms of HeLa extract was diluted in low-salt buffer 20 mM HEPES pH 7.9, 150 mM KCl, 10 mM MgCl2 , and 1× final transamidation assay buffer was added 25 mM Tris-HCl, pH 8, 10 mM CaCl2, 10 mM DTT, 10 mM KCl . Biotin-cadaverine Millipore Sigma A5348 was added to a final concentration of 1 mM, and reactions incubated at 30 °C for 2 h. A western blot was run as described above, using a Streptavidin-Alexa Fluor 488 Conjugate ThermoFisher S32354 to measure incorporation of biotin-cadaverine. Statistics Statistical analyses for behavioural and immunoassay data were conducted using Prism software GraphPad, v10.4.1 . Data distribution was assessed for normality. Data that met assumptions of normality were analysed using parametric tests, while non-normally distributed data were analysed using non-parametric alternatives. For experiments involving multiple conditions, one-way or two-way ANOVAs were performed, followed by post hoc analyses when appropriate. For time-course analyses where multiple measurements were taken from the same animal, repeated measures ANOVAs were performed. Two-tailed Student’s t -tests were used for comparisons between two conditions. Behavioural data derived from manual observations and pregnancy/litter outcomes were analysed using chi-square tests. Grubb’s test alpha = 0.05 was applied to detect outliers where necessary. Statistical significance was defined as P ≤ 0.05. Bulk RNA-seq and analysis RNA isolation and library preparation Total mRNA was extracted from frozen brain tissues after homogenization in Trizol Reagent Thermo Fisher and cleaned using RNeasy Microcolumns Qiagen following the manufacturer’s instructions. For RNA-seq library preparation, 150 ng of mRNA per sample was used with either the Illumina Stranded mRNA Prep Kit Illumina, 20040534 or the TruSeq RNA Library Prep Kit v2 Illumina, RS-122-2001 , according to the manufacturer’s protocols. Library quality was assessed using a Qubit Fluorometer 2.0 Thermo Fisher and a High Sensitivity D5000 TapeStation assay Agilent before sequencing on a NovaSeq 6000 or NovaSeq X system. Differential expression analysis Raw fastq files, containing an average of 20–30 million reads per sample, were processed for pseudoalignment and abundance quantification using Kallisto v0.46.1 against the Ensembl Mus musculus reference v79 67. To filter lowly expressed genes, only those with a total read count of at least ten across all samples were retained. To account for unwanted variation among samples within each sequencing experiment that could arise from technical or biological factors unrelated to the conditions of interest including litter size, oestrous stage, day of sample collection, etc. , RUVs v1.32.0 was applied with a negative control gene set derived from the total genes identified per sequencing experiment, after ensuring that unwanted variation did not correlate with covariates of interest, as described previously . Differential expression analysis was performed using DESeq2 v1.38.3 68 /articles/s41586-026-10509-4 ref-CR68 , 69 /articles/s41586-026-10509-4 ref-CR69 , with significant genes defined by an adjusted 70 /articles/s41586-026-10509-4 ref-CR70 P value < 0.05. For brain-wide transcriptome comparisons in which samples were processed across multiple sequencing runs and stemming from separate cohorts, pairwise comparisons were performed independently for each brain region. For all other experiments, where individuals came from the same cohort and were processed in a single sequencing run, all groups were analysed together to maintain consistent normalization within the experiment. Gene expression time-course analyses examining the periods before, during, and after pregnancy and postpartum were performed on normalized count data using the ImpulseDE2 package v0.99.10 for each brain region. Significant genes exhibiting transient regulation or monotonous changes in expression were identified using case-only differential expression analysis, with a Q-value threshold of 0.05. To adjust for cell-type composition in the DESeq2 model while avoiding multicollinearity and preserving degrees of freedom, principle components analysis was performed on the scaled surrogate proportion values obtained from BRETIGEA for the entire dataset 71 /articles/s41586-026-10509-4 ref-CR71 . The first two principle components PC1 and PC2 were incorporated into the DESeq2 design formula, with the number of RUVSeq factors 72 /articles/s41586-026-10509-4 ref-CR72 , 73 /articles/s41586-026-10509-4 ref-CR73 k reduced by two to accommodate these additional covariates. WGCNA To identify brain-wide gene co-expression networks, normalized count data for all brain regions were compiled and analysed using the WGCNA package v1.73 74. Co-expression networks were constructed from the 7,500 most variable genes, determined by ranking gene variance. A soft-threshold power of 12 was identified with the pickSoftThreshold function to ensure scale-free network properties. Modules were identified based on dissimilarity of a signed topological overlap matrix, and named with an arbitrary colour. The ‘grey’ module encompassed genes that did not segregate into any specific module, and was therefore removed from further analyses. To assess the enrichment of DEGs within each brain region per module, Fisher’s exact tests were performed using the fisher.test function in R v4.3.0 . To analyse the correlation between gene modules and brain regional sensitivity to parity, brain regions were categorized into two groups based on the top five regions with significant DEG overlap. These regions were designated as high-sensitivity regions, while all other regions were classified as low-sensitivity. Pearson correlation coefficients were calculated to examine the relationship between each module and the trait regional sensitivity × parity status , with statistical significance determined by Student’s t -tests P < 0.05 . Heat maps were generated to visualize these module–trait correlations, and module-specific gene lists were exported for pathway enrichment analysis. Pathway and predicted upstream regulator analyses Functional annotation of DEGs was conducted using ShinyGO v0.81 75, with all protein-coding genes in the mm10 genome used as the background. Pathways with an FDR < 0.05 were considered significant. Relevant pathways were selected from the top 10 or one-third of significant terms, ranked by FDR, to emphasize processes consistent with hypotheses informed by published literature. For GO term selection, Revigo was used to reduce redundancy of overlapping GO terms when needed. Ingenuity Pathway Analysis IPA; Qiagen v01-23-01 was used to predict upstream regulators for DEG lists 76 /articles/s41586-026-10509-4 ref-CR76 . For high-sensitivity and low-sensitivity regions see ‘WGCNA’ section , DEGs shared by at least two or three brain regions were extracted, irrespective of the direction of change. The IPA software was used to identify upstream regulators associated with each DEG list, with statistical significance defined as 77 /articles/s41586-026-10509-4 ref-CR77 P < 0.05. To investigate the mechanisms underlying parity programming of dHF plasticity, significant upstream regulators were systematically prioritized based on one or more of the following criteria: 1 multiple molecules involved in the same signalling pathways for example, progesterone and its receptor, PGR ; 2 structurally similar molecules that engage analogous signalling cascades for example, levodopa and dopamine ; 3 molecules known to be influenced by reproductive exposures; 4 molecules demonstrated to play a role in hippocampal plasticity; 5 molecules identified as significant in both high-sensitivity and low-sensitivity regulator lists; and 6 molecules with high statistical significance values. Following selection, upstream regulators were categorized into general molecular classes for example, ‘hormone’, ‘transcription factor’ or ‘lipid’ and grouped for visualization. Data are presented as a bubble plot, with marker shape representing number of brain regions sharing DEGs and marker size corresponding to significance. Gene expression overlap analyses Jaccard indices were calculated for overlapping gene lists using the GeneOverlap package v1.36.0 , with significance defined by p < 0.05. Fisher’s exact tests were calculated using the base fisher.test function in R v4.3.0 following construction of contingency tables for each comparison. Transcriptome-wide, threshold-free gene expression overlap was visualized using RRHO heat maps generated with the RRHO2 package v1.0 78. Gene lists were ranked by signed P values, calculated as the log 10-transformed nominal P value multiplied by the sign of the fold change, without applying differential expression thresholds. Cell-type deconvolution Brain cell-type proportion was estimated from normalized expression data using the BRETIGEA package v1.0.4 79 with default markers. For comparisons between groups, surrogate proportion variables were normalized using the Normalize function in Prism GraphPad v10.4.1 . Each sub column was normalized separately, with the smallest value set to 0% and the largest value set to 100%. To calculate log 2 fold change between groups, normalized values per cell type per group were averaged and expressed as the log 2-transformed ratio of the group means. snRNA-seq Nuclei isolation and library preparation For each mouse, 2 mm dHF-enriched tissue micropunches were collected bilaterally from consecutive 1 mm slices from −0.80 to −2.80 mm relative to bregma for a total of 4 punches Ted Pella . Samples were processed in batches of 4–6, with each group represented within each batch. Nuclei were isolated using a modified version of a sucrose density gradient isolation protocol 80. In brief, thawed tissues were placed in 1 ml of lysis buffer 0.32 M sucrose, 5 mM CaCl 2, 3 mM magnesium acetate, 0.1 mM EDTA, 10 mM Tris-HCl pH 8, 1 mM DTT, 0.1% Triton X-100 with 50 μl of 25 U ml −1RNase inhibitor Takara 2313B in a dounce homogenizer Wheaton 357538 . Homogenization was performed with 20 strokes using a tight pestle. Another 1 ml of lysis buffer was added, followed by an additional 10 strokes. The resulting 2 ml homogenate was transferred to a 15 ml Open-Top Thinwall Polypropylene Tube Beckman-Coulter 361707 . The homogenizer and pestle were rinsed with 2 ml of lysis buffer, and this wash was combined with the homogenate for a total of 4 ml. The homogenate was carefully underlaid with 9 ml of sucrose solution 1.8 M sucrose, 3 mM magnesium acetate, 1 mM DTT, 10 mM Tris-HCl pH 8 and ultracentrifuged at 24,000 rpm for 1 h at 4 °C using a Sorvall WX+ centrifuge. After centrifugation, the supernatant and debris at the interface were gently removed. The nuclear pellet was resuspended in 1 ml of resuspension buffer 0.02% bovine serum albumin in DPBS with 25 μl of 25 U ml −1RNase inhibitor and incubated on ice for 10 min. The suspension was passed through a 35 μm nylon mesh filter Corning 352235 into a 1.5 ml RNase/DNase-free microcentrifuge tube and centrifuged at 2,600 g for 10 min at 4 °C. Supernatants were discarded, and nuclei were resuspended in 200 μl of resuspension buffer. A 10 μl aliquot of the nuclei suspension was stained with Trypan Blue to assess quality and concentration using a haemocytometer. Nuclei suspensions were loaded onto a Chromium Single Cell 3′ chip 10X Genomics, v3 and processed according to the manufacturer’s protocol, targeting 10,000 nuclei. Single-nuclei libraries were generated using the 10X Chromium Next GEM Single Cell 3′ v3.1 Dual Index protocol CG000315 Rev A . Libraries were pooled, loaded onto a single 10B 100 Cycle Flowcell, and sequenced using an Illumina NovaSeq 6000 system to generate 25,000–30,000 paired-end 2× 100 bp reads per cell. Data analysis FastQ files were processed with the 10X Genomics Cell Ranger pipeline v7.1.0 to demultiplex reads, align them to the mouse genome mm10-2020-A , remove PCR duplicates, and generate gene expression matrices. Cell Ranger filtered outputs were analysed using Seurat v4.3.0 81, and mitochondrial RNA content per cell was calculated using the GRCm39 mm10 genome annotation and regressed out using SCTransform normalization protocol included in the Seurat toolkit with 20 principal components and a resolution of 0.1. To estimate ambient RNA and correct for background contamination, the SoupX v1.6.2 package was used for each sample using raw and filtered feature matrices from the Cell Ranger output. Heterotypic doublets were identified and removed using DoubletFinder v2 82 /articles/s41586-026-10509-4 ref-CR82 to ensure the integrity of singlet datasets. Filtered singlet datasets were then re-normalized and integrated using the same Seurat SCTransform v2 workflow mentioned above. Cell clusters were annotated using a combination of expert curation based on published marker genes 83 /articles/s41586-026-10509-4 ref-CR83 , and label transfer from hippocampal reference datasets, including the Allen Brain Map and Broad Institute resources 84 ref-CR84 , 85 ref-CR85 , 86 ref-CR86 , 87 ref-CR87 , 88 ref-CR88 , 89 ref-CR89 , 90 /articles/s41586-026-10509-4 ref-CR90 . Clusters with contaminant cell populations expressing markers for choroid plexus 89 /articles/s41586-026-10509-4 ref-CR89 , 91 /articles/s41586-026-10509-4 ref-CR91 Ttr , ependymal Tmem212 , and vascular leptomeningeal cells Vtn and Col1a2 were removed from the analysis . Additionally, as the sequential 2 mm micropunches encompassed portions of cortical, thalamic and vHF regions, clusters characterized by enrichment of published non-dHF neuronal markers using Seurat’s FindMarkers function layer 5/6 cortical: 90 /articles/s41586-026-10509-4 ref-CR90 , 92 ref-CR92 , 93 ref-CR93 , 94 /articles/s41586-026-10509-4 ref-CR94 Rorb and Foxp2 ; ventral granule neurons: 94 /articles/s41586-026-10509-4 ref-CR94 Tox3 were also removed from the analysis. Cell cluster proportion analyses were conducted using the scProportionTest package 86 /articles/s41586-026-10509-4 ref-CR86 , which employs a Monte Carlo permutation test to evaluate whether observed differences result from random sampling. Proportional differences between conditions were compared to a null distribution generated by resampling, and statistical significance was determined by permutation-based 95 /articles/s41586-026-10509-4 ref-CR95 P values with confidence intervals estimated via bootstrapping. Differential expression analysis was conducted using pseudobulk analysis, where gene counts were summed across all cells within each sample for each cell-type cluster using the AggregateExpression function. DESeq2 was then applied at the sample level to conduct differential expression. Note that differential expression analysis could not be performed for the GABA.4 cluster due to insufficient sample representation in multiple groups. To explore pathways underlying cluster-specific differences across conditions, pathway analysis was conducted using ShinyGO on genes meeting the following criteria: log 2 fold change 1.5 and P < 0.05. CUT&RUN-seq The procedure was adapted from established protocols 33,96. For mouse brain tissue, two 1.5-mm dHF punches were used for each biological replicate and split across two reactions H3K4me3Q5dop and IgG . Samples were not pooled. For human brain tissue, ~10 mg of tissue dissected from individual frozen subiculum samples was collected from each subject, and split across the indicated antibodies. For HeLa cells, samples were split across H3K4me3Q5dop, H3K4me3 and IgG reactions. Samples were dounce homogenized in nuclear extract NE buffer 20 mM HEPES-KOH pH 7.9, 10 mM KCl, 0.5 mM spermidine, 0.1% Triton X-100, 20% glycerol with protease inhibitors and passed through a 21 gauge needle 10 times. Nuclei were pelleted at 1,100 g for 5 min at 4 °C in a swinging-bucket rotor, passed through a 40-μm strainer pluriSelect , washed again in 500 μl NE buffer and counted. BioMag Plus Concanavalin A beads Polysciences were prepared per reaction by washing three times with binding buffer 20 mM HEPES-KOH pH 7.9, 10 mM KCl, 1 mM CaCl 2, 1 mM MnCl 2 . Beads 15 μl were aliquoted into 1.7 ml DNA low-bind tubes Eppendorf containing 500 μl NE buffer and 100,000 nuclei per reaction. Samples were rotated at room temperature for 10 min, then bead-bound nuclei were washed three times with wash buffer 20 mM HEPES pH 7.5, 150 mM NaCl, 0.1% Triton X-100, 0.1% Tween-20, 0.5 mM spermidine, 0.1% BSA, and protease inhibitors , resuspended in 100 μl antibody buffer wash buffer with 2 mM EDTA , and mixed. 2 μl antibodies were added to the corresponding tubes: H3K4me3 Active Motif, 39159 , H3K4me3Q5dopaminyl Millipore, ABE2590 , or rabbit IgG Invitrogen, 10500c . Samples were incubated overnight at 4 °C on a rotating mixer angled upward at 20 degrees. The next day, nuclei were washed twice with cold wash buffer, and incubated with 2.5 μl pAG-MNase Epicypher, 15-1016 for 1 h at 4 °C. After 4 washes with cold wash buffer and one with low-salt rinse buffer 20 mM HEPES pH 7.5, 0.5 mM spermidine, 0.1% Tween-20, 0.1% Triton X-100 , nuclei were resuspended in calcium incubation buffer 3.5 mM HEPES pH 7.5, 10 mM CaCl 2, 0.1% Tween-20, 0.1% Triton X-100 and placed into an ice-cold block at 4 °C. MNase digestion was stopped by adding 100 μl of 2× Stop Buffer 340 mM NaCl, 20 mM EDTA, 5 mM EGTA, 0.1% Tween-20, 0.1% Triton X-100, 25 μg ml −1RNase A, and 0.05 ng per 100 μl Escherichia coli spike-in DNA , followed by a 15 min incubation at 37 °C without shaking. Beads were then placed on a magnet and 200 μl of supernatant was collected. DNA was purified using the Zymo ChIP DNA Clean & Concentrator kit D5205 , eluted in 30 μl, and stored at −20 °C for library preparation. Libraries were generated using the NEBNext Ultra II DNA library kit, quantified using a Qubit fluorometer with the HS DNA kit, checked for size distribution on the Agilent TapeStation, pooled equimolarly, and sequenced on an Illumina NovaSeq X. Data analysis Raw fastq files were aligned to the hg19 or mm10 genome using bowtie2 v2.5.0 97. Low-quality reads were filtered using Samtools v1.9 with a MAPQ cut-off score of 30 . Only unique, deduplicated reads were retained for further processing. Bigwig files were produced using the deepTools package v3.5.1 , using an ENCODE hg19 or mm10 v2 blacklist file to discard regions with consistently non-specific signal, and scaled using 98 /articles/s41586-026-10509-4 ref-CR98 E. coli spike-in controls to normalize sequencing depth. To determine normalization factors based on E. coli reads, each sample was aligned to the E. coli genome MG1655 , and the unique deduplicated reads were compared across groups per antibody per experiment. The sample with the lowest number of E. coli reads was determined minimum , and all samples were scaled by dividing their corresponding E. coli read count by this minimum number . For each group, bigwig files were merged and peak calling was conducted using MACS2 v2.1.0 with the corresponding merged IgG file as control, and filtered for peaks with FDR < 0.05 99 /articles/s41586-026-10509-4 ref-CR99 . Peak annotation was conducted using HOMER v4.1.1 100 /articles/s41586-026-10509-4 ref-CR100 . Heat maps were made either using the DiffBind v3.8.4 or deepTools v3.5.5 packages 101 /articles/s41586-026-10509-4 ref-CR101 . For deepTools, heat maps were made by merging DEGs from RNA-seq data with TSSs downloaded from the UCSC genome browser using the canonically annotated transcript for each gene. Profiles were generated and statistically analysed using the deepStats package 102 /articles/s41586-026-10509-4 ref-CR102 by using the dsCompareCurves function to perform Wilcoxon Rank-sum tests per-bin. For DiffBind analysis, heat maps were made for peaks identified by DiffBind’s differential peak algorithm, where differential peaks were first filtered using a log 103 /articles/s41586-026-10509-4 ref-CR103 2 fold change threshold 0.1 and defined at p < 0.05, where log 2 fold change was calculated as log 2 parity − log 2 NP , based on prior empirical observations used to define thresholds for differential peaks . ChEA analysis on annotated loci was conducted using EnrichR with a significance threshold of adjusted 36 /articles/s41586-026-10509-4 ref-CR36 P < 0.05 ref. . 104 /articles/s41586-026-10509-4 ref-CR104 Reporting summary Further information on research design is available in the Nature Portfolio Reporting Summary /articles/s41586-026-10509-4 MOESM2 linked to this article. Data availability The genomics data generated in this study have been deposited in the National Center for Biotechnology Information Gene Expression Omnibus GEO database under the SuperSeries GSE298544 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE298544 . Reference genomes used in this study include the mouse genome mm10 and human genome hg19. Publicly available reference datasets used for annotation included resources from the Allen Brain Map and Broad Institute, with gene annotations obtained from the UCSC Genome Browser and ENCODE blacklist regions applied where appropriate. All data supporting the findings of this study are available within the article and Supplementary Information /articles/s41586-026-10509-4 MOESM1 . Related data, including raw microscopy images, are available from the corresponding authors upon reasonable request. Source data /articles/s41586-026-10509-4 Sec53 are provided with this paper. Code availability Data were analysed using publicly available software and standard pipelines, as described in the Methods /articles/s41586-026-10509-4 Sec10 . No custom code was generated for this study; however, any related analysis scripts or details are available from the corresponding authors upon reasonable request. References Hoekzema, E. et al. Pregnancy leads to long-lasting changes in human brain structure. Nat. Neurosci. 20 , 287–296 2017 .Martínez-García, M. et al. Do pregnancy-induced brain changes reverse? The brain of a mother six years after parturition. Brain Sci. 11 , 168 2021 .Voldsbekk, I. et al. A history of previous childbirths is linked to women’s white matter brain age in midlife and older age. Hum. Brain Mapp. 42 , 4372–4386 2021 .Pritschet, L. et al. Neuroanatomical changes observed over the course of a human pregnancy. Nat. Neurosci. 27 , 2253–2260 2024 .de Lange, A.-M. G. et al. Population-based neuroimaging reveals traces of childbirth in the maternal brain. Proc. Natl Acad. Sci. USA 116 , 22341–22346 2019 .de Lange, A.-M. G. et al. The maternal brain: region-specific patterns of brain aging are traceable decades after childbirth. Hum. Brain Mapp. 41 , 4718–4729 2020 .Pawluski, J. L. & Galea, L. A. M. Hippocampal morphology is differentially affected by reproductive experience in the mother. J. Neurobiol. 66 , 71–81 2006 .Duarte-Guterman, P. et al. Cellular and molecular signatures of motherhood in the adult and ageing rat brain. Open Biol. 13 , 230217 2023 .Workman, J. L. et al. Parity modifies the effects of fluoxetine and corticosterone on behavior, stress reactivity, and hippocampal neurogenesis. Neuropharmacology 105 , 443–453 2016 .Macbeth, A. H., Scharfman, H. E., Maclusky, N. J., Gautreaux, C. & Luine, V. N. Effects of multiparity on recognition memory, monoaminergic neurotransmitters, and brain-derived neurotrophic factor BDNF . Hormones Behav. 54 , 7–17 2008 .Kinsley, C. H. & Lambert, K. G. Reproduction-induced neuroplasticity: natural behavioural and neuronal alterations associated with the production and care of offspring. J. Neuroendocrinol. 20 , 515–525 2008 .Barha, C. K., Lieblich, S. E., Chow, C. & Galea, L. A. M. Multiparity-induced enhancement of hippocampal neurogenesis and spatial memory depends on ovarian hormone status in middle age. Neurobiol. Aging 36 , 2391–2405 2015 .Love, G. et al. Maternal experience produces long-lasting behavioral modifications in the rat. Behav. Neurosci. 119 , 1084–1096 2005 .Ammari, R. et al. Hormone-mediated neural remodeling orchestrates parenting onset during pregnancy. Science 382 , 76–81 2023 .Bridges, R. S. Long-term effects of pregnancy and parturition upon maternal responsiveness in the rat. Physiol. Behav. 14 , 245–249 1975 .Pawluski, J. L., Vanderbyl, B. L., Ragan, K. & Galea, L. A. M. First reproductive experience persistently affects spatial reference and working memory in the mother and these effects are not due to pregnancy or ‘mothering’ alone. Behav. Brain Res. 175 , 157–165 2006 .Pawluski, J. L., Walker, S. K. & Galea, L. A. M. Reproductive experience differentially affects spatial reference and working memory performance in the mother. Horm. Behav. 49 , 143–149 2006 .Lemaire, V. et al. Motherhood-induced memory improvement persists across lifespan in rats but is abolished by a gestational stress. Eur. J. Neurosci. 23 , 3368–3374 2006 .Gatewood, J. D. et al. Motherhood mitigates aging-related decrements in learning and memory and positively affects brain aging in the rat. Brain Res. Bull. 66 , 91–98 2005 .Macbeth, A. H., Gautreaux, C. & Luine, V. N. Pregnant rats show enhanced spatial memory, decreased anxiety, and altered levels of monoaminergic neurotransmitters. Brain Res. 1241 , 136–147 2008 .Tomizawa, K. et al. Oxytocin improves long-lasting spatial memory during motherhood through MAP kinase cascade. Nat. Neurosci. 6 , 384–390 2003 .Paris, J. J. & Frye, C. A. Estrous cycle, pregnancy, and parity enhance performance of rats in object recognition or object placement tasks. Reproduction 136 , 105–115 2008 .McArdle, H. J. The metabolism of copper during pregnancy—a review. Food Chem. 54 , 79–84 1995 .Leuner, B. & Sabihi, S. The birth of new neurons in the maternal brain: Hormonal regulation and functional implications. Front. Neuroendocrinol. 41 , 99–113 2016 .Le, A. A. et al. Prepubescent female rodents have enhanced hippocampal LTP and learning relative to males, reversing in adulthood as inhibition increases. Nat. Neurosci. 25 , 180–190 2022 .Gangarossa, G. et al. Characterization of dopamine D 1and D2receptor-expressing neurons in the mouse hippocampus. Hippocampus 22 , 2199–2207 2012 .Puighermanal, E. et al. Anatomical and molecular characterization of dopamine D 1receptor-expressing neurons of the mouse CA1 dorsal hippocampus. Brain Struct. Funct. 222 , 1897–1911 2017 .Godino, A. et al. Dopamine D 1–D2 signalling in hippocampus arbitrates approach and avoidance. Nature 643 , 448–457 2025 .Hansen, S., Bergvall, ÅH. & Nyiredi, S. Interaction with pups enhances dopamine release in the ventral striatum of maternal rats: A microdialysis study. Pharmacol. Biochem. Behav. 45 , 673–676 1993 .Champagne, F. A. et al. Variations in nucleus accumbens dopamine associated with individual differences in maternal behavior in the rat. J. Neurosci. 24 , 4113–4123 2004 .Farrelly, L. A. et al. Histone serotonylation is a permissive modification that enhances TFIID binding to H3K4me3. Nature 567 , 535–539 2019 .Lepack, A. E. et al. Dopaminylation of histone H3 in ventral tegmental area regulates cocaine seeking. Science 368 , 197–201 2020 .Zheng, Q. et al. Bidirectional histone monoaminylation dynamics regulate neural rhythmicity. Nature 637 , 974–982 2025 .Chan, J. C. & Maze, I. Nothing is yet set in Hi stone: novel post-translational modifications regulating chromatin function. Trends Biochem. Sci 45 , 829–844 2020 .Al-Kachak, A. et al. Histone serotonylation in dorsal raphe nucleus contributes to stress- and antidepressant-mediated gene expression and behavior. Nat. Commun. 15 , 5042 2024 .Chan, J. C. et al. Serotonin transporter-dependent histone serotonylation in placenta contributes to the neurodevelopmental transcriptome. J. Mol. Biol. 436 , 168454 2024 .Sardar, D. et al. Induction of astrocytic Slc22a3 regulates sensory processing through histone serotonylation. Science 380 , eade0027 2023 .Stewart, A. F., Lepack, A. E., Fulton, S. L., Safovich, P. & Maze, I. Histone H3 dopaminylation in nucleus accumbens, but not medial prefrontal cortex, contributes to cocaine-seeking following prolonged abstinence. Mol. Cell Neurosci. 125 , 103824 2023 .Fulton, S. L. et al. Histone H3 dopaminylation in ventral tegmental area underlies heroin-induced transcriptional and behavioral plasticity in male rats. Neuropsychopharmacology 47 , 1776–1783 2022 .Chen, M. et al. Quinone reductase 2 reads H3 serotonylation to support neuronal maturation. Preprint at bioRxiv https://doi.org/10.64898/2026.03.17.712426 https://doi.org/10.64898/2026.03.17.712426 2026 .Pawluski, J. L., Lambert, K. G. & Kinsley, C. H. Neuroplasticity in the maternal hippocampus: Relation to cognition and effects of repeated stress. Hormones Behav. 77 , 86–97 2016 .Kinsley, C. H. et al. Motherhood induces and maintains behavioral and neural plasticity across the lifespan in the rat. Arch. Sex Behav. 37 , 43–56 2008 .Jeffery, K. J. Integration of the sensory inputs to place cells: What, where, why, and how?. Hippocampus 17 , 775–785 2007 .Lisman, J. E. & Grace, A. A. The hippocampal–VTA loop: controlling the entry of information into long-term memory. Neuron 46 , 703–713 2005 .Kimble, D. P., Rogers, L. & Hendrickson, C. W. Hippocampal lesions disrupt maternal, not sexual, behavior in the albino rat. J. Comp. Physiol. Psychol. 63 , 401–407 1967 .Terlecki, L. J. & Sainsbury, R. S. Effects of fimbria lesions on maternal behavior in the rat. Physiol. Behav. 21 , 89–97 1978 .Numan, M. et al. The effects of D 1or D2dopamine receptor antagonism in the medial preoptic area, ventral pallidum, or nucleus accumbens on the maternal retrieval response and other aspects of maternal behavior in rats. Behav. Neurosci. 119 , 1588–1604 2005 .Neukel, C. et al. The maternal brain in women with a history of early-life maltreatment: an imagination-based fMRI study of conflictual versus pleasant interactions with children. J. Psychiatry Neurosci. 43 , 273–282 2018 .Kim, P. et al. Perceived quality of maternal care in childhood and structure and function of mothers’ brain. Dev. Sci. 13 , 662–673 2010 .Deems, N. P. & Leuner, B. Pregnancy, postpartum and parity: resilience and vulnerability in brain health and disease. Front. Neuroendocrinol. 57 , 100820 2020 .Dumais, A. et al. Risk factors for suicide completion in major depression: a case-control study of impulsive and aggressive behaviors in men. Am. J. Psychiatry 162 , 2116–2124 2005 .Dumais, A. et al. Is violent method of suicide a behavioral marker of lifetime aggression?. Am. J. Psychiatry 162 , 1375–1378 2005 .Mai, J. K., Majtanik, M. & Paxinos, G. eds . Atlas of the Human Brain , 4th Edn Elsevier, 2016 .Rice, C. J., Sandman, C. A., Lenjavi, M. R. & Baram, T. Z. A novel mouse model for acute and long-lasting consequences of early life stress. Endocrinology 149 , 4892–4900 2008 .Gilles, E. E., Schultz, L. & Baram, T. Z. Abnormal corticosterone regulation in an immature rat model of continuous chronic stress. Pediatr. Neurol. 15 , 114–119 1996 .Peña, C. J. et al. Early life stress confers lifelong stress susceptibility in mice via ventral tegmental area OTX2. Science 356 , 1185–1188 2017 .Fleming, A. S. & Rosenblatt, J. S. Maternal behavior in the virgin and lactating rat. J. Comp. Physiol. Psychol. 86 , 957–972 1974 .Leuner, B., Glasper, E. R. & Gould, E. Parenting and plasticity. Trends Neurosci. 33 , 465–473 2010 .Bridges, R. S. Neuroendocrine regulation of maternal behavior. Front. Neuroendocrinol. 0 , 178–196 2015 .Krishnan, V. et al. Molecular adaptations underlying susceptibility and resistance to social defeat in brain reward regions. Cell 131 , 391–404 2007 .Champagne, F. A., Curley, J. P., Keverne, E. B. & Bateson, P. P. G. Natural variations in postpartum maternal care in inbred and outbred mice. Physiol. Behav. 91 , 325–334 2007 .Stolzenberg, D. S. & Rissman, E. F. Oestrogen-independent, experience-induced maternal behaviour in female mice. J. Neuroendocrinol. 23 , 345–354 2011 .Carcea, I. et al. Oxytocin neurons enable social transmission of maternal behaviour. Nature 596 , 553–557 2021 .Wolcott, N. S. et al. Automated classification of estrous stage in rodents using deep learning. Sci. Rep. 12 , 17685 2022 .Nagai, Y. et al. Deschloroclozapine, a potent and selective chemogenetic actuator enables rapid neuronal and behavioral modulations in mice and monkeys. Nat. Neurosci. 23 , 1157–1167 2020 .Secci, M. E., Reed, T., Quinlan, V., Gilpin, N. W. & Avegno, E. M. Quantitative analysis of gene expression in RNAscope-processed brain tissue. Bio. Protoc. 13 , e4580 2023 .Bray, N. L., Pimentel, H., Melsted, P. & Pachter, L. Near-optimal probabilistic RNA-seq quantification. Nat. Biotechnol. 34 , 525–527 2016 .Risso, D., Ngai, J., Speed, T. P. & Dudoit, S. Normalization of RNA-seq data using factor analysis of control genes or samples. Nat. Biotechnol. 32 , 896–902 2014 .Peixoto, L. et al. How data analysis affects power, reproducibility and biological insight of RNA-seq studies in complex datasets. Nucleic Acids Res. 43 , 7664–7674 2015 .Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15 , 550 2014 .Fischer, D. S., Theis, F. J. & Yosef, N. Impulse model-based differential expression analysis of time course sequencing data. Nucleic Acids Res. 46 , e119 2018 .Quinn, T. P., Erb, I., Richardson, M. F. & Crowley, T. M. Understanding sequencing data as compositions: an outlook and review. Bioinformatics 34 , 2870–2878 2018 .Gural, B. et al. Novel insights into post-myocardial infarction cardiac remodeling through algorithmic detection of cell-type composition shifts. PLoS Genet. 21 , e1011807 2025 .Langfelder, P. & Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9 , 559 2008 .Ge, S. X., Jung, D. & Yao, R. ShinyGO: a graphical gene-set enrichment tool for animals and plants. Bioinformatics 36 , 2628–2629 2020 .Supek, F., Bošnjak, M., Škunca, N. & Šmuc, T. REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS ONE 6 , e21800 2011 .Krämer, A., Green, J., Pollard, J. Jr & Tugendreich, S. Causal analysis approaches in Ingenuity Pathway Analysis. Bioinformatics 30 , 523–530 2014 .Cahill, K. M., Huo, Z., Tseng, G. C., Logan, R. W. & Seney, M. L. Improved identification of concordant and discordant gene expression signatures using an updated rank-rank hypergeometric overlap approach. Sci. Rep. 8 , 9588 2018 .McKenzie, A. T. et al. Brain cell type specific gene expression and co-expression network architectures. Sci. Rep. 8 , 8868 2018 .Matevossian, A. & Akbarian, S. Neuronal nuclei isolation from human postmortem brain tissue. J. Vis. Exp. https://doi.org/10.3791/914 https://doi.org/10.3791/914 2008 .Hao, Y. et al. Integrated analysis of multimodal single-cell data. Cell 184 , 3573–3587.e29 2021 .Young, M. D. & Behjati, S. SoupX removes ambient RNA contamination from droplet-based single-cell RNA sequencing data. Gigascience 9 , giaa151 2020 .McGinnis, C. S., Murrow, L. M. & Gartner, Z. J. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. 8 , 329–337.e4 2019 .Artegiani, B. et al. A single-cell RNA sequencing study reveals cellular and molecular dynamics of the hippocampal neurogenic niche. Cell Rep. 21 , 3271–3284 2017 .Zhou, Z., Zhi, C., Chen, D., Cai, Z. & Jiang, X. Single-nucleus RNA sequencing reveals cell type-specific transcriptome alterations of Down syndrome hippocampus using the Dp16 mouse model. Genes Genomics 45 , 1305–1315 2023 .Cembrowski, M. S., Wang, L., Sugino, K., Shields, B. C. & Spruston, N. Hipposeq: a comprehensive RNA-seq database of gene expression in hippocampal principal neurons. eLife 5 , e14997 2016 .Harris, K. D. et al. Classes and continua of hippocampal CA1 inhibitory neurons revealed by single-cell transcriptomics. PLoS Biol. 16 , e2006387 2018 .Haslinger, A., Schwarz, T. J., Covic, M. & Lie, D. C. Expression of Sox11 in adult neurogenic niches suggests a stage-specific role in adult neurogenesis. Eur. J. Neurosci. 29 , 2103–2114 2009 .Yao, Z. et al. A taxonomy of transcriptomic cell types across the isocortex and hippocampal formation. Cell 184 , 3222–3241.e26 2021 .Bandler, R. C. et al. Single-cell delineation of lineage and genetic identity in the mouse brain. Nature 601 , 404–409 2022 .Habib, N. et al. Div-seq: single-nucleus RNA-seq reveals dynamics of rare adult newborn neurons. Science 353 , 925–928 2016 .Olney, K. C. et al. Widespread choroid plexus contamination in sampling and profiling of brain tissue. Mol. Psychiatry 27 , 1839–1847 2022 .Xiang, M. et al. A single-cell transcriptional roadmap of the mouse and human lymph node lymphatic vasculature. Front. Cardiovasc. Med. 7 , 52 2020 .Tasic, B. et al. Adult mouse cortical cell taxonomy by single cell transcriptomics. Nat. Neurosci. 19 , 335–346 2016 .Miller, S. A. et al. LSD1 and aberrant DNA methylation mediate persistence of enteroendocrine progenitors that support BRAF mutant colorectal cancer. Cancer Res. 81 , 3791–3805 2021 .Skene, P. J. & Henikoff, S. An efficient targeted nuclease strategy for high-resolution mapping of DNA binding sites. eLife 6 , e21856 2017 .Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9 , 357–359 2012 .Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25 , 2078–2079 2009 .Qasim, M. N. et al. Genome-wide profiling of transcription factor–DNA binding interactions in Candida albicans : a comprehensive CUT&RUN method and data analysis workflow. J. Vis. Exp. https://doi.org/10.3791/63655 https://doi.org/10.3791/63655 2022 .Zhang, Y. et al. Model-based analysis of ChIP–seq MACS . Genome Biol. 9 , R137 2008 .Heinz, S. et al. Simple combinations of lineage-determining transcription factors prime cis -regulatory elements required for macrophage and B cell identities. Mol. Cell 38 , 576–589 2010 .Ramírez, F., Dündar, F., Diehl, S., Grüning, B. A. & Manke, T. deepTools: a flexible platform for exploring deep-sequencing data. Nucleic Acids Res. 42 , W187–W191 2014 .Gautier, R. gtrichard/deepStats: deepStats 0.3.1. Zenodo https://doi.org/10.5281/zenodo.3361799 https://doi.org/10.5281/zenodo.3361799 2019 .Lachmann, A. et al. ChEA: transcription factor regulation inferred from integrating genome-wide ChIP-X experiments. Bioinformatics 26 , 2438–2444 2010 . Acknowledgements We thank members of the Maze laboratory for their helpful discussion on this study; members of the Nestler, Russo and Kenny laboratories for providing materials and/or access to behavioural equipment; and M. Chen and E. Andraka for experimental assistance. I.M. discloses support for the research of this work from the National Institutes of Health R01 MH116900 and the Howard Hughes Medical Institute. J.C.O. discloses support for the research of this work from the National Institutes of Health F32 MH126534, K01 MH139990 and the Brain and Behavior Research Foundation. A.M.C. discloses support from the National Institutes of Health F31 NS132558, F99 NS139541 . B.H.W. discloses support from the National Institutes of Health F32 MH140478 . W.C. discloses support for the research of this work from the National Institutes of Health F30 AG090096 . All other authors declare no relevant funding. Author information Authors and Affiliations Contributions J.C.O. and I.M. conceived of the study, designed the experiments and interpreted the data. J.C.O., G.D.S., A.M.C. and S.D. performed behavioural experiments. J.C.O., G.D.S., A.M.C., S.D. and W.C. performed surgeries. J.C.O., G.D.S., R.R.I. and C.Z. performed oestrous cycle analyses. J.C.O., G.D.S., A.M.C. and S.D. conducted immunofluorescence and RNAscope analyses. J.C.O., G.D.S. and B.H.W. performed molecular experiments. J.C.O. and E.B. performed single-cell sequencing experiments. J.C.O., G.D.S., E.W. and B.H.W. performed bioinformatic analyses. N.M. and G.T. provided human tissues. J.C.O. and I.M. wrote the manuscript. Corresponding authors Ethics declarations Competing interests The authors declare no competing interests. Peer review Peer review information Nature thanks Tatiana Kutateladze and the other, anonymous, reviewer s for their contribution to the peer review of this work. Peer review reports /articles/s41586-026-10509-4 MOESM4 are available. Additional information Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Extended data figures and tables Extended Data Fig. 1 Brain-wide transcriptomic analysis of cell markers and gene expression networks. /articles/s41586-026-10509-4/figures/6 a Brain regions selected for transcriptional profiling: olfactory bulb OB , medial prefrontal cortex mPFC , nucleus accumbens NAc , medial preoptic area mPOA , paraventricular nucleus PVN , dorsal hippocampal formation dHF , basolateral amygdala BLA , ventral hippocampus vHF , ventral tegmental area VTA , dorsal raphe nucleus DRN , and locus coeruleus LC . Yellow circles indicate the micropunched area. Created in BioRender; Lab, M. https://biorender.com/4f8zyrj https://biorender.com/4f8zyrj 2026 . b-f Linear regression analyses between fold change in normalized surrogate proportion variables generated from cell-type deconvolution analyses for b astrocytic, c microglia, d endothelial, e mature oligodendrocyte, and f oligodendrocyte precursor cell markers with the number of DEGs identified from each brain region. g DEGs per brain region following adjustment for cell composition changes padj <0.05 . h-i Enriched GO Biological Processes in dHF between NP vs . RE, ranked by top 10 FDR values, for h baseline unadjusted analyses and i cell composition-adjusted analyses; with the number of genes in each GO term indicated. j Heatmap of dHF DEGs between NP vs . RE, stratified by litter size. k Scale independence plot for soft-threshold power selection for WGCNA analysis, illustrating the relationship between β and the scale-free topology fit index, identifying the optimal β value = 12. l Mean connectivity plot illustrating the average network density for each evaluated soft-thresholding power. m Dendrogram of gene clustering based on topological overlap. n Dendrogram of eigengene network adjacency depicting similarity of modules. o-p Significantly enriched o GO Biological Processes and p KEGG pathways, sorted by FDR values, from module gene sets. Bars are color-coded to indicate modules. The total number of significant terms for each analysis is provided. Extended Data Fig. 2 RE-associated behavioral adaptations are not influenced by locomotion, anxiety-/depressive-like behaviors, or estrous stage. /articles/s41586-026-10509-4/figures/7 a-d In the open field test, there was no change in total distance, time spent in the center, or latency to the center. While there was a significant increase in time spent in corners for RE dams unpaired two-sided Student’s t-test, p = 0.0382 , this did not impact any other outcomes on the open field test. e-g In the light-dark box, there was no difference in time spent in light vs . dark zones, latency to enter the light zone, or transitions between zones. For panels a-g, RE n = 7 , NP n = 8 . h There was no difference in time spent immobile on the forced swim test. RE n = 8 , NP n = 10 . i-k Stratification of pup retrieval behaviors by estrous stage Fisher’s exact test, p = 0.1002 . RE n = 7 , NP n = 6 . l-n Stratification of contextual fear conditioning by estrous stage Fisher’s exact test, p = 0.6041 . RE n = 7 , NP n = 10 . Error bars represent mean ± SEM. All p -values from Fisher’s exact tests comparing estrous stage distributions across groups for all behavioral experiments are provided in Supplementary Table 19 /articles/s41586-026-10509-4 MOESM3 . Extended Data Fig. 3 Postpartum experiences modulate the extent of RE-associated dHF transcriptomic and behavioral adaptations. /articles/s41586-026-10509-4/figures/8 a Experimental design to examine the contribution of discrete female reproductive events. Created in BioRender; Lab, M. https://biorender.com/4tj0uq7 https://biorender.com/4tj0uq7 2026 . b Chi square analyses showing significant shifts towards pup-directed behaviors after four days of pup exposure in Pup Sensitized virgin females χ² 7 = 33.52, p = 2.11 × 10−5 . c Unsupervised principle components analysis of NP vs . RE DEGs DESeq2, padj < 0.05 from dHF. d Heatmap of differential expression profiles of the top 500 genes ranked by ascending nominal p -values, DESeq2 between NP vs . RE. e Heatmaps of DEGs for NP vs . Mating + No Pregnancy, f NP vs . Mating + Pregnancy + Birth, g NP vs . Pup sensitized dHF transcriptomes. h Odds ratio analysis of DEG overlap for all comparisons vs . NP . Insert numbers indicate respective p values for each association following Fisher’s exact test N.S., p 0.05 . i Overlap of DEGs vs . NP, padj < 0.05 for each gene co-expression module Fisher’s exact test, p < 0.05, p < 0.01, p < 0.001, p < 0.0001 . High RE sensitivity modules are in bold. j Threshold-free comparison of indicated comparisons by rank-rank hypergeometric overlap. The lower left and upper right quadrants represent concordant gene regulation. k Shared ontologies following individual reproductive exposures vs . processes exclusive to the combined experiences in RE FDR < 0.05 . Extended Data Fig. 4 RE adaptations in the dHF transcriptome integrate pregnancy and postpartum experiences. /articles/s41586-026-10509-4/figures/9 a Experimental design to examine the trajectory of RE dHF gene expression changes across pregnancy and postpartum periods dpc: days post-conception; dpp: days postpartum . b Enrichment of DEGs vs . NP with RE-sensitive gene co-expression modules across reproductive time points in high left, dHF , moderate middle, vHF , and low right, mPFC transcriptional sensitivity to RE status Fisher’s exact test, p < 0.05, p < 0.01, p < 0.001, p < 0.0001 . c Time course analysis of dHF gene expression identifying patterns of transient vs . persistent up- and downregulation across pregnancy and postpartum FDR < 0.05 . d-g Select GO terms and KEGG pathways enriched for dHF genes exhibiting d persistent downregulation, e persistent upregulation, f transient downregulation, and g transient upregulation across reproductive stages FDR < 0.05 . Extended Data Fig. 5 Postpartum stress attenuates spatial learning enhancements independent of metabolic, locomotor, or hormonal measures. /articles/s41586-026-10509-4/figures/10 a Pup separation increases maternal corticosterone levels over 3-h two-way rmANOVA, group: p = 0.0006, Tukey’s multiple comparisons test: Control vs . Stress RE, p < 0.01 . n = 6 animals/group. b No effect of postpartum stress on litter size. n = 12/group. c Postpartum stress did not alter maternal weights. d Unsupervised principal components analysis of dHF transcriptomes. e All groups discriminated the novel location following 10-min of training on the object location task one way ANOVA, p = 0.9378 n = 6 animals/group. f Control RE n = 10 , Stress RE n = 11 , and NP n = 12 on the object location task following 5-min of training one-way ANOVA, p = 0.0278, Holm-Šídák multiple comparisons test; p < 0.05, p < 0.1 . g Stratification of object location test score by estrous stage Fisher’s exact test, p = 0.5551 . Control RE n = 10 , Stress RE n = 11 , NP n = 12 . h There was no difference in the total distance traveled in the open field test across groups. Control RE n = 12 , Stress RE n = 11 , NP n = 11 . Error bars represent mean ± SEM. Extended Data Fig. 6 Early postpartum stress impact on maternal dHF plasticity. /articles/s41586-026-10509-4/figures/11 a Early postpartum stress paradigm. Dams were provided with limited nesting and subjected to pup separation 3 h daily between 2–12 dpp. Created in BioRender; Lab, M. https://biorender.com/mvh9p5o https://biorender.com/mvh9p5o 2026 . b Unsupervised principal components analysis of dHF transcriptomes. c Heatmap of all significant DEGs DESeq2, padj < 0.05 with hierarchical clustering. d Venn diagrams showing overlap of DEGs altered in Control RE vs . NP and exhibiting opposing directionality in Early Stress RE vs . Control RE . e DEGs between Control RE vs . Early Stress RE left and Control RE vs . Late Stress RE right , shown as percentages of total DEGs. Slashes denote the fraction reversed in Control RE relative to NP. f RE-dependent enrichment in select Gene Ontology pathways Cellular Component are significantly reversed in Early Stress RE dHF FDR < 0.05 . g Control RE, but not Early Stress RE, dams froze more than NP during the acquisition phase of contextual fear conditioning two-way rmANOVA, interaction p = 0.0179, Tukey’s multiple comparisons test: Control RE vs . NP, p < 0.05 . h During the context recall test, NP froze significantly less than Control RE one-way ANOVA, p = 0.0219; Tukey’s multiple comparisons test, p < 0.05 . i Control RE dams spent more time investigating the moved object in the object location task compared to NP and Early Stress RE dams after 5-min training and a 24-h delay before testing one-way ANOVA, p = 0.001; Tukey’s multiple comparisons test: p < 0.01 For panels g-i, Control RE n = 11 , Early Stress RE n = 11 , NP n = 14 . Error bars represent mean ± SEM. Extended Data Fig. 7 Validation of snRNA-seq analyses. /articles/s41586-026-10509-4/figures/12 Violin plots showing a number of detected genes, b UMI number, and c mitochondrial gene percentage, by group and d-f annotated clusters. Unsupervised clustering of snRNA-seq data with cells colored for g neuronal, excitatory, inhibitory markers; h hippocampal subregional; and i non-neuronal markers. j-k Point-range plots depicting proportions of cell-type clusters between j NP vs . Control RE, and k Stress RE vs . Control RE scProportionTest, FDR < 0.05 and |log2Differerence| log2 1.5 . n = 6 NP, 4 Control RE, 5 Stress RE. l Number of DEGs following pseudobulk analysis p < 0.05 and FC 1.5 . m-n Ontology analyses for m Stress vs . Control RE and n NP vs . Control RE pseudobulk DEGs FDR < 0.05 . o Overlap of bulk and pseudobulk DEGs by directionality. p-q KEGG pathways enriched from overlapping DEGs from o, p < 0.05 . r UMAPs of Drd1 and Drd2 expression. s Quantification of Drd1 and Drd2 puncta in CA1 and DG three-way LMM; group x region x probe p = 2.466 × 10−9 . Box plots show median center line , interquartile range box , and minimum-to-maximum values whiskers . t-u Chi-square tests of Drd1 and Drd2 puncta distributions across expression bins for t CA1 χ² 4 = 30.69, p = 3.5 × 10−6 and u DG χ² 4 = 63.12, p = 6.4 × 10−13 . v Basal dopamine in NAc unpaired two-sided Student’s t-test, p = 0.0371 . Acute n = 5 , Chronic n = 4 . w Dopamine levels in VTA n : t0 = 6; t30 = 8; t180 = 7; t210 = 6 and LC n : t0 = 5; t30 = 8; t180 = 8; t210 = 6 during acute pup separation two-way ANOVA followed by Dunnett’s multiple comparisons test: VTA t0 vs . VTA t180, p = 0.0476 . Error bars represent mean ± SEM. Extended Data Fig. 8 H3K4me3Q5dop corresponds with transcription. /articles/s41586-026-10509-4/figures/13 a-b Validation of TG2 knockout cells using a western blot, and b transamidation activity assay. c Heatmaps showing H3K4me3 Q5dop enrichment in TG2 KO vs . WT. d Heatmaps of H3K4me3Q5dop peaks vs . IgG at all enriched loci in NP dHF. e, n Distribution of H3K4me3Q5dop across genic loci. f, o Overlap of H3K4me3Q5dop at promoter-TSS and gene body regions. g, p H3K4me3Q5dop enrichment within 2kB of TSS. h, q Boxplots of log10-transformed gene expression across quartiles of H3K4me3Q5dop signal for h NP top, n = 4 , Control RE middle, n = 4 , Stress RE bottom, n = 4 , and q NP left, n = 3 , parous subjects right, n = 5 . Differences were assessed within each group Kruskal-Wallis test, p < 2.2 × 10−16 . Boxes indicate interquartile range with median center line, and whiskers extend to 1.5x IQR. i, r Two-sided Spearman’s rank correlations between H3K4me3Q5dop enrichment vs . mean gene expression log10-transformed for i NP top; ρ = 0.55 , Control RE middle; ρ = 0.52 , Stress RE bottom; ρ = 0.56 and r NP left; ρ = 0.41 , Parous right; ρ = 0.40 . For all, p < 2.2 × 10−16. j Heatmaps of differential peaks Diffbind, p < 0.05; log2FoldChange ≥ |0.1| , by directionality. k, s Odds ratio analysis of differential H3K4me3Q5dop peaks and DEGs Fisher’s exact test; inset p -values; N.S. p 0.05 . l Statistical comparison of H3K4me3Q5dop profiles across DEG-associated loci showing –log10 Bonferroni-corrected p -values from bin-wise two-sided Wilcoxon rank-sum tests. Gray bars indicate non-significant bins after multiple-testing correction. Orange lines indicate statistical significance of signal differences across the TSS window, with higher values indicating greater confidence in differential enrichment between groups. Mouse tissues are shown in d-l and human tissues in m-s. For gel source data, see Supplementary Fig. 1 /articles/s41586-026-10509-4 MOESM1 . Extended Data Fig. 9 Assessment of VTA-dHF projection inhibition. /articles/s41586-026-10509-4/figures/14 a Representative images scale bar, 20 μm and b quantification of c-Fos in VTA TH+ neurons unpaired two-sided Student’s t-test; p = 0.03 . n = 3 animals/group. c NP-mCherry n = 14 , NP-hM4Di n = 15 , RE-mCherry n = 12 , and RE-hM4Di n = 11 weights post-DCZ injections two-way ANOVA; Tukey’s multiple comparisons corrections; p < 0.0001 . d NP-mCherry n = 11 , NP-hM4Di n = 12 , RE-mCherry n = 10 , and RE-hM4Di n = 9 during context recall two-way ANOVA , p = 0.0109; Fisher’s LSD, p = 0.05 . e-f NP-mCherry and RE-mCherry binned by high vs . low estrogen E during e pup retrieval two-way ANOVA, group: p < 0.0001, E: p = 0.8699 ; Tukey’s multiple comparisons test, p = 0.0001; RE-mCherry Low-E n = 4 , RE-mCherry High-E n = 8 , NP-mCherry Low-E n = 5 , NP-mCherry High-E n = 8 , and f contextual fear acquisition three-way rmANOVA, group: p = 0.0208, preshock: p < 0.0001, E: p = 0.3199; RE-mCherry Low-E n = 3 , RE-mCherry High-E n = 7 , NP-mCherry Low-E n = 4 , NP-mCherry High-E n = 6 . Error bars represent mean ± SEM. g Drd1 and Drd2 puncta quantification three-way LMM; group x region: p = 2.559 × 10−9 . Box plots show median center line , interquartile range box , minimum-to-maximum values whiskers . n = 3 animals/group. h-i Chi-square tests of Drd1 and Drd2 puncta distributions across expression bins for h CA1 χ2 2 = 4.33, p = 0.1145 and i DG χ² 2 = 48.67, p = 2.7 × 10−11 . j Heatmap of top 500 genes ranked by NP-mCherry vs . RE-mCherry ascending p -values . k DEG overlap with RE-sensitive modules p < 0.0001 . l-m H3K4me3Q5dop profiles across DEG-associated loci, showing –log10 Bonferroni-corrected p -values from bin-wise two-sided Wilcoxon rank-sum tests across the TSS window. Extended Data Fig. 10 Effect of viral H3.3Q5A transduction in mouse dHF. /articles/s41586-026-10509-4/figures/15 a Heatmaps of all differential peaks Diffbind, p < 0.05; log2FoldChange ≥ |0.1| between Stress RE H3.3WT vs . Stress RE H3.3Q5A, separated by directionality. b Overlap of DEGs with gene co-expression modules Fisher’s exact test; p < 0.01, p < 0.1 . High sensitivity RE modules are indicated in bold. c GO term enrichment for overlapping differential H3 dopaminylation peaks and genes in Control RE H3.3WT vs . Stress RE H3.3WT purple and Stress RE H3.3WT vs . Stress RE H3.3Q5A yellow comparisons FDR < 0.05 . d Normalized counts of select transcription factors and chromatin remodelers predicted as upstream regulators of transcriptional changes p < 0.05 between Stress RE H3.3WT vs . Stress RE H3.3Q5A. Differences between groups were assessed using multiple two-sided unpaired Student’s t-tests with FDR correction Benjamini-Krieger-Yekutieli; q < 0.05 . n = 5/group. Error bars represent median ± interquartile range. e Schematic of working model. Panel e created in BioRender; Lab, M. https://biorender.com/4f8zyrj https://biorender.com/4f8zyrj 2026 . Supplementary information Supplementary Information download PDF https://static-content.springer.com/esm/art%3A10.1038%2Fs41586-026-10509-4/MediaObjects/41586 2026 10509 MOESM1 ESM.pdf This file contains Supplementary Figs. 1–3 and the legends for Supplementary Tables 1–19. Supplementary Tables download XLSX https://static-content.springer.com/esm/art%3A10.1038%2Fs41586-026-10509-4/MediaObjects/41586 2026 10509 MOESM3 ESM.xlsx This file contains Supplementary Tables 1–19. Rights and permissions Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author s and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/ http://creativecommons.org/licenses/by/4.0/ . About this article Cite this article O’Chan, J.C., Di Salvo, G., Cunningham, A.M. et al. Dopamine drives persistent remodelling of the maternal brain. Nature 2026 . https://doi.org/10.1038/s41586-026-10509-4 Received: Accepted: Published: Version of record: DOI: https://doi.org/10.1038/s41586-026-10509-4