Introduction
Epigenetic modifications of DNA and chromatin-associated histone proteins establish and maintain the unique patterns of gene expression in maturing and adult neurons (Kundakovic and Champagne, 2015). Neuron development requires the reconfiguration of epigenetic modifications, including methylation of genomic cytosine (DNA methylation, or mC) (Guo et al., 2014; Lister et al., 2013; Stroud et al., 2017) as well as covalent histone modifications associated with active or repressed gene transcription (Fagiolini et al., 2009; Putignano et al., 2007). While mC primarily occurs at CG dinucleotides (mCG) in mammalian tissues, neurons also accumulate a substantial amount of non-CG methylation (mCH) during postnatal brain development in the first 2–3 weeks of life in mice and the first two decades in humans (Lister et al., 2013). Accumulation of mCH, and the gain of mCG at specific sites, depend on the activity of the de novo DNA methyltransferase DNMT3A (Gabel et al., 2015). In mice, the abundance of
One challenge in investigating the developmental role of
To address the role of
Results
Previous studies of
Figure 1.
(A) An experimental model of the conditional loss of
Figure 1—figure supplement 1.
(A) Expression quantification of in situ hybridization data of gene
Figure 1—figure supplement 2.
(A–C) Cell-type-specific expression of nuclear membrane tag (Sun1-sfGFP-myc) in non-inhibitory NEUROD6+ neurons in mouse brain. (A) Overview (10× magnification) of a sagittal section stained with anti-GFP antibody. Note that the Sun1-sfGFP nuclear tag expression is region-specific, as described for the characterization of the NexCre mouse line (Goebbels et al., 2006). The numbered insets in (A) are enlarged and shown in (B), as example regions where the nuclear membrane tag is expressed (mPFC and CA1 hippocampal region) and where it is not (caudate putamen and thalamus). (C) Confocal images (60× magnification, maximum intensity projection) of mouse brain slices (40 µm thick) triple-labeled with antibodies anti-GFP (Sun1-sfGFP-tag, green channel), anti-GAD67 (inhibitory cells marker, red channel), and anti-NeuN (neuronal marker, gray). Note that Sun1-sfGFP-tag is only expressed in neurons that are not inhibitory. (D) Expression of pan-neuronal, excitatory, and inhibitory neuron marker genes in our RNA-seq data. TPM, transcripts per million.
Figure 1—figure supplement 3.
(A) Genome browser tracks of mRNA-seq data show confirmation of the deletion of
Figure 1—figure supplement 4.
The conditional ablation of
(A) and (D)
Figure 1—figure supplement 5.
(A) The increased PPI accompanied the impairment in startle responses to a 120 dB tone played at three time points during the recording session (HAB1 – beginning of the session; HAB2 – middle of the session; HAB3 – end of the session) (Wilcoxon test p=0.0027, 0.0019, and 0.0035 in male HAB1, HAB2, HAB3, respectively, and not significant in female). The habituation to the 120 dB auditory tone (i.e. the relative reduction in startle response throughout the experiment) was not significantly different between genotypes. (B) The percentage of PPI at prepulse intensity of 69, 73, and 81 dB (4, 8, and 16 dB above the 65 dB background, respectively) was increased in male, but not female mice (Wilcoxon test, ***, p=0.00076; *, p=0.016; n.s, not significant). (C–F) Fear learning and extinction were tested over 4 consecutive days. N=14–15 per group. (C) Fear acquisition to three tone-shock pairings occurred on day 1; (D) contextual fear in relation to the acquisition context (8 min, Block = 2 min) was measured on day 2; (E) cued fear recall and extinction training occurred on day 3 (Block = 4 tone trials); and (F) extinction recall (Block = 4 tone trials) occurred on day 4. Wilcoxon test reported no significant changes between the
We focused on cognitive domains associated with neurodevelopmental illness, including working memory and sensorimotor gating (Habib et al., 2019) and social interest (Dodell-Feder et al., 2015).
To test whether these deficits in specific neurocognitive domains reflect generalized impairment in brain function, we assessed long-term memory using a fear conditioning paradigm. There were no significant differences between
Loss of
To test the impact of
Figure 2.
Immature spine morphology and reduced excitability of layer 2 excitatory neurons following
(A) Example dendritic segments of layer 2 pyramidal neurons in the prelimbic region labeled with DiI and visualized using a 63× objective coupled to an Airyscan confocal microscope. Arrowheads show filopodia, which were more abundant in
Figure 2—figure supplement 1.
The membrane protrusions in the
Related to Figure 2C. Here, each line/box represents data from one dendrite fragment. (A) and (C) The cumulative distribution of the length (A) and width (C) of membrane protrusions in the
Figure 2—figure supplement 2.
The conditional ablation of
(A) Action potentials were initiated at a similar membrane potential in both genotypes (Wilcoxon test, n.s., not significant). (B) While miniature excitatory postsynaptic currents (mEPSCs) amplitude was not significantly different between genotypes (Wilcoxon test, n.s., not significant), it was slightly, yet significantly, more variable in the
To test how the loss of
Altered gene expression in
To investigate the impact of epigenetic disruption on gene expression, we compared the transcriptomes of cKO and control excitatory neuron nuclei in mature mice (postnatal day 39) (Supplementary file 2). We isolated nuclei from excitatory neurons in frontal cortex by backcrossing the
Figure 3.
Loss of
(A) 70 genes were differentially expressed (DE) (false discovery rate [FDR] < 0.05) in P39 pyramidal neurons in
Figure 3—figure supplement 1.
RNA-seq data showed transcriptomic disruption in P39
(A) Correlation matrix of gene expression (log10TPM) in the biological replicates of the control and
Figure 3—figure supplement 2.
The expression of transposable elements (TEs) was not affected by
TE subfamily expression was estimated in fragments per kilobase per million (FPKM) from P39 control and
Figure 3—figure supplement 3.
Genome-wide reduction of DNA methylation was observed in
(A) Heatmap to show the Spearman correlations across the control and
Figure 3—figure supplement 4.
Reduction of DNA methylation cannot fully explain the disruption in the transcriptome after
(A) Correlation of gene expression and gene body mCH level for up-regulated genes (red), down-regulated genes (blue), and non-differentially expressed (DE) genes (black, false discovery rate [FDR] ≥ 0.05 and fold-change <1.1) in the P39 control samples. For each gene group, genes are stratified by their expression in the control sample by 15 bins, and the mean gene body mCH levels are plotted. The shaded ribbon areas indicate the standard error of the mean. TPM, transcripts per million. (B) Violin plots to show the expression levels of non-DE genes (gray) selected to match baseline expression as those in significantly (FDR < 0.05) up-regulated genes (left, red) and down-regulated genes (right, blue). (C) CG DNA methylation (mCG) in P39 pyramidal cells in 1 kb bins in the region around the transcription start (TSS) and end site (TES) of DE genes and non-DE genes with matched expression levels. The lines denote the means across genes in each gene set, and the shared areas represent the 95% confidence intervals of the means. (D) Density scatter plots show the relationship between changes of gene body methylation (delta mCG or mCH) and the gene expression fold-changes for expressed genes (14,754 genes) between P39
Figure 3—figure supplement 5.
Differentially methylated regions (DMRs) in P39
(A) 2D distribution of mCG levels in P39 control vs.
Deletion of
Reduced DNA methylation does not fully explain altered transcription in
We investigated whether the altered gene expression in
We also examined the pattern of mCG at the promoter and gene body of DE genes. In contrast with the pattern of mCH, mCG was not significantly different between DE and control genes (Figure 3—figure supplement 4C).
The difference in gene body methylation (cKO – control) was negatively correlated with gene expression changes, consistent with repressive regulation (Gabel et al., 2015; Lavery et al., 2020; Figure 3E). This correlation accounted for 0.46% of the variance of differential gene expression, whereas the total explainable variance (R2 between biological replicates) was 1.30% (Figure 3—figure supplement 4D). The strength of the association between mCH and mRNA changes may be limited by the use of only two biological replicates in our dataset.
We next sought to determine if up- and down-regulated genes differ in ways that could explain their different responses to the loss of
In addition to promoters and gene bodies, distal regulatory elements such as enhancers are major sites of dynamic DNA methylation where epigenetic regulation can activate or repress the expression of genes over long genomic distances through 3D chromatin interactions (Malik et al., 2014). We investigated gene regulatory elements by identifying differentially methylated regions (DMRs) where mCG is altered in cKO compared to control neurons. We found a limited number of DMRs in newborn mice (P0: 1087 DMRs with lower, 164 with higher mCG in cKO; ≥30 difference in %mCG; FDR < 0.01). In mature neurons (P39), by contrast, we found 222,006 DMRs with substantially lower mCG in cKO compared with controls (Figure 3—figure supplement 5A; Supplementary file 3). Only 89 DMRs had ≥30 higher %mCG in cKO. To illustrate, we found five DMRs in an ~40 kb region around the promoter of the DE gene
The majority of P39 cKO DMRs (68.0%) were distal (≥10 kb) from the annotated transcription start sites. These DMRs were significantly enriched in both active enhancers and repressed chromatin, suggesting they have a regulatory role (Figure 3—figure supplement 5C; see also below). We found 113,557 developmental DMRs that gain mCG between P0 and P39 in control neurons (≥30 difference in %mCG, FDR < 0.01). These DMRs strongly overlapped (83.2%) with the cKO DMRs in mature neurons (Figure 3G–H, Supplementary file 3). Moreover, the P39 cKO DMRs were enriched in DNA sequence motifs of multiple transcription factors (TFs) associated with neuronal differentiation, including
Increased PRC2-associated repressive histone modification H3K27me3 in
Given that the cKO loses DNA methylation throughout the genome, we were surprised that the expression level of many genes was not disrupted. We, therefore, explored other potential epigenetic regulators which could contribute to maintaining gene expression. To identify TFs and chromatin regulators with experimental evidence of binding at cis-regulatory regions of the DE genes, we performed Binding Analysis for Regulation of Transcription (BART) (Wang et al., 2018). Chromatin regulators associated with PRC2, including
Figure 4.
Polycomb repressive complex 2 (PRC2) associated histone modification H3K27me3 is up-regulated following the loss of DNA methylation.
(A) Transcription factors (TFs) predicted to regulate P39
Figure 4—figure supplement 1.
Increased signal of the repressive histone mark H3K27me3 after
(A) Transcription factors (TFs) associated with chromatin organization were predicted to regulate differentially expressed (DE) genes in the P39
Figure 4—figure supplement 2.
The increased H3K27me signal in P39 was generally not observed in E14 or P0.
(A–B) Correlations of the differences of H3K27me3 signal between
Figure 4—figure supplement 3.
H3K27me3 differentially modified (DM) regions in P39
(A) H3K27me3 signals in DM regions (left) compared to random shuffles across the genome (middle; keeping peak sizes and the chromosome same as DM regions) and random shuffles within H3K27me3 peaks (right; H3K27me3 peaks are the union set of cKO and control peaks; restrict shuffles within the said peaks and keeping peak sizes and the chromosome same as DM regions). ****, Wilcoxon rank-sum test p<0.0001; n.s., not significant. (B) Top 10 Gene Ontology terms of the Biological Process ontology enriched in the genes associated with regions with up-regulated H3K27me3 signal in the P39
To experimentally address this, we performed chromatin immunoprecipitation sequencing (ChIP-seq) in excitatory neurons at embryonic day 14 (E14) and postnatal days 0 and 39 to measure trimethylation of histone H3 lysine 27 (H3K27me3), a repressive mark whose deposition is catalyzed by PRC2 and is important for transcriptional silencing of developmental genes. In P39 neurons, we also measured two histone modifications associated with active chromatin: H3K4me3 (trimethylation of histone H3 lysine 4, associated with promoters) and H3K27ac (acetylation of histone H3 lysine 27, associated with active promoters and enhancers) (Heinz et al., 2015). For each mark, we performed sequencing on two independent samples, each of which used pooled tissue from two mice. The active and repressive marks had positive and negative correlations with mRNA expression, respectively (Figure 4—figure supplement 1B). In addition, we noted regions where increased H3K27me3 concurred with the decreased DNA methylation. For example, at the
Using a conservative strategy to call ChIP-seq peaks (Zang et al., 2009), we found that marks associated with transcriptional activity (H3K4me3 and H3K27ac) were largely conserved in the P39 cKO and control (Figure 4D). By contrast, we found 51.9% more H3K27me3 peaks in mature (P39) cKO than in control neurons (Figure 4D, Supplementary file 6). When we directly identified differentially modified (DM) regions, we found no DM for H3K4me3 and H3K27ac between the cKO and control in P39 neurons (Figure 4—figure supplement 1C). By contrast, we found 4040 regions with significantly increased H3K27me3 in the P39 cKO, covering ~31.05 MB of the genome (Figure 4—figure supplement 1D, FDR < 0.05, Supplementary file 7). Differential H3K27me3 appears late during brain maturation: only 3 DM regions were found in earlier development stages (E14 or P0; Figure 4—figure supplement 1D), and the signal differences of cKO and control did not show clear correlations between P39 and earlier stages (Figure 4—figure supplement 2). These DM regions have a medium but non-zero level of H3K27me3 in the P39 control (higher than random shuffles across the whole genome but lower than random shuffles within the peak regions, Figure 4—figure supplement 3A), and hence presumably fine-tunable after
The increase in H3K27me3 was closely associated with regions that lost CG DNA methylation in the
Figure 5.
Increased H3K27me3 correlates with the loss of postnatal DNA methylation.
(A) Most of the P39 H3K27me3 peaks (57.1% and 66.2% of control and conditional knockout [cKO] peaks), but only some of the P39 H3K27ac peaks (28.9% and 33.1%), overlap with P39
Figure 5—figure supplement 1.
Differentially methylated regions (DMRs) are particularly unique in showing increased H3K27me3 signal.
(A) Histone modification chromatin immunoprecipitation sequencing (ChIP-seq) signal around the center of DMRs. Two sets of controlled regions were selected as comparisons: The first shuffled set (dotted lines) is regions generated by shuffling the DMRs randomly across the same chromosome (excluding blacklist regions) as the observed and keeping the same region size. The second shuffled set (dashed lines) is more stringent shuffles that meet the same criteria as the first set but with an additional restriction that the shuffles must reside within the P39 H3K27me3 peak regions (the union of P39 control and P39 cKO H3K27me3 peaks). RPKM, reads per kilobase per million. (B) The differences of H3K27me3 signal between
Figure 5—figure supplement 2.
Relationships between gene expression, H3K27me3 signal and gene body mCH.
(A) Correlations of gene expression fold-change and gene body H3K27me3 signal difference in P39
We further examined whether the increased H3K27me3 could account for the reduced expression of some genes in the
We next analyzed how changes in H3K27me3 related to the loss of mCH. In all non-DE genes (FDR ≥ 0.05), we observed a larger increase of H3K27me3 signal in the gene body of genes that lost most mCH in the cKO (Figure 5—figure supplement 2B). Indeed, when we grouped all genes by the extent of the loss of mCH, we found that genes that lost the most mCH had a negative correlation between the changes in gene body H3K27me3 signal and gene expression fold-change. Such correlation was not observed in genes that did not lose mCH (Figure 5—figure supplement 2C). These results suggest that H3K27me3 may have a role in repressing expression specifically in regions that lose repressive non-CG DNA methylation.
Developmental changes in H3K27me3 are not affected by
Our finding that
Figure 6.
Developmental dynamics of H3K27me3.
(A) Heatmap of developmentally regulated H3K27me3 regions in E14 and P39 control samples. CPM – counts per million; R1/2 – replicates. (B) Bar plots show the numbers of developmental differentially modified H3K27me3 regions (E14 vs. P39) that overlap developmental differentially methylated regions (DMRs) (P0 vs. P39, left panel), and the numbers of developmental DMRs that overlap developmental differentially modified H3K27me3 regions (right panel). (C–E) Normalized H3K27me3 signal (fold-changes compared to P39 control), and mCG, mCH differences (compared to the average of the two replicates from P39 control) in peaks that overlap with E14 vs. P39 developmental loss-of-H3K27me3 regions (C), developmental gain-of-H3K27me3 regions (D), or increased H3K27me3 in P39
Figure 6—figure supplement 1.
Regions prone to alteration of H3K27me3 by
(A) Differentially modified H3K27me3 regions during development in control samples. (B) Top 10 most enriched biological process terms from the Gene Ontology in the genes associated with differentially modified H3K27me3 regions during development in control samples. Gene-region association was predicted by GREAT (McLean et al., 2010). FDR, false discovery rate; the vertical dashed line shows the threshold of FDR = 0.05. (C) Correlations of P39 H3K27me3 signal fold-changes of
Figure 6—figure supplement 2.
Correlations of chromatin immunoprecipitation sequencing (ChIP-seq) signals across replicates.
Correlations were computed with ChIP-seq counts in 1 kb genomic bins (A) and in peaks (B) across replicates. r, Spearman correlation coefficient. n, number of peaks. RPKM, reads per kilobase million.
The P39
To further stratify the joint distribution of developmental and
Novel DNA methylation valleys with increased H3K27me3 signal in the
DNA methylation and H2K27me3 have complementary roles at DNA methylation valleys (DMVs), that is, large regions (≥5 kb) with low mCG (≤15%) that occur around key transcriptional regulators of development in human and mouse tissues (Mo et al., 2015; Xie et al., 2013). Previous studies comparing the epigenetic profile of DMVs across tissues identified multiple categories, including constitutive DMVs present in all tissues as well as tissue-specific DMVs (Li et al., 2018). We found more than twice as many DMVs in P39 cKO (1838) compared with control (881) neurons (Supplementary file 10), covering a greater genomic territory (16.02 Mbp in cKO, 7.91 Mbp in control). The majority of these P39 DMVs were either expanded or unique in cKO (Figure 7A), while DMVs identified in P0 samples were mostly not altered (Figure 7—figure supplement 1). Most P39 DMVs had active histone marks (H3K27ac+, H3K27me3-), while some had repressed or bivalent profiles consistent with PRC2-associated gene silencing (H3K27me3+) (Figure 7B).
Figure 7.
Distinct clusters of DNA methylation valleys (DMVs) were associated with the increased H3K27me3 signal in the
(A) Number of DMVs identified in the P39
Figure 7—figure supplement 1.
Number of DNA methylation valleys (DMVs) identified in the P0
P0 DMVs were categorized by whether they appear in one or both groups or change size in the cKO, as in Figure 7A.
By clustering the DMVs using their pattern of DNA methylation, chromatin modifications, and gene expression, we found eight distinct categories (Figure 7C–D). Whereas most DMVs lack H3K27me3 (clusters C2,3,5,7,8), we found some groups of DMVs associated with moderate (C1, C4) or high (C6) levels of H3K27me3. Cluster C6 DMVs, such as the promoter of
The remaining clusters lack H3K27me3 and are instead marked by low mCG and either high H3K27ac (cluster C8) or both H3K27ac and H3K4me3 (clusters C7). Cluster C5 DMVs (e.g.
Discussion
To examine the role of DNA methylation in cortical excitatory neurons after their birth, we developed a mouse model where the loss of DNA methylation occurs in postmitotic neurons, prior to the postnatal increase in
There is a strong antagonistic relationship between DNA methylation and the PRC2-associated histone mark, H3K27me3 (Brinkman et al., 2012; Jermann et al., 2014; Lynch et al., 2012; Reddington et al., 2013; Wu et al., 2010). Switching between polycomb- and DNA methylation-mediated repression has been observed during development and in cancer (Mohn et al., 2008; Schlesinger et al., 2007; Widschwendter et al., 2007). Severe depletion of mCG can lead to redistribution of H3K27me3, causing derepression of developmental regulators such as the
Previous work showed that the early embryonic deletion of
Our findings highlight the critical and interconnected roles in brain development and cognitive function of two major modes of epigenetic repression of gene expression: DNA methylation and PRC2-mediated repression. The loss of DNA methylation in excitatory neurons has effects on gene expression, synaptic function, and cognitive behavior. Moreover, loss of DNA methylation leads to a gain of the PRC2-associated repressive mark H3K27me3. Our cKO is a restricted manipulation of one neuron type, yet it directly impacts DNA methylation throughout the genome at millions of sites. PRC2-mediated repression may compensate for the loss of mCG and/or mCH, acting as an alternative repressive mechanism when DNA methylation is disrupted. Future work focusing on earlier developmental stages, and using targeted methods to manipulate epigenetic marks in local genomic regions (Liu et al., 2016), may help elucidate the causal interactions among epigenetic modifications that are critical for neuronal maturation and function.
Materials and methods
Generation of the
All animal procedures were conducted in accordance with the guidelines of the American Association for the Accreditation of Laboratory Animal Care and were approved by the Salk Institute for Biological Studies Institutional Animal Care and Use Committee (protocol number 18-00006). For behavior, slice physiology, and spine analyses,
Frontal cortex dissection, nuclei isolation, and flow cytometry
Frontal cortex tissue was produced as described (Lister et al., 2013; Luo et al., 2017) from postnatal day 0 and 39 (P39)
Western blot
Frontal cortex proteins were obtained by homogenization in RIPA buffer of the following composition: 150 mM NaCl, 10 mM Na2HPO4, 1% NaDOC, 1% NP-40, 0.5% SDS, 1 mM DTT, 1 mM PMSF in DMSO, supplemented with protease inhibitor (Sigma-Aldrich #11836153001) and phosphatase inhibitor (Pierce #A32957) cocktails. After centrifugation at 15,000×
Patch-clamp electrophysiology
Male and female mice (6–9 weeks) were anesthetized with isoflurane and decapitated. The brains were quickly removed and coronal slices of the frontal cortex containing the prelimbic region (~2 mm anterior to Bregma) were cut in an ice-cold slicing medium of the following composition (in mM): 110 sucrose, 2.5 KCl, 0.5 CaCl2, 7 MgCl2, 25 NaHCO3, 1.25 NaH2PO4, and 10 glucose (bubbled with 95% O2 and 5% CO2). The slices were then transferred to artificial CSF (aCSF) containing (in mM): 130 NaCl, 2.5 KCl, 1.25 NaH2PO4, 23 NaHCO3, 1.3 MgCl2, 2 CaCl2, and 10 glucose, equilibrated with 95% O2 and 5% CO2 at 35°C for 30 min and afterward maintained at room temperature (22–24°C) for at least 1 hr (patch-clamp recording) before use. Brain slices were then transferred to a recording chamber and kept minimally submerged under continuous superfusion with aCSF at a flow rate of ~2 ml/min. Whole-cell recordings were obtained from putative prelimbic layer 2 (L2) pyramidal cells (identified by their pyramidal-shaped cell bodies and long apical dendrite using an upright microscope equipped with differential interference contrast optics). In acute mPFC slices, the prelimbic L2 is clearly distinguishable from L1 and L3 as a thin dark band that is densely packed with neuron somata. Pipettes had a tip resistance of 4–8 MΩ when filled with an internal solution of the following composition (in mM): 125 K-gluconate, 15 KCl, 8 NaCl, 10 HEPES, 2 EGTA, 10 Na2 phosphocreatine, 4 MgATP, 0.3 NaGTP (pH 7.25 adjusted with KOH, 290–300 mOsm). Access resistance (typically 15–35 MΩ) was monitored throughout the experiment to ensure stable recordings.
After obtaining the whole-cell configuration in voltage-clamp mode, cells were switched from a holding potential of –70 mV to current-clamp mode and the bridge-balance adjustment was performed. Passive electrical properties were quantified from recordings with hyperpolarizing current injections that evoked small ~5 mV deflections in membrane potential from resting. Responses to stepwise current injections (10–300 pA in increments of 10 pA; duration, 1 s) were recorded at 20 kHz in order to calculate input-output curves and rheobase – the minimal current necessary to trigger the first action potential. Miniature excitatory postsynaptic currents (mEPSCs) were recorded for 5 min in voltage-clamp mode (Vh=−70 mV) in the presence of the Na+ channel blocker, TTX (0.5 μM), to prevent the generation of action potentials, and picrotoxin (50 μM), an antagonist of GABAA receptors, to minimize inhibitory responses. In these conditions, mEPSCs could be blocked by the AMPA receptor antagonist, CNQX (25 μM). Single events larger than 6 pA were detected offline using the Minianalysis program (Synaptosoft Inc Decatur, GA). All data were acquired using a Multiclamp 700B amplifier and pCLAMP 9 software (Molecular Devices, LLC, San Jose, CA).
Fluorescent labeling of dendritic spines
Coronal brain slices containing the mPFC of 10–13 weeks’ female mice were prepared as for electrophysiological recordings and placed in a beaker for 3 hr at room temperature (24°C) to allow functional and morphological recovery. One slice was then transferred to a recording chamber and kept minimally submerged under continuous superfusion with aCSF bubbled with carbogen (95% O2/5% CO2) at a flow rate of ~2 ml/min. Previously sonicated crystals of the fluorescent marker DiI were placed next to the somata of layer 2 neurons in the prelimbic cortex, identified with the aid of an upright microscope equipped with differential interference contrast optics. The mice in these experiments had a Thy1-YFP background to help rule out non-specific labeling of deeper layer neurons. The neurons were exposed to the DiI crystals for 60 min. The slices were then gently removed from the incubation chamber with a transfer pipette and immersed in fixative (4% PFA) for 30 min. Then, the slices were rinsed three times with PBS for 5–10 min each, after which they were mounted on slides with prolonged gold antifade mounting medium (Life Technologies – Molecular Probes). The slides were kept in a dark box for 24 hr at room temperature to allow the liquid medium to form a semi-rigid gel. Imaging took place 24–48 hr from the time of the initial staining.
Confocal imaging
Dendritic spines were imaged by an investigator blind to the genotype using a Zeiss AiryScan confocal laser scanning microscope. All images were taken using the Zeiss Plan-APOCHROMAT 63× oil-immersion lens (N/A 1.4). A 543 nm laser was used to visualize the fluorescence emitted by DiI. Serial stack images with a 0.2 μm step size were collected, and then projected to reconstruct a three-dimensional image that was post-processed by the AiryScan software. Dendritic segments in layer 1, which were derived from layer 2 pyramidal neurons retrogradely labeled with DiI and that were well separated from neighboring neural processes, were randomly sampled and imaged. Each dendritic segment imaged for quantification belonged to a different neuron.
Dendritic spine quantification
The z-stack series were imported into the Reconstruct software (https://synapseweb.clm.utexas.edu/software-0/), with which a second investigator also blind to the genotype performed the identification of dendritic spines and their morphometric analysis. By scrolling through the stack of different optical sections, individual spine heads could be identified with greater certainty. All dendritic protrusions with a clearly recognizable stalk were counted as spines. Spine density was determined by summing the total number of spines per dendritic segment length (30–40 μm) and then calculating the average number of spines per μm. Individual dendritic spines were classified in the following order according to pre-established criteria: protrusion longer than 3 μm, filopodia; head wider than 0.6 μm, mushroom; protrusion longer than 2 μm and head narrower than 0.6 μm, long-thin; protrusion longer than 1 μm and head narrower than 0.6 μm, long-thin; the remaining spines were labeled stubby. Branched spines (with more than one neck) were counted separately.
Behavioral testing
Phenotypic characterization was initiated when the animals reached 9 weeks of age using cohorts of 10–15 male or female mice per genotype, according to the order described below.
Open field test
The open field test was performed using MED Associates hardware and the Activity Monitor software according to the manufacturer’s instructions (MED Associates Inc, St Albans, VT). Animals were individually placed into clear Plexiglas boxes (43.38 × 43.38 × 30.28 cm3) surrounded by multiple bands of photo beams and optical sensors that measure horizontal and vertical activity. Movement was detected as breaks within the beam matrices and automatically recorded for 60 min.
Light/dark transfer test
The light/dark transfer procedure was used to assess anxiety-like behavior in mice by capitalizing on the conflict between exploration of a novel environment and the avoidance of a brightly lit open field (150–200 lux in our experiments). The apparatus were Plexiglas boxes as for the open field test (43.38 × 43.38 × 30.28 cm3) containing dark box inserts (43.38 × 12.8 × 30.28 cm3). The compartments were connected by an opening (5.00 × 5.00 cm2) located at floor level in the center of the partition. The time spent in the light compartment was used as a predictor of anxiety-like behavior, that is, a greater amount of time in the light compartment was indicative of decreased anxiety-like behavior. Mice were placed in the dark compartment (4–7 lux) at the beginning of the 15 min test.
Elevated plus maze
The maze consisted of four arms (two open without walls and two with enclosed walls) 30 cm long and 5 cm wide in the shape of a plus sign. The apparatus was elevated approximately 33 cm over a table. At the beginning of each trial, one animal was placed inside a cylinder located at the center of the maze for 1 min. The mouse was then allowed to explore the maze for 5 min. The session was video-recorded by an overhead camera and subjected to automated analysis using ANY-maze software. The apparatus was wiped down with sani-wipes between trials to remove traces of odor cues. The percentage of time spent in open or closed arms was scored and used for analysis.
Y-maze test for spontaneous alternations
Spontaneous alternations between three 38 cm long arms of a Y-maze were taken as a measure of working memory. Single 6 min trials were initiated by placing each mouse in the center of the Y-maze. Arm entries were recorded with a video camera and the total number of arm entries, as well as the order of entries, was determined. The apparatus was wiped down with sani-wipes between trials to remove traces of odor cues. Spontaneous alternations were defined as consecutive triplets of different arm choices and % spontaneous alternation was defined as the number of spontaneous alternations divided by the total number of arm entries minus 2.
Social approach
The apparatus consisted of a Plexiglas box (60 × 38 × 23.5 cm3) divided into three compartments by Plexiglas partitions containing openings through which the mice could freely enter the three chambers. The test was conducted in two 10 min phases. In phase I, the test mouse is first allowed to explore the chambers for 10 min. Each of the two outer chambers contained an empty, inverted stainless steel wire cup. In phase II, the test mouse is briefly removed, and a sex-matched unfamiliar mouse was placed under one of the wire cups, and plastic blocks were placed under the other wire cup. The test mouse was then gently placed back in the arena and given an additional 10 min to explore. An overhead camera and video tracking software (ANY-maze, Wood Dale, IL) were used to record the amount of time spent in each chamber. The location (left or right) of the novel object and novel mouse alternates across subjects.
Acoustic startle responses and PPI of the acoustic startle response
Acoustic startle responses were tested inside SR-LAB startle apparatus (San Diego Instruments, San Diego, CA), consisting of an inner chamber with a speaker mounted to the wall and a cylinder mounted on a piezoelectric sensing platform on the floor. At the beginning of testing, mice were placed inside the cylinder and then were subjected to background 65 dB white noise during a 5 min acclimation period. The PPI session began with the presentation of six pulse-alone trials of 120 dB, 40 ms. Then, a series of pulse-alone trials and prepulse trials (69, 73, or 81 dB, 20 ms followed by 100 ms pulse trial, 120 dB) were each presented 12 times in a pseudorandom order. The session concluded with the presentation of six pulse-alone trials. The apparatus was wiped down with sani-wipes between trials to remove traces of odor cues. The startle amplitude was calculated using arbitrary units, and the acoustic startle response was the average startle amplitude of pulse-alone trials. The percent PPI was calculated as follows: [100 − (mean prepulse response/mean pulse alone response) × 100].
Cued and contextual fear conditioning
Fear conditioning experiments were performed using automated fear conditioning chambers (San Diego Instruments, San Diego, CA), similar to previous studies (Gresack et al., 2010; Risbrough et al., 2014). On day 1, after a 2 min acclimation period, mice were presented with a tone conditioned stimulus (75 dB, 4 kHz) for 20 s that co-terminated with a foot shock unconditioned stimulus (1 s, 0.5 mA). A total of three tone-shock pairings were presented with an inter-trial interval of 40 s. To assess acquisition, freezing was quantified during foot shock presentations. Mice were returned to their home cages 2 min after the final shock. These moderate shock parameters were previously found suitable to detect both increases and decreases in fear-conditioned behavior (Risbrough et al., 2014). Twenty-four hr later, on day 2, mice were re-exposed to the conditioning chamber to assess context-dependent fear retention. This test lasted 8 min during which time no shocks or tones were presented and freezing was scored for the duration of the session. Time freezing was quantified across four 2 min blocks. Day 3: 24 hr after the context fear-retention test, mice were tested for CS-induced fear retention and extinction. The context of the chambers was altered across several dimensions (tactile, odor, visual) for this test in order to minimize generalization from the conditioning context. After a 2 min acclimation period, during which time no tones were presented (‘pre-tone'), 32 tones were presented for 20 s with an inter-trial interval of 5 s. Freezing was scored during each tone presentation and quantifications were done in eight blocks of four tones. Mice were returned to their home cage immediately after termination of the last tone. On day 4, after a 2 min acclimation period, during which time no tones were presented (‘pre-tone’), a shorter session of 16 tones was used to assess extinction. Time freezing was quantified across four blocks of four tones.
RNA extraction, RNA-seq library construction, and sequencing
Nuclei (between 50,000 and 60,000) were used to isolate RNA using Single-Cell RNA Purification Kit (Norgen, catalog# 51800). In brief, aliquots of nuclei were resuspended in 350 µl of RL buffer (Norgen) and passed through an 18 G syringe five times. RNA extraction, including DNase digestion, followed manufacturers’ instructions. RNA was eluted in 20 µl of Elution Solution A (Norgen). The nuclear RNA concentration was determined using TapeStation (Agilent). RNA was diluted to 1 ng/µl and a total of 5 ng was processed for RNA-seq library preparation. RNA libraries were prepared using NuGen Ovation RNA-Seq System V2 (#7102–32) for cDNA preparation following the product manual. cDNA purification was done using Zymo Research DNA Clean & Concentrator-25 with modification from the Ovation protocol. cDNA, eluted in 30–40 µl of TE (1 µg per sample), was fragmented at 300 bp using Covaris S2 (Sonolab S-series V2), followed by library preparation according to KAPA LTP Library Preparation Kit (KK8232), using Illumina indexed adapters. Libraries were sequenced on NovaSeq 6000.
DNA extraction
DNA extraction was performed using the Qiagen DNeasy Blood and Tissue kit (catalog #69504) and eluted into 50–100 µl AE.
Genomic DNA library construction and sequencing
1.5 µg of genomic DNA was fragmented with a Covaris S2 (Covaris, Woburn, MA) to 400 bp, followed by end repair and addition of a 3’ A base. Cytosine-methylated adapters provided by Illumina (Illumina, San Diego, CA) were ligated to the sonicated DNA at 16°C for 16 hr with T4 DNA ligase (New England Biolabs). Adapter-ligated DNA was isolated by two rounds of purification with AMPure XP beads (Beckman Coulter Genomics, Danvers, MA). Half of the adapter-ligated DNA molecules were enriched by 6 cycles of PCR with the following reaction composition: 25 µl of Kapa HiFi Hotstart Readymix (Kapa Biosystems, Woburn, MA) and 5 µl TruSeq PCR Primer Mix (Illumina) (50 µl final). The thermocycling parameters were: 95°C 2 min, 98°C 30 s, then 6 cycles of 98°C 15 s, 60°C 30 s, and 72°C 4 min, ending with one 72°C 10 min step. The reaction products were purified using AMPure XP beads and size selection was done from 400 to 600 bp. Libraries were sequenced on NovaSeq 6000.
MethylC-seq library construction and sequencing
MethylC-seq libraries were prepared as previously described (Urich et al., 2015). All DNA obtained from the extraction was spiked with 0.5% unmethylated Lambda DNA. The DNA was fragmented with a Covaris S2 (Covaris, Woburn, MA) to 300 bp, followed by end repair and addition of a 3’ A base. Cytosine-methylated adapters provided by Illumina (San Diego, CA) were ligated to the sonicated DNA at 16°C for 16 hr with T4 DNA ligase (New England Biolabs). Adapter-ligated DNA was isolated by two rounds of purification with AMPure XP beads (Beckman Coulter Genomics, Danvers, MA). Adapter-ligated DNA (≤450 ng) was subjected to sodium bisulfite conversion using the EZ methylation Direct kit (Zymo, D5021) as per the manufacturer’s instructions. The bisulfite-converted, adapter-ligated DNA molecules were enriched by 8 cycles of PCR with the following reaction composition: 25 µl of Kapa HiFi Hotstart Uracil +Readymix (Kapa Biosystems, Woburn, MA) and 5 µl TruSeq PCR Primer Mix (Illumina) (50 µl final). The thermocycling parameters were: 95°C 2 min, 98°C 30 s, then 8 cycles of 98°C 15 s, 60°C 30 s and 72°C 4 min, ending with one 72°C 10 min step. The reaction products were purified using AMPure XP beads. Up to two separate PCR reactions were performed on subsets of the adapter-ligated, bisulfite-converted DNA, yielding up to two independent libraries from the same biological sample. MethylC-seq libraries were sequenced on NovaSeq 6000.
ChIP-seq library construction and sequencing
Sorted nuclei were crosslinked for 15 min in 1% formaldehyde solution and quenched afterward with glycine at a final concentration of 0.125 M. After crosslinking, nuclei were sonicated in Lysis buffer (50 mM Tris HCl pH 8, 20 mM EDTA, 1% SDS, 1× EDTAfree protease inhibitor cocktail). ChIP assays were conducted with antibodies against H3K27me3 (39156, Active Motif), H3K27ac (39133, Active Motif), and H3K4me3 (04-745, Millipore Sigma). Mouse IgG (015-000-003, Jackson ImmunoResearch) served as a negative control. H3K4me3 ChIP-seq assays were conducted with 100 K nuclei and 500 K nuclei were used for H3K27me3 and H3K27ac ChIP-seq assays. The respective antibodies and IgG were coupled for 4–6 hr to Protein G Dynabeads (50 µl, 10004D, Thermo Fisher Scientific). Equal amounts of sonicated chromatin were diluted with 9 volumes of Binding buffer (1% Triton X-100, 0.1% Sodium Deoxycholate, 2× EDTA free protease inhibitor cocktail) and subsequently incubated overnight with the respective antibody-coupled Protein G beads. Beads were washed successively with low salt buffer (50 mM Tris HCl pH 7.4, 150 mM NaCl, 2 mM EDTA, 0.5% Triton X-100), high salt buffer (50 mM Tris HCl pH 7.4, 500 mM NaCl, 2 mM EDTA, 0.5% Triton X-100) and wash buffer (50 mM Tris HCl pH 7.4, 50 mM NaCl, 2 mM EDTA) before de-crosslinking, proteinase K digestion, and DNA precipitation. Libraries were generated with the Accel-NGS 2S Plus DNA Library Kit (21024, Swift Biosciences) and sequenced on the Illumina HiSeq 4000 Sequencing system.
RNA-seq data processing
RNA-seq reads first went through quality control using FastQC (Andrews et al., 2012) (v0.11.8, https://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and then were trimmed to remove sequencing adapters and low-quality sequences (minimum Phred score 20) using Trim Galore (v0.5.0, https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/, a wrapper tool powered by Cutadapt [Martin, 2011] v1.16) in the paired-end mode. Clean reads were then mapped to the mouse mm10 (GRCm38) genome and the GENCODE annotated transcriptome (release M10) with STAR (Dobin et al., 2013) (Spliced Transcripts Alignment to a Reference, v2.5.1b). Gene expression was estimated using RSEM (Li and Dewey, 2011) (RNA-Seq by Expectation Maximization, v1.2.30). Gene-level ‘expected count’ from the RSEM results were rounded and fed into edgeR (Robinson et al., 2010) (v3.24.1) to call DE genes. Only genes that were expressed (with counts per million [CPM] > 2) in at least two samples were kept. These counts were then normalized using the TMM method (Robinson and Oshlack, 2010), and DE genes were then called in the quasi-likelihood F-test mode, requiring FDR < 0.05 and FC > 20% (Supplementary file 2). To quantify TE expression, we used the Fetch, Clean, Map, and Count modules from SQuIRE (Yang et al., 2019) with the RepeatMasker annotation from UCSC. To evaluate the differences in TE family and class abundance between
Enrichment test of GO terms in DE genes
Gene Ontology (GO) enrichment analysis was performed using clusterProfiler (Yu et al., 2012) (v3.10.0). Only ‘Biological Process’ terms with no less than 10 genes and no more than 250 genes were considered. Terms with FDR < 0.05 were considered significantly enriched.
Genomic DNA sequencing data processing and SNP calling
To estimate the completeness of the inbreeding of the mouse strains, and to avoid incorrect cytosine context assignment in the following MethylC-seq data processing, we used the genomic DNA sequencing data of both the
MethylC-seq data processing
MethylC-seq reads were processed using the methylpy pipeline (v1.3.2, https://github.com/yupenghe/methylpy; He, 2021) as previously described (Lister et al., 2013; Mo et al., 2015). Briefly, a computationally bisulfite-converted genome index was built using the aforementioned substituted genome file appended with the lambda phage genomic sequence. MethylC-seq raw reads were first trimmed to remove sequencing adapters and low-quality sequences (minimum Phred score 10) using Cutadapt in paired-end mode. To acquire higher mappability, we treated the two ends of the clean reads as they were sequenced in single-end mode, and mapped them to the converted genome index with bowtie2 (Langmead and Salzberg, 2012) (v2.3.0) as aligner in the single-end pipeline of methylpy. Only reads uniquely mapped were kept, and clonal reads were removed. The bisulfite non-conversion rate (NCR) was estimated using the spiked-in unmethylated lambda phage DNA. For each cytosine, a binomial test was performed to test whether the methylation levels are significantly greater than 0 with an FDR threshold of 0.01.
For a particular genomic region, the raw methylation level for a given cytosine context (CG or CH) was defined as:
where m is the total number of methylated based calls within the region, and h is the total number of covered based calls within the region. Methylation levels were then corrected for NCR using the following maximum likelihood formula:
When profiling the methylation landscapes around DE genes, we selected control genes with comparable gene expression using the R package MatchIt (Ho et al., 2011) (v3.0.2). These control genes were defined with nearest neighbor matching of the expression (in the unit of TPM) using logistic link propensity score as a distance measure, requiring the standard deviation of the distance to be less than 0.01.
DMRs calling
CG DMRs were identified using a previously reported method (Ma et al., 2014; Schmitz et al., 2013; Schultz et al., 2015), which is implemented in the DMRfind function in methylpy. We required at least three differentially methylated sites within a particular DMR, and significant sites located within 500 bp of each other were merged into the same DMR. With an FDR cutoff of 0.01 and a post-filtering cutoff of methylation levels change greater than 30 (in the unit of %mCG), we found 222,006
Enrichment test of DMRs and other genomic regions
To test whether DMRs were significantly enriched in certain genomic features, we use methods adapted from a recent report (Rizzardi et al., 2019). Briefly, for each genomic feature, we constructed a 2 by 2 contingency table of (n11, n12, n21, n22), where:
n11 is the number of CG sites in DMRs that were inside the feature;
n12 is the number of CG sites in DMRs that were outside of the feature;
n21 is the number of CG sites not in DMRs that were inside in feature;
n22 is the number of CG sites not in DMRs that were outside of the feature.
The total number of CG sites in consideration was the number of autosomal and chromosome X CG in the reference genome. Counting the number of CG rather than the number of DMRs or bases accounts for the non-uniform distribution of CG along the genome and avoids double-counting DMRs that are both inside and outside of the feature.
With this contingency table, we estimated the enrichment log odd ratio (OR) along with its standard error (se) and 95% confidence interval (ci) with the following formulas:
p-Value from performing Fisher’s exact test for testing the null of independence of rows and columns in the contingency table (the null of no enrichment or depletion) was computed using the fisher.test() function in R.
The genomic regions/features used in these enrichment tests include: a list of developmental DMRs that gain or lose methylation during development (Lister et al., 2013); gene features (genic, exonic, intronic, promoter, 5’UTR, 3’UTR, intergenic) based on the GENCODE vM10 annotation (promoters were defined as the ±2 kb regions around transcription start sites); CpG island (CGI)-related features based on the ‘cpgIslandExt’ annotation from the UCSC genome browser (Karolchik et al., 2004; Kent et al., 2002) (http://genome.ucsc.edu/index.html), where CGI shores were defined as CGI ± 2 kb, CGI shelves were defined as ±2–4 kb of CGI and open seas were defined as regions that were at least 4 kb away from any CGI; the 12 states of the chromatin states map in mouse embryonic stem cell (Pintacuda et al., 2017) (https://github.com/guifengwei/ChromHMM_mESC_mm10) generated by ChromHMM (Ernst and Kellis, 2017; Ernst and Kellis, 2012) using ChIP-seq data from the ENCODE project (Davis et al., 2018; ENCODE Project Consortium, 2012); the H3K4me3, H3K27ac, and H3K27me3 peaks and the H3K27me3 differentially binding regions generated with our ChIP-seq data (see the ‘ChIP-seq data processing’ section).
Predicting functional TFs regulating the DE genes with BART
The lists of up-regulated and down-regulated genes were fed into BART (Wang et al., 2016; Wang et al., 2018) (BART) separately to predict functional TFs and chromatin regulators that bind at cis-regulatory regions of the DE genes. To make a better visualization, we transformed the relative rank (a metric generated by BART to represent the average rank of Wilcoxon p-value, z-score, and max AUC for each factor divided by the total number of factors) into the functional TF rank score (which is simply 1 minus the relative rank) so that the higher the rank score the more possible the TF regulates the DE genes. The integrative rank significance was estimated with the Irwin-Hall p-value. TFs with Irwin-Hall p-value < 0.05 were considered significant.
Enrichment test of known TF binding motifs in DMRs
We used the ‘findMotifsGenome.pl’ tool in HOMER (Heinz et al., 2010) (Hypergeometric Optimization of Motif EnRichment, v4.8.3) to find known TF binding motif in the
ChIP-seq data processing
ChIP-seq reads were pre-processed with the ENCODE Transcription Factor and Histone ChIP-Seq processing pipeline (https://github.com/ENCODE-DCC/chip-seq-pipeline2, v1.1.6; Jin wook, 2022). Briefly, paired-end reads were mapped to the mm10 genome with BWA (Li and Durbin, 2009) (v0.7.13-r1126). Reads were then filtered using samtools (Li et al., 2009) (v1.2) to remove unmapped, mate unmapped, not primary alignment and duplicate reads (-F 1804). Properly paired reads were retained (-f 2). Multi-mapped reads (MAPQ < 30) were removed. PCR duplicates were removed using the MarkDuplicates tool in Picard (Broad Institute, 2018) (v2.10.6). Reads mapped to the blacklist regions (ENCODE Project Consortium, 2012) in the mouse mm10 genome (http://mitra.stanford.edu/kundaje/akundaje/release/blacklists/mm10-mouse/mm10.blacklist.bed.gz) were also removed.
Peak calling was performed using epic2 (Stovner and Sætrom, 2019) (v0.0.16), a reimplementation of SICER (Zang et al., 2009). For H3K4me3 and H3K27ac, we used the following parameters: ‘
DM regions of the histone modification ChIP-seq data were called using DiffBind (Ross-Innes et al., 2012) (v2.10.0) in DESeq2 (Love et al., 2014) (v1.22.1) mode. Regions with FDR < 0.05 were considered significant. Genes associated with these regions were identified using GREAT (McLean et al., 2010) (v3.0.0) with the ‘Basal plus extension’ association rule with default parameters. GO enrichment analysis of the associated genes was performed with GREAT, and we considered the GO biological process terms with hypergeometric test FDR < 0.05 as significant.
To select a set of non-DM control regions to match the base levels of H3K27me3 in DM regions, we started with the union peaks of the control and cKO samples and removed any peaks that overlap the DM regions. From these non-DM peaks, we used the R package MatchIt (Ho et al., 2011) (v3.0.2) to select regions with matched peak lengths and H3K27me3 levels (in the unit of RPKM) as those in the DM regions, with the greedy nearest neighbor matching using logistic link propensity score as a distance measure (requiring the standard deviation of the distance to be less than 0.01).
Definition of bivalent and active CGI promoters
CGI promoters were defined as CpG islands (downloaded from the UCSC genome browser) that overlapping with promoters (±2 kb regions around transcription start sites annotated in GENCODE vM10). These CGI promoters were further tested to see whether they overlapped with the ChIP-seq peaks of H3K4me3 and H3K27me3. Bivalent CGI promoters were defined as CGI promoters that overlapped with both the H3K4me3 and H3K27me3 peaks, whereas active CGI promoters were defined as CGI promoters that overlapped with only the H3K4me3 peaks but not the H3K27me3 peaks.
Identification of DMVs
To find DMVs, we first identified UMRs (undermethylated regions) and LMRs (low methylated regions) using MethylSeekR (Burger et al., 2013) (v1.22.0) with m=0.3 (for P39) or 0.5 (for P0), n=7 and FDR < 0.05. In P39 samples PMDs (partially methylated domains) were excluded from further consideration, and no PMDs were found in P0 samples. DMVs were then defined as UMRs with length ≥5 kb and mean methylation level ≤15%. To compare DMVs identified in the P39
For the clustering visualization in Figure 7C, we sorted DMVs according to the following criteria: (1) whether they overlap an H3K27me3 DM region; whether they overlap an up- or down-regulated DE gene (with absolute log2(fold-change) > 0.2); mean mCG > 0.3 in P39 control, P39 cKO; P0 control, or P0 cKO; and finally by the average level of H3K27me3, H3K3me3, and H3K27ac. The ‘Fraction DM H3K27me3’ shows what fraction of the length of the DMV overlaps a DM H3K27me3 region (called using DiffBind). And we also plotted the mean mRNA logFC for all genes contained within each DMV.
DMR enrichment around CGI promoter
We used regioneR (Gel et al., 2016) (v1.14.0) to test whether two sets of genomic regions had significantly higher numbers of overlaps compared to expected by chance. We used permTest() to perform the permutation test, and used the randomizeRegions() function to generate the shuffled control for 5000 times, where the query regions were randomly placed along the genome independently while maintaining their size. The strength of the association of the two sets of regions was estimated using z-score, the distance (measured in standard deviation) between the expected overlaps in the shuffled control and the observed overlaps, and the p-value was reported. To check if the association was specifically linked to the exact position of the query regions, we used the localZscore() function with window = 5000 and step = 50, which shifted the query regions and estimated how the value of the z-score changed when moving the regions.
Other tools used in the data analysis
Browser representations were created using AnnoJ (Lister et al., 2009). All analyses were conducted in R (v3.5.0), MATLAB 2017a and Python 3. Genomic ranges manipulation was done either with bedtools (Quinlan and Hall, 2010) or GenomicRanges (Lawrence et al., 2013). To generate randomly shuffled regions used in Figure 4—figure supplement 3A and Figure 5—figure supplement 1A, we used the sub-command ‘shuffle’ in the bedtools suite, with the ‘-excl’ parameter to exclude the blacklist regions, the ‘-noOverlapping’ parameter to prevent overlapping shuffled regions, and the optional ‘-incl’ parameter to limit the shuffle regions to reside within the P39 H3K27me3 peak regions (the union of control and cKO peaks). Multiple comparison correction for p-values was performed with the Benjamini-Hochberg FDR method (Benjamini and Hochberg, 1995). Results with FDR < 0.05 were considered significant except where stated otherwise. The smoothed lines in Figure 5D and Figure 6—figure supplement 1C were fitted with a generalized additive model using the ‘gam’ function in the ‘mgcv’ R package, with formula = y ~ s(x, bs = ‘cs’).
Data access
All sequencing data are available in the Gene Expression Omnibus under accession GSE141587. A genome browser displaying the sequencing data is available at https://brainome.ucsd.edu/annoj/mm_dnmt3a_ko/.
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer
© 2022, Li, Pinto-Duarte et al. This work is published under https://creativecommons.org/licenses/by/4.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.
Abstract
Two epigenetic pathways of transcriptional repression, DNA methylation and polycomb repressive complex 2 (PRC2), are known to regulate neuronal development and function. However, their respective contributions to brain maturation are unknown. We found that conditional loss of the de novo DNA methyltransferase
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer