Human participants
Brain tissues used in this study were provided by the Douglas Brain Bank (RRID: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 individuals51,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 described51. 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 atlas53 (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.
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 described54,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 described16,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 behaviours58,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 described60.
Observation of pup-directed behaviours in the home cage
Confirmation of maternal behaviour in pup sensitized females was conducted based on prior studies61. 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 exposures62,63. 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.
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 EstrousNet64, 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. 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 described31,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 effects65—at 1 μg kg−1 in 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 μm2 nuclear 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.
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.
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 previously68,69. Differential expression analysis was performed using DESeq2 (v1.38.3)70, with significant genes defined by an adjusted 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)71 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 dataset72,73. The first two principle components (PC1 and PC2) were incorporated into the DESeq2 design formula, with the number of RUVSeq factors (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, Revigo76 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 lists77. 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 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 log10-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 log2(fold change) between groups, normalized values per cell type per group were averaged and expressed as the log2-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 protocol80. In brief, thawed tissues were placed in 1 ml of lysis buffer (0.32 M sucrose, 5 mM CaCl2, 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−1 RNase 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−1 RNase 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,600g 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.081, 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) package82 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)83 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 genes84,85,86,87,88,89,90, and label transfer from hippocampal reference datasets, including the Allen Brain Map and Broad Institute resources89,91. Clusters with contaminant cell populations expressing markers for choroid plexus (Ttr), ependymal (Tmem212), and vascular leptomeningeal cells (Vtn and Col1a2) were removed from the analysis90,92,93,94. 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: Rorb and Foxp294; ventral granule neurons: Tox386) were also removed from the analysis. Cell cluster proportion analyses were conducted using the scProportionTest package95, 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 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: log2(fold change) > 1.5 and P < 0.05.
CUT&RUN-seq
The procedure was adapted from established protocols33,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,100g 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 CaCl2, 1 mM MnCl2). 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 CaCl2, 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−1 RNase 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 3098. 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 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 number99. 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.05100. Peak annotation was conducted using HOMER (v4.1.1)101. Heat maps were made either using the DiffBind (v3.8.4) or deepTools (v3.5.5) packages102. 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 package103 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 log2(fold change) threshold > 0.1 and defined at p < 0.05, where log2(fold change) was calculated as log2(parity) − log2(NP), based on prior empirical observations used to define thresholds for differential peaks36. ChEA analysis on annotated loci was conducted using EnrichR with a significance threshold of adjusted P < 0.05 (ref. 104).
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
