Dnmt3a knockout in excitatory neurons impairs postnatal synapse maturation and increases the repressive histone modification H3K27me3

  1. Junhao Li
  2. Antonio Pinto-Duarte
  3. Mark Zander
  4. Michael S Cuoco
  5. Chi-Yu Lai
  6. Julia Osteen
  7. Linjing Fang
  8. Chongyuan Luo
  9. Jacinta D Lucero
  10. Rosa Gomez-Castanon
  11. Joseph R Nery
  12. Isai Silva-Garcia
  13. Yan Pang
  14. Terrence J Sejnowski
  15. Susan B Powell
  16. Joseph R Ecker  Is a corresponding author
  17. Eran A Mukamel  Is a corresponding author
  18. M Margarita Behrens  Is a corresponding author
  1. Department of Cognitive Science, University of California, San Diego, United States
  2. Computational Neurobiology Laboratory, Salk Institute for Biological Studies, United States
  3. Genomic Analysis Laboratory, Salk Institute for Biological Studies, United States
  4. Bioinformatics and Systems Biology Graduate Program, University of California, San Diego, United States
  5. Waitt Advanced Biophotonics Core, Salk Institute for Biological Studies, United States
  6. Howard Hughes Medical Institute, Salk Institute for Biological Studies, United States
  7. Department of Psychiatry, University of California, San Diego, United States

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 Dnmt3a in mouse excitatory neurons altered expression of synapse-related genes, stunted synapse maturation, and impaired working memory and social interest. At the genomic level, loss of Dnmt3a abolished postnatal accumulation of CG and non-CG DNA methylation, leaving adult neurons with an unmethylated, fetal-like epigenomic pattern at ~222,000 genomic regions. The PRC2-associated histone modification, H3K27me3, increased at many of these sites. Our data support a dynamic interaction between two fundamental modes of epigenetic repression during postnatal maturation of excitatory neurons, which together confer robustness on neuronal regulation.

Editor's evaluation

In this manuscript the authors conditionally knock out the DNA methyltransferase Dnmt3a in developing excitatory cortical neurons to determine the consequences for chromatin regulation, gene expression, and neuron function. As expected they find widespread loss of DNA methylation but also an increase in histone methylation (H3K27me3) at many similar regions of the genome, which they propose may be a mechanism of functional compensation. Overall this study offers new insights into the gene regulatory and neuronal cellular functions of an important chromatin regulatory process.

https://doi.org/10.7554/eLife.66909.sa0

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 Dnmt3a mRNA and protein peaks during the second postnatal week (Lister et al., 2013; Stroud et al., 2017), a time of intense synaptogenesis and neuronal maturation. Despite evidence for a unique role of Dnmt3a and mCH in epigenetic regulation of developing neurons, the long-term consequences of Dnmt3a-mediated methylation on brain function remain largely unknown (Stroud et al., 2017).

One challenge in investigating the developmental role of Dnmt3a has been the lack of adequate animal models. Deleting Dnmt3a around embryonic day 7.5 (E7.5) driven by the Nestin promoter dramatically impaired neuromuscular and cognitive development and led to early death (Nguyen et al., 2007). This early loss of Dnmt3a specifically affects the expression of long genes with high levels of gene body mCA (Boxer et al., 2020; Kinde et al., 2016). By contrast, deletion of Dnmt3a starting around postnatal day 14 driven by the Camk2a promoter caused few behavioral or electrophysiological phenotypes (Feng et al., 2010), with only subtle alterations in learning and memory depending on genetic background (Morris et al., 2014). These results suggest that Dnmt3a may play a critical role during a specific time window between late gestation and early postnatal life. During these developmental stages, regulated gains and losses of DNA methylation throughout the genome establish unique epigenomic signatures of neuronal cell types (Qian et al., 2017; Luo et al., 2017; Mo et al., 2015).

To address the role of Dnmt3a-dependent epigenetic regulation in the functional maturation of cortical excitatory neurons, we created a mouse line using the Neurod6 promoter (Schwab et al., 2000) (Nex-Cre) to delete exon 19 of Dnmt3a (Okano et al., 1999). In this conditional knockout (cKO), Dnmt3a is functionally ablated in excitatory neurons in the neocortex and hippocampus starting in mid-to-late gestation (embryonic day E13–15) (Goebbels et al., 2006). In Dnmt3a cKO animals, DNA methylation was substantially disrupted in excitatory neurons, leading to altered behavior and synaptic physiology without early life lethality or overt brain morphological alterations. We generated deep DNA methylome, transcriptome, and histone modification data in Dnmt3a cKO and control pyramidal cells of the frontal cortex (Supplementary file 1), enabling a detailed assessment of the molecular basis of neurophysiological and behavioral phenotypes. We found that the polycomb repressive complex 2 (PRC2)-associated chromatin modification, H3K27me3, increases during postnatal development following the loss of mCH DNA methylation in Dnmt3a cKO neurons primarily at sites where mCG is also depleted. Our data support a dynamic interaction between two fundamental modes of epigenetic repression in developing brain cells.

Results

Dnmt3a cKO in pyramidal neurons during mid-gestation specifically impairs working memory, social interest, and acoustic startle

Previous studies of Dnmt3a KO mice yielded results ranging from little cognitive or health effects (Feng et al., 2010; Morris et al., 2014) to profound impairment and lethality (Dura et al., 2022; Nguyen et al., 2007). These studies suggest that the developmental timing and cell type of Dnmt3a loss may determine the extent of subsequent phenotypes. Here, we took a targeted approach by functionally ablating Dnmt3a in cortical pyramidal cells starting during mid-gestation (Figure 1A). We took advantage of the developmental onset of Neurod6 expression between embryonic day E11 and E13, after the onset of Nestin expression (Thompson et al., 2014) but well before the major postnatal wave of reprogramming of the neuronal DNA methylome (Figure 1—figure supplement 1; Lister et al., 2013). We confirmed that Neurod6-dependent Cre recombination occurred only in excitatory neurons using the INTACT (isolation of nuclei tagged in specific cell types) (Mo et al., 2015) mouse (Figure 1—figure supplement 2). The deletion of Dnmt3a exon 19 was faithfully captured in cortical excitatory neurons of cKO animals (Figure 1B and Figure 1—figure supplement 3A), which was known to produce a deletion of 50 amino acids in the methyltransferase domain of the DNMT3A protein (Lyko, 2018). We also verified the reduction in DNMT3A protein in whole tissue extracts at early postnatal time points (P5 and P13) when Dnmt3a mediated accumulation of mCH normally begins in the frontal cortex (Lister et al., 2013; Figure 1—figure supplement 3B). Dnmt3a cKO animals survived and bred normally, without overt morphological alterations in the brain (Figure 1—figure supplement 3C). We found no impairments in gross motor function in an open field test: cKO mice traveled more during the first 5 min of testing (Figure 1—figure supplement 4A) while performing fewer rearings associated with exploratory interest (Figure 1—figure supplement 4B). Moreover, cKO mice had no signs of increased anxiety-like behavior on three separate behavioral tests (Figure 1—figure supplement 4C-E), in contrast with the reported anxiogenic effects of Dnmt3a knockdown in the mPFC of adult mice (Elliott et al., 2016). The absence of major impairments in overall health, motor function, or anxiety-like behavior established a baseline for investigating the role of Dnmt3a in specific cognitive and social behaviors.

Figure 1 with 5 supplements see all
Dnmt3a conditional knockout (cKO) in cortical pyramidal neurons during mid-gestation impaired working memory, social interest, and acoustic startle responses.

(A) An experimental model of the conditional loss of Dnmt3a in excitatory neurons. P0 and P39, postnatal day 0 and 39. FANS, fluorescence-activated nuclei sorting. (B) RNA-seq confirmation of the deletion of Dnmt3a exon 19 in P39 excitatory neurons. RPKM, reads per kilobase per million. R1/2, replicate 1/2. *, t-test p=0.014. (C) Dnmt3a cKO mice made fewer spontaneous alternations in the Y-maze test of working memory (Wilcoxon test, **, p=0.0079; *, p=0.011; n=15 male control, 15 male cKO, 11 female control, 10 female cKO). (D) Male Dnmt3a cKO mice spent less time interacting with an unfamiliar mouse, indicating reduced social interest (Wilcoxon test; *, p=0.01048; **p=0.006833; n=14 male control, 15 male cKO, 11 female control, 10 female cKO). (E) Male Dnmt3a cKO mice had decreased startle response to a 120 dB acoustic pulse (Wilcoxon test, **, p=0.0019; n.s., not significant; n=14 male control, 15 male cKO, 11 female control, 10 female cKO).

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). Dnmt3a cKO mice did not alternate spontaneously between the arms of a Y-maze (p=0.0079 for males, p=0.011 for females, Figure 1C), indicating impaired spatial working memory. Moreover, when tested in a three-chamber box in which one of the sides contained a novel mouse and the opposite a novel object, male Dnmt3a cKO animals spent less time in the former, opting, instead, to remain significantly longer in the center (empty compartment), which is suggestive of reduced exploration and social interest (Figure 1D, left panel; p=0.01048). Male Dnmt3a cKO mice also had significantly attenuated acoustic startle reflex (p=0.0019, Figure 1E and Figure 1—figure supplement 5A). We observed increased prepulse inhibition (PPI) in male Dnmt3a cKO mice, but this may be driven by the reduced startle reflex (Figure 1—figure supplement 5B, p<0.05). It is noteworthy that the observed deficits in startle response were not due to impaired hearing, since Dnmt3a cKO mice displayed intact hearing in tests of PPI and fear conditioning.

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 Dnmt3a cKO and control male mice in acquisition or recall of fear memory following re-exposure to the context or conditioned stimulus after 24–48 hr (Figure 1—figure supplement 5C-E), or in extinction (Figure 1—figure supplement 5F). Altogether, these behavioral results indicate that Dnmt3a cKO in excitatory neurons impairs working memory, social interest, and acoustic startle, without generalized cognitive disruption.

Loss of Dnmt3a impairs synapse maturation and attenuates neuronal excitability

To test the impact of Dnmt3a cKO on dendritic morphology, we quantified the number and structure of 1278 DiI-labeled dendritic spines (NexCre/C57: n=701 from 5 mice; cKO: n=577 from 4 mice) of layer 2 pyramidal neurons of the mouse prelimbic region (~2 mm anterior to Bregma) (Figure 2A; Materials and methods), a region critical for working memory (Yang et al., 2014) and social approach behavior (Lee et al., 2016). While the overall density of dendritic spines was equivalent in control and Dnmt3a cKO neurons (Figure 2B), the spines were significantly longer (mean length 2.219±0.052 µm in cKO, 1.852±0.034 µm in control) and narrower (mean width 0.453±0.008 µm in cKO, 0.519±0.008 µm in control) in Dnmt3a cKO neurons (Figure 2C, KS test p<0.001; Figure 2—figure supplement 1). Consistent with this, a larger proportion of spines in Dnmt3a cKO mice were classified as immature filopodia, and fewer were mushroom-shaped mature spines (Figure 2D) according to pre-established morphometric criteria (see Materials and methods). The proportion of spines with other morphologies, including branched spines with more than one neck (data not shown), was not significantly different between genotypes (Figure 2D). These data indicate a role for Dnmt3a in dendritic spine maturation.

Figure 2 with 2 supplements see all
Immature spine morphology and reduced excitability of layer 2 excitatory neurons following Dnmt3a conditional knockout (cKO).

(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 Dnmt3a cKO mice. (B) The density of membrane protrusions was unchanged in the Dnmt3a cKO (Wilcoxon test, n.s., not significant). (C) Membrane protrusions were significantly longer and narrower in the Dnmt3a cKO (KS test, p<0.001). (D) More spines were classified as immature filopodia, and fewer as mature mushroom-shaped spines with large postsynaptic densities in the Dnmt3a cKO (Wilcoxon test, **, p=0.0015; *, p=0.046 and 0.011 for thin and mushroom, respectively). (E) Example whole-cell patch-clamp recordings from prelimbic layer 2 pyramidal neurons following 60 pA current injections. (F) The median rheobase (i.e. the minimal current necessary to elicit an action potential) was significantly higher in the Dnmt3a cKO (t-test, **, p=0.0042). (G) Action potential frequency vs. injected current (mean ± SEM) showed reduced excitability in Dnmt3a cKO (Wilcoxon test, *, p<0.05). (H) and (I) Dnmt3a cKO neurons were slightly hyperpolarized at Vrest when compared to control (Wilcoxon test, *, p=0.049) and had lower membrane resistance (Wilcoxon test, *, p=0.023).

To test how the loss of Dnmt3a and subsequent stunting of spine maturation affect intrinsic neuronal excitability and synapse sensitivity, we next performed patch-clamp experiments in visually identified layer 2 pyramidal neurons from the prelimbic region (Figure 2E). Whole-cell current-clamp recordings showed that Dnmt3a cKO neurons (n=22 cells from 11 mice) required greater current injections than control (n=17 cells from 12 mice) to trigger an action potential (higher rheobase, t-test p=0.0042, Figure 2F), though there was no difference in membrane potential at the firing threshold (Figure 2—figure supplement 2A). Dnmt3a cKO neurons also produced fewer spikes in response to injected current (Figure 2G). These neurons were slightly hyperpolarized at rest (mean –70.90±0.8 mV in cKO, n=22 cells from 11 mice vs. –67.22±1.4 mV in control, n=17 cells from 12 mice, Wilcoxon test p=0.049, Figure 2H), which could reflect differential expression of ion channels at the plasma membrane. Consistent with this, Dnmt3a cKO neurons had a lower input resistance (Wilcoxon test p=0.023, Figure 2I), suggesting increased expression of functional transmembrane ion channels. Whole-cell voltage-clamp recordings of miniature excitatory postsynaptic currents (mEPSCs) showed slight, yet significant, increased amplitude variability in Dnmt3a cKO mice (Figure 2—figure supplement 2, 8.77±0.32 pA in cKO, n=8 cells from 4 mice vs. 8.67±0.089 pA in control, n=8 cells from 5 mice, F-test, p=0.0032), consistent with disruption at postsynaptic sites. However, we found no alteration in the mean amplitude (Figure 2—figure supplement 2B) or frequency (Figure 2—figure supplement 2C) of mEPSCs recorded at the soma.

Altered gene expression in Dnmt3a cKO excitatory neurons

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 Dnmt3a cKO animals into the INTACT mouse (Mo et al., 2015) on a C57BL/6J background, followed by fluorescence-activated nuclei sorting (FANS) and RNA-seq. Although sorted nuclei contain only a subset of the cell’s total mRNA and are enriched in immature transcripts, nuclear RNA-seq is nevertheless a quantitatively accurate assay of gene expression that is robust with respect to neural activity-induced transcription (Bakken et al., 2018; Lacar et al., 2016). All of the molecular assays were applied to at least two independent biological samples, with each sample being derived from a pool of tissue from 1 or 2 mice (Supplementary file 1). Nuclear RNA abundance was highly consistent across independent replicates within the same group (Spearman correlation r=0.93–0.94, Figure 3—figure supplement 1A). Given the repressive role of DNA methylation (mCG and mCH) in regulating gene expression in neurons (Guo et al., 2014; Lister et al., 2013), we expected to find increased gene expression in the cKO. Consistent with this, we detected 46 differentially expressed (DE) genes with higher expression in the cKO (false discovery rate [FDR] < 0.05, Figure 3A–B and Figure 3—figure supplement 1B, Supplementary file 2). We also detected significantly lower expression of 24 genes in cKO neurons (Figure 3A). Several of the DE genes had annotated roles in dendrite morphogenesis (Elavl4, Hecw2, Ptprd), as well as in the regulation of Na+ (Hecw2, Scn3b) and Ca2+ levels (Cacnb3) (Figure 3B). Early life experiences were shown to alter the expression of Dnmt3a and hence the DNA methylation levels in some transposable elements (TEs) in the hippocampus (Bedrosian et al., 2018). We examined the expression of transposon families or classes in our data but found no significant changes between cKO and control in the cortical excitatory neurons (two-sample t-test, FDR > 0.05; Figure 3—figure supplement 2).

Figure 3 with 5 supplements see all
Loss of Dnmt3a leaves thousands of genomic regions in a fetal-like demethylated state.

(A) 70 genes were differentially expressed (DE) (false discovery rate [FDR] < 0.05) in P39 pyramidal neurons in Dnmt3a conditional knockout (cKO) vs. control. Top, fold-change (FC) between cKO and control; bottom, heatmap showing normalized expression of the DE genes in each sample. Z-scores were computed using mRNA counts per million (CPM) for each DE gene. (B) Differential gene expression in control vs. Dnmt3a cKO excitatory neurons at P39. Significant up-regulated and down-regulated DE genes are shown in red and blue, respectively. Differentially expressed (DE) genes associated with dendrite morphogenesis (Elavl4, Hecw2, Ptprd), and regulation of Na+ (Hecw2, Scn3b) and Ca2+ levels (Cacnb3) are labeled. (C) Non-CG DNA methylation (mCH) is eliminated, and mCG is reduced, in P39 Dnmt3a cKO pyramidal cells, while mCG and mCH levels are not changed in P0 (t-test: *, p<0.05; n.s., not significant). P0 and P39, postnatal days 0 and 39, respectively. Each bar represents the methylation level in one replicate. (D) Non-CG DNA methylation (mCH) in P39 pyramidal cells in control samples, in 1 kb bins in the flanking 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. (E) The difference in gene body methylation vs. fold-change of gene expression between P39 Dnmt3a cKO and control. The plots show mean ± SEM gene expression fold-change for genes in 10 non-overlapping bins (deciles of mC difference). (F) The Nedd4 promoter locus contains five differentially methylated regions (DMRs, yellow horizontal rectangles and shaded in blue boxes) with naive, fetal-like mCG in P39 Dnmt3a cKO. Ticks show mCG at CG sites. Four out of the five P39 Dnmt3a cKO DMRs overlapping developmental gain-of-methylation DMRs (red horizontal rectangles) are marked with arrows. CGI, CpG island. R1 and R2, replicates 1 and 2. (G) Overlap of P39 Dnmt3a cKO DMRs and developmental DMRs. (H) P39 Dnmt3a cKO hypo-DMRs are significantly enriched (depleted) in DMRs that normally gain (lose) methylation during development (Fisher’s test, p<0.05).

Dnmt3a cKO abolishes postnatal DNA methylation

Deletion of Dnmt3a during mid-gestation should disrupt the subsequent gain of DNA methylation at specific genomic sites during development (Qian et al., 2017; Lister et al., 2013), without affecting sites that maintain or lose methylation after E14.5. Using single base resolution, whole-genome MethylC-seq (Lister et al., 2008) in biological replicates with strong consistency (Figure 3—figure supplement 3A), we confirmed that non-CG DNA methylation (mCH) in excitatory neurons is absent at birth (<0.1% of all CH sites at P0), and accumulates by postnatal day 39 (1.98% at P39) (Lister et al., 2013). The cKO all but eliminated mCH (<0.1% at P0 and P39) (Figure 3C and Figure 3—figure supplement 3B). While mCH increases in neurons during postnatal life, the genome-wide level of mCG in the brain remains high throughout the lifespan (Lister et al., 2013). We found that the genome-wide mCG level was 12.5% lower in mature (P39) cKO neurons (60.1% in cKO vs. 72.6% in control). There was no difference in mCG in newborn mice (P0, 73.1% in both cKO and control, Figure 3C). mCG at P39 was reduced in 92.2% of all genomic bins (10 kb resolution), and was significantly lower in introns, 3’ UTR, and intergenic regions (Figure 3—figure supplement 3C). The reduction in mCG was strongly correlated with reduced mCH (Spearman correlation r=0.805, p<10–3, Figure 3—figure supplement 3D). These data support a role for Dnmt3a in postnatal de novo CG and CH DNA methylation across genomic compartments (Lister et al., 2013; Stroud et al., 2017).

Reduced DNA methylation does not fully explain altered transcription in Dnmt3a cKO

We investigated whether the altered gene expression in Dnmt3a cKO neurons correlated with loss of DNA methylation at specific sites. We first analyzed DNA methylation around DE genes in mature neurons (P39). The simple model of DNA methylation as a repressive regulator of gene expression predicts that genes that lose the most mC should be most transcriptionally up-regulated in the cKO. Consistent with this, we found that mCH was strongly enriched in the gene body of up-regulated genes in control neurons (Figure 3D and Figure 3—figure supplement 4A). By contrast, genes with similar expression levels in the control neurons which were not transcriptionally up-regulated in the cKO (Figure 3—figure supplement 4B) had significantly lower gene body mCH (Figure 3D). Moreover, down-regulated genes had low gene body mCH. These data support a causal role for gene body mCH in repressing gene expression. The relatively lower mCH level in down-regulated genes could make them less sensitive to the loss of Dnmt3a. The dysregulation of their expression may be due to secondary effects subsequent to the direct loss of DNA methylation.

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 Dnmt3a. Up-regulated genes were on average longer than down-regulated genes (Figure 3—figure supplement 4E, Wilcoxon rank-sum test p<10–4) and non-DE genes (p<0.01), consistent with the reported enrichment of mCA and MeCP2-dependent gene repression in long genes (Boxer et al., 2020; Gabel et al., 2015; Kinde et al., 2016). However, there was a broad distribution of gene lengths for both up- and down-regulated genes (Lavery et al., 2020).

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 Nedd4 (Figure 3F). Four of these DMRs were also unmethylated in excitatory neurons in newborn mice, showing that the loss of Dnmt3a blocked the developmental gain of mCG at these sites. However, the density of P39 cKO DMRs around the DE genes was not significantly different in the up- and down-regulated genes and the non-DE genes (Figure 3—figure supplement 5B).

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 Rest, Lhx2, Pou3f2(Brn2), and Pax6 (FDR < 0.05, Figure 3—figure supplement 5D, Supplementary file 4). Notably, the DMRs in Dnmt3a cKO neurons represent only a part of the global reduction in mCG that we observed throughout the genome. Indeed, we found that mCG is reduced by ~10% in all genomic compartments, even after excluding P39 DMRs (Figure 3—figure supplement 3C). These results suggest that Dnmt3a is essential for the methylation and subsequent repression of neuronal enhancers that are active during prenatal brain development.

Increased PRC2-associated repressive histone modification H3K27me3 in Dnmt3a cKO

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 Ezh2, Suz12, Eed, and Jarid2, were among the top DNA binding proteins enriched near the promoters of both up- and down-regulated DE genes (Figure 4A). Several TFs associated with chromatin organization, including the histone deacetylase (Hdac) and demethylase (Kdm) families, and Ctcf, were also enriched (Figure 4—figure supplement 1, Supplementary file 5). These results suggest that Dnmt3a cKO impacts the chromatin landscape in excitatory neurons, potentially via altered PRC2 activity.

Figure 4 with 3 supplements see all
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 Dnmt3a conditional knockout (cKO) differentially expressed genes include many proteins associated with PRC2. The functional TF rank score was assigned by Binding Analysis of Regulation of Transcription (BART; Wang et al., 2018). PRC2-associated TFs are labeled and highlighted in gold circles. (B) Browser view of the Mab21l2 locus, where increased H3K27me3 (differentially modified regions, bottom red bars and highlighted in blue shaded box) coincides with the loss of DNA methylation (Dnmt3a cKO DMRs, orange bars under the ‘P39 cKO’ track and highlighted in blue shaded box) in P39 Dnmt3a cKO. This region loses H3K27me3 during normal development in control pyramidal neurons (blue bars, P39 < E14). DMR, differentially methylated region; E14, embryonic day 14; P0 and P39, postnatal days 0 and 39. (C) Quantification of the increase in H3K27me3 chromatin immunoprecipitation sequencing (ChIP-seq) signal in each replicate at the H3K27me3 differentially modified region between P39 control and cKO at the Mab21l2 locus shown in (B). Each bar shows the DEseq2 normalized counts in each replicate, and the triple asterisks denote a significant increase (false discovery rate [FDR] = 1.33e-4, fold-change=2.41). (D) Histone modification ChIP-seq peaks for active marks (H3K4me3, H3K27ac) are largely preserved in the Dnmt3a cKO, while repressive H3K27me3 peaks expand. The Venn diagrams denote numbers of peaks that overlap between cKO (yellow) and control (black) (numbers in the center), and numbers of peaks that are unique to one of the conditions (numbers on the edges).

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 Mab21l2 locus (Figure 4B), we observed a 2.41-fold increase (Figure 4C) in H3K27me3 in P39 Dnmt3a cKO neurons, coinciding with the loss of CG methylation at multiple DMRs spanning the gene body and surrounding region. The Mab21l2 gene was lowly expressed (TPM < 3) in both control and cKO neurons, consistent with a role for the gain of H3K27me3 in maintaining the repression of this gene.

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 Dnmt3a cKO. Genes associated with these DM regions were enriched in development-related functions (Figure 4—figure supplement 3B, Supplementary file 8). These results suggest that the increase of H3K27me3 in Dnmt3a cKO excitatory neurons occurred postnatally, following the major impact of the loss of Dnmt3a on neuronal DNA methylation.

The increase in H3K27me3 was closely associated with regions that lost CG DNA methylation in the Dnmt3a cKO. The majority (66.2%) of the regions marked by H3K27me3 in the cKO overlapped with P39 cKO DMRs (Figure 5A, Fisher’s test p<1e-100). Likewise, P39 cKO DMRs were significantly enriched in peaks (22.4% overlapped H3K27me3 peaks in cKO, p<1e-100, Figure 5B) and DM regions of H3K27me3 (Figure 4—figure supplement 3D). The DM regions of H3K27me3 had significantly more overlaps with DMRs compared to the non-DM control regions (Fisher’s exact test p<2.2e-16, OR: 2.70, Figure 4—figure supplement 3C-D). The DM regions also have larger decreases in both mCG and mCH when compared to non-DM regions (Wilcoxon rank-sum test p<0.0001, Figure 4—figure supplement 3E). Conversely, the DMRs were depleted in regions marked by H3K27ac (12.7%, p<1e-100) or H3K4me3 (0.98%, p=1.73e-149) (Figure 5A–B). Moreover, H3K27me3 was more abundant at the center of DMRs in cKO compared to control neurons at P39 (∆=0.15 in the unit of fold enrichment vs. IgG, 5.04% increases, Figure 5C). There were much smaller changes of H3K27me3 at these DMRs in newborn (P0, ∆=0.0075, 0.27% changes) or fetal (E14, ∆=–0.028, –1.02% changes) neurons (Figure 5C), and the increases were not seen in randomly shuffled regions (Figure 5—figure supplement 1A). The signal of H3K4me3 and H3K27ac at the DMRs was also elevated in the cKO, to a lesser extent (H3K4me3: ∆=0.068, 2.92% changes; H3K27ac: ∆=0.13, 4.64% changes). Going beyond overlaps of regions, we found a quantitative association between the changes in DNA methylation and H3K27me3 in mature (P39) neurons (Figure 5D). At H3K27me3 peaks, the ChIP-seq signal intensity fold-change between cKO and control correlated with the loss of mCG in cKO (Spearman r=–0.26, p<2.2e-16). When accessing the H3K27me3 changes as a function of baseline H3K27me3 levels in the control across the genome tiled in 1 kb bins, we observed bigger H3K27me3 increases in bins with overlapping DMRs compared to bins without overlapping DMRs (Wilcoxon rank-sum test p<0.001 in each bin, Figure 5—figure supplement 1B). These results indicate that DMRs are particularly unique in showing increased H3K37me3 signal.

Figure 5 with 2 supplements see all
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 Dnmt3a cKO differentially methylated regions (DMRs). (B) Significant enrichment (red) or depletion (blue) of P39 Dnmt3a cKO DMRs in the histone modification chromatin immunoprecipitation sequencing (ChIP-seq) peaks (Fisher’s test, p<0.05). E14, embryonic day 14; P0 and P39, postnatal days 0 and 39. (C) Histone modification ChIP-seq signal around the center of DMRs. RPKM, reads per kilobase per million. ***, Wilcoxon rank-sum test of the differences at the center, p<0.001. (D) Correlation of P39 H3K27me3 signal fold-changes and P39 CG methylation levels differences between Dnmt3a cKO and control in H3K27me ChIP-seq peaks. The smoothed line is fitted using a generalized additive model, and the shaded area shows the 95% confidence interval of the fit. r, Spearman correlation coefficient. (E) P39 Dnmt3a cKO down-regulated DE genes (false discovery rate [FDR] < 0.05) with overlapping P39 Dnmt3a cKO DMRs show small but significant increases of H3K27me3 in P39 cKO (upper panel). *, Wilcoxon rank-sum test against zero, p<0.05. No such differences were observed in differentially expressed (DE) genes that do not overlap with P39 cKO DMRs (bottom panel). H3K27me3 signal was calculated as the RPKM fold-change between H3K27me3 and IgG.

We further examined whether the increased H3K27me3 could account for the reduced expression of some genes in the Dnmt3a cKO neurons. We found that down-regulated genes that overlap DMRs and the non-DE genes selected with matched expressions of up-regulated DE genes that overlap DMRs showed a small but significant median increase of H3K27me3 in the cKO (Wilcoxon rank-sum test p-value < 0.05; Figure 5E upper panel). When we considered a larger set of DE genes with a more relaxed threshold (FDR < 0.2), we observed that down-regulated DE genes containing DMRs accumulate significantly more H3K27me3 than non-DE genes containing DMRs and up-regulated genes containing DMRs (Wilcoxon rank-sum test p=0.0029, Figure 5—figure supplement 2A). By contrast, up-regulated genes had no significant accumulation of H3K27me3 and were not significantly different from the control genes (Figure 5—figure supplement 2A). No such differences were observed between DE and non-DE genes without overlapping DMRs (Figure 5E lower panel and Figure 5—figure supplement 2A). These results suggest that the effect of H3K27me3 is likely specific to genes containing DMRs and the effect is stronger in the down-regulated DE genes, which may partially explain the fact that 24 genes were significantly down-regulated after the loss of repressive DNA methylation in the Dnmt3a cKO (Figure 2A–B).

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 Dnmt3a cKO

Our finding that Dnmt3a cKO disrupts the normal developmental gain of DNA methylation prompted us to ask whether the changes in H3K27me3 in the cKO are likewise associated with developmental regulation of H3K27me3. Indeed, our ChIP-seq data from E14, P0, and P39 excitatory neurons revealed striking developmental dynamics in H3K27me3. We identified 12,994 developmentally regulated H3K27me3 DM regions between E14 and P39, with a similar number of regions that gain (6774) and lose (6220) H3K27me3 (Figure 6A and Figure 6—figure supplement 1A, Supplementary file 9). We also examined DM regions at P0 vs. E14 and P0 vs. P39 (Figure 6—figure supplement 1A). However, due to greater biological variability at the perinatal time point, the two replicate ChIP-seq datasets at P0 were less consistent than the E14 and P39 samples (Figure 6—figure supplement 2). As a result, our data were not well powered to detect changes at P0 and we chose to focus on the E14 vs. P39 DM regions as sites of developmental chromatin remodeling. Genes associated with these DM regions were enriched in biological processes involved in development such as nervous system development and neurogenesis (Figure 6—figure supplement 1B). The developmental H3K27me3 DM regions overlapped developmental DMRs, with a notable overlap of regions gaining both mCG and H3K27me3 (Figure 6B). Both modifications may thus act together to repress thousands of genomic regions during development.

Figure 6 with 2 supplements see all
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 Dnmt3a cKO (E).

The P39 Dnmt3a cKO changes of H3K27me3 signal were weakly correlated with the developmental changes of H3K27me3 signal (r = –0.12, p<2.2e-16, Figure 6—figure supplement 1C). Moreover, there was little appreciable difference between the H3K27me3 signal fold-changes between P39 and E14 cKO samples compared with those of control samples (Spearman r=0.67, p<2.2e-16, Figure 6—figure supplement 1D).

To further stratify the joint distribution of developmental and Dnmt3a cKO-dependent changes in both H3K27me3 and DNA methylation, we assigned peaks to three groups (Figure 6C–E and Figure 6—figure supplement 1E-G). ‘Group DevLoss’ and ‘Group DevGain’ peaks lose or gain H3K27me3 during development, respectively (Figure 6C–D and Figure 6—figure supplement 1E-F). ‘Group cKO’ peaks have higher H3K27me3 in the Dnmt3a cKO compared to control at P39 (Figure 6E and Figure 6—figure supplement 1G). We found that developmental peaks (DevLoss and DevGain) were relatively unaffected by the cKO (ΔH3K27me3=0.02, –0.03 respectively, in units of log10(CPM +1)), whereas Group cKO had 4.5-fold larger mean effect (ΔH3K27me3=0.09). Group cKO peaks also experienced greater loss of mCG (ΔmCG = –23.5%) than Group DevLoss (–13.0%) or DevGain (18.6%) (Figure 6C–E and Figure 6—figure supplement 1E-G, middle and right panels). These results suggest that regions prone to alteration of H3K27me3 by Dnmt3a cKO are distinct from the regions affected by developmentally dynamic H3K27me3. We observed a very similar pattern when we examined sites with differential H3K27me3 between P0 and P39 (data not shown).

Novel DNA methylation valleys with increased H3K27me3 signal in the Dnmt3a cKO

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 with 1 supplement see all
Distinct clusters of DNA methylation valleys (DMVs) were associated with the increased H3K27me3 signal in the Dnmt3a conditional knockout (cKO).

(A) Number of DMVs identified in the P39 Dnmt3a cKO and the control samples, categorized by whether they appear in one or both groups or change size in the cKO. (B) Overlap of DMVs with the H3K4me4 and/or the H3K27me3 chromatin immunoprecipitation sequencing (ChIP-seq) peaks. (C) Heatmap of DMVs clustered by their methylation levels and histone modifications. The last two columns show the enrichments of differentially modified (DM) peaks of H3K27me3 and differentially expressed (DE) genes. R1/2, replicate 1/2; RPKM, reads per kilobase per million. (D) Browser tracks show examples of unique DMVs in the Dnmt3a cKO samples and the increased H3K27me3 signal in their flanking regions.

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 Tfap2c (Figure 7C–D), had high H3K27me3 and low mCG in both control and cKO neurons, and were not strongly affected by the loss of mC in the Dnmt3a cKO. By contrast, cluster C1 and C4 DMVs, including the promoters of Lhx2, Foxp1, Foxp2, Slc17a6 (encoding vesicular glutamate transporter, Vglut2) and Sema3f, gain mC during normal development (Figure 7C–D). Cluster C1 DMVs gained H3K27me3 in the P39 cKO compared to control animals. The loss of mC in these regions in the Dnmt3a cKO did not lead to strong activation of gene expression, potentially due to compensatory PRC2-mediated repression.

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. Nr1d1, Pde4d, Figure 7D) have high mCG at P0 and lose methylation in excitatory neurons during brain development. This demethylation is not affected in the Dnmt3a cKO, and we found little difference between the control and cKO neurons at these sites. Finally, clusters C2 and C3 were enriched for up- and down-regulated genes, respectively.

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 Dnmt3a expression and non-CG methylation (Lister et al., 2013). Neurod6 conditional deletion of Dnmt3a in cortical excitatory neurons abolished non-CG methylation and reduced CG methylation throughout the genome (Figure 3C). We found that Dnmt3a cKO neurons had altered expression of dozens of genes, including both up- and down-regulated genes (Figure 3A–B). Similar complex patterns in gene expression were reported when Dnmt3a was deleted in inhibitory neurons (Lavery et al., 2020) and following manipulation of the DNA methylation reader MeCP2 (Boxer et al., 2020; Johnson et al., 2017; Lavery et al., 2020). These gene expression changes may partly reflect the direct effect of lower gene body CH methylation and loss of CG and CH methylation at gene promoters and distal enhancers (Boxer et al., 2020; Clemens et al., 2020). In addition, some gene expression changes could result from the disruption of other regulatory processes, such as TF expression or chromatin modification. Indeed, we found that the Dnmt3a cKO DMRs overlapped with regions that gain methylation during normal postnatal development (Figure 3G–H). These DMRs had increased H3K27ac in the Dnmt3a cKO (Figure 5C), consistent with observations in adult animals lacking Dnmt3a specifically in GABAergic Sst- or Vip-expressing interneurons (Stroud et al., 2020). In those experiments, embryonic gene regulatory elements had lower cytosine methylation and increased H3K27ac and H3K4me1 (Stroud et al., 2020). This suggests an essential role for Dnmt3a and DNA methylation in shaping the transcriptome during development in part via the inactivation of embryonic enhancers, with potentially long-lasting effects on the gene expression pattern of mature neurons. Although we did not detect individual peaks with a statistically significant difference in H3K27ac between cKO and control (Figure 4—figure supplement 1C), this could be due to a small number of replicates and the subtle magnitude of changes in gene expression relative to the large number of tested peaks. Future experiments using more replicates or orthogonal techniques like CUT&RUN and CUT&Tag could help to further investigate the relationship between the Dnmt3a-dependent DNA methylation and the activation of enhancers.

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 Hox gene clusters (Reddington et al., 2013). We did not observe ectopic expression of these genes, possibly due to the relatively modest reduction in mCG in our model compared with cells lacking Dnmt1. Instead, we found that in cortical excitatory neurons, thousands of sites gained H3K27me3 following the loss of mCG in the Dnmt3a cKO (Figure 4—figure supplement 1D). These sites, which normally gain DNA methylation during postnatal development, were left unmethylated in cKO neurons (Figures 4B and 6E). These regions were largely distinct from the sites that gain or lose H3K27me3 during normal development (Figure 6C–E), and had an intermediate level of H3K27me3 in the control samples (Figure 4—figure supplement 3A). A subset of these regions formed large-scale DMVs spanning key regulatory genes (Figure 7C–D). Overall, our results suggest that when DNA methylation is disrupted, H3K27me3 might partially compensate for the loss of mCG and/or mCH and act as an alternative mode of epigenetic repression. Nevertheless, we did not find differential expression in any of the four core components of PRC2 (Ezh2, Suz12, Eed, and Rbbp4) in adult Dnmt3a cKO animals. It is possible that the increased H3K27me3 was mediated by transient expression of PRC2 components during development in the cKO. Furthermore, the predictions from BART (Figure 4A) were derived from various cell lines and tissues from the ENCODE project (Davis et al., 2018; ENCODE Project Consortium, 2012), suggesting that the potential PRC2 binding at our DEGs may normally happen in systems other than the brain or pyramidal neurons, or at other time points during development. Additional experiments which directly manipulate components of the PRC2 system are required to further test the potential compensation mechanism.

Previous work showed that the early embryonic deletion of Dnmt3a in excitatory neurons caused gross motor deficits and a shortened lifespan (Nguyen et al., 2007). Mid-gestation Neurod6-driven Dnmt3a ablation, however, did not cause such alterations, allowing us to test how the epigenetic and transcriptional alterations affected the morphology and function of excitatory synapses, neurons’ passive and active properties, and behavior. We found that prelimbic layer 2 neurons of Dnmt3a cKO animals had more immature dendritic spines and were less sensitive to somatic injections of depolarizing current compared to control neurons (Figure 2). This is consistent with several of the observed DE genes having annotated roles in dendrite morphogenesis (Elavl4, Hecw2, Ptprd, Figure 3B) and Na+ influx/transport (Hecw2, Scn3b, Figure 3B). Indeed, the down-regulation of the latter set of genes could be expected to increase the action potential threshold (thus dampening neuronal excitability) and trigger neurodevelopmental delays (Berko et al., 2017). As well, Cacnb3, a gene that encodes a regulatory beta subunit of the voltage-dependent calcium channel, was found among the set of down-regulated genes after loss of Dnmt3a. This gene has been linked with schizophrenia (Maycox et al., 2009), ADHD, and bipolar disorder (van Hulzen et al., 2017) in humans. Our study thus suggests that the disruption of methylation patterns established by Dnmt3a during infancy might have far-reaching mechanistic relevance for multiple neurodevelopmental disorders, complementing previous genetic risk association studies linking Dnmt3a with autism (C Yuen et al., 2017; Sanders et al., 2015). Although our experiments were not designed to make quantitative comparisons between cohorts of different sexes, we observed a wider range of behavioral impairments in male Dnmt3a cKO mice, which included, for example, decreased social interest (Figure 1). Those data correlate well with several human neurodevelopmental disorders, in which males are disproportionately more affected than females, and thus provide important support for follow-up investigations into the underlying causes of those differences.

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 Dnmt3a cKO mice line

Request a detailed protocol

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, Dnmt3a-floxed animals (Okano et al., 1999) (backcrossed to C57BL/6 for at least seven generations) were crossed to Neurod6-Cre (Nex-Cre) (Goebbels et al., 2006, backcrossed to C57BL/6 J for >10 generations) mice to generate Dnmt3a-KO animals carrying the deletion only in pyramidal cells. To be able to isolate pyramidal neuron nuclei for DNA methylation, transcription, and ChIP analyses, the mouse lines (Dnmt3a-KO and Neurod6-Cre) were crossed to a mouse line carrying the INTACT background (B6.129-Gt(ROSA)26Sortm5(CAG-Sun1/sfGFP)Nat/MmbeJ.Strain 030952, Jackson laboratories). The deletion of Dnmt3a from pyramidal cells was confirmed by RNA-seq (deletion of exon 19) and Western blot (Figure 1B and Figure 1—figure supplement 3A-B). For both backgrounds, Nex-Cre hemizygous mice were used as controls.

Frontal cortex dissection, nuclei isolation, and flow cytometry

Request a detailed protocol

Frontal cortex tissue was produced as described (Lister et al., 2013; Luo et al., 2017) from postnatal day 0 and 39 (P39) Dnmt3a cKO and control animals, in an INTACT background. The nuclei of GFP-expressing NeuN-positive excitatory neurons were isolated and collected using FANS as described (Lister et al., 2013; Luo et al., 2017) with the following modification: prior to FANS, nuclei were labeled with anti-NeuN-AlexaFluor647 and anti-GFP-AlexaFluor488. Nuclei were sorted as described (Lister et al., 2013). Double positive nuclei were retained for RNA-seq, ChIP-seq, and MethylC-seq library preparation and sequencing.

Western blot

Request a detailed protocol

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× g, supernatants were preserved and protein concentration was determined by the BCA method (Pierce). Protein bands were separated in 8% PAGE gels and transferred to nitrocellulose membranes. After blocking in TBS-tween with 5% milk, DNMT3A was detected by the use of anti-DNMT3A antibody (Abcam) and chemiluminescence. DNMT3A bands were normalized to ACTIN content in each sample.

Patch-clamp electrophysiology

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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 Dnmt3a knockout and control samples, we grouped SQuIRE’s TE subfamily expression estimates (FPKM) into their respective families and classes and performed two-sample t-tests.

Enrichment test of GO terms in DE genes

Request a detailed protocol

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

Request a detailed protocol

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 Dnmt3a cKO and the control animals to call SNPs against the mouse mm10 genome. We followed the GATK ‘best practices for germline SNPs and indels in whole genomes and exomes’ pipeline (DePristo et al., 2011; McKenna et al., 2010; Van der Auwera et al., 2013). Briefly, raw data were first trimmed to remove sequencing adapters and low-quality sequences (minimum Phred score 20) using Trim Galore in the paired-end mode. Clean data were then mapped to the mm10 genome using BWA (Li and Durbin, 2009) (v0.7.13-r1126). Duplicates reads were marked with Picard (Broad Institute, 2018). Then the analysis-ready reads were fed into GATK (v3.7) to perform two rounds of joint genotyping and base recalibration. Variants were then filtered using the following criteria: QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < –12.5 || ReadPosRankSum < –8.0 || SOR > 4.0. By that, we identified 548,530 and 507,669 SNPs (relative to mm10) in the Dnmt3a cKO and the control animals, respectively. At last, we created a substituted genome to mask out all these SNPs (replaced with Ns) with the ‘maskfasta’ tool in the BEDTools suite (Quinlan and Hall, 2010) (v2.27.1). This substituted genome was used in the following MethylC-seq data processing pipeline.

MethylC-seq data processing

Request a detailed protocol

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:

%mC = 100×mh,

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:

%mC_adj=%mc%NCR100%NCR,where %mC_adj[0, 100].

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

Request a detailed protocol

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 Dnmt3a cKO hypo-DMRs in P39 (Supplementary file 3). Note that with these criteria we also found 89 Dnmt3a cKO hyper-DMRs, which we thought were noise and/or SNPs that failed to be detected by the masking pipeline described earlier. Therefore, we removed these hyper-DMRs from further consideration.

Enrichment test of DMRs and other genomic regions

Request a detailed protocol

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:

log2OR = log2n11 + log2n12 log2n21 log2n22
se(log2OR) = 1/n11 + 1/n12  + 1/n21+1/n22
ci(log2OR) = [log2OR  2×se(log2OR),  log2OR + 2×se(log2OR)]

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

Request a detailed protocol

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

Request a detailed protocol

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 Dnmt3a cKO DMRs. The parameters used are as follows: ‘-size 500 -len 8,10,12S 25 -fdr 100p 10 -mset vertebrates -bits -gc -nlen 3 -nomotif’. A set of non-neural DMRs (Hon et al., 2013) was used as the background.

ChIP-seq data processing

Request a detailed protocol

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: ‘--bin-size 200 --gaps-allowed 1’. For H3K27me3, we used the following parameters: ‘--bin-size 200 --gaps-allowed 3’. The IgG sample was used as a control.

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

Request a detailed protocol

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

Request a detailed protocol

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 Dnmt3a cKO and the control samples, we further grouped these DMVs into six categories, namely consistent (exact same DMV in the two conditions), overhang, cKO unique, control unique, expanded (wider in the cKO), and shrunken (wider in the control) (see illustrations in Figure 7A).

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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/.

Appendix 1

Appendix 1—key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Gene (Mus musculus)Dnmt3aGENCODEGENCODE:ENSMUSG00000020661GENCODE vM10
Genetic reagent (Mus musculus)Dnmt3a-floxedPMID:10555141Backcrossed to C57BL/6 for at least seven generations
Genetic reagent (Mus musculus)Nex-Cre; Neurod6-CrePMID:17146780Jackson Laboratories
Backcrossed to C57BL/6J for 10 generations
Genetic reagent (Mus musculus)INTACTPMID:26087164(B6.129-
Gt(ROSA)26Sor
tm5(CAG-Sun1/sfGFP)Nat/MmbeJ)
Genetic reagent (Mus musculus)Dnmt3a-KO; Dnmt3a cKOThis paperSee ‘Materials and Methods’, section ‘Generation of the Dnmt3a cKO mice line’
AntibodyAnti-NeuN (Mouse monoclonal, Clone A60)MilliporeMAB377(1:1000)
AntibodyAnti-DNMT3A (Rabbit polyclonal)AbcamAb2850(1:250)
AntibodyAnti-H3K27me3 (Rabbit polyclonal)Active MotifCat:#39156; RRID:AB_2636821(5 µl)
AntibodyAnti-H3K27ac (Rabbit polyclonal)Active MotifCat:#39133; RRID:AB_2561016(5 µl)
AntibodyAnti-H3K4me3 (Rabbit monoclonal)Millipore SigmaCat:#04–745(5 µl)
AntibodyAnti-IgG (Mouse unknown clonality)Jackson ImmunoResearch LabsCat#015-000-003; RRID:AB_2337188(2 µl)
Sequence-based reagentCytosine-methylated adaptersIlluminaAD001, AD005
Sequence-based reagentTruSeq PCR Primer MixIllumina20015960, 20015961
Peptide, recombinant proteinT4 DNA ligaseNew England BiolabsM0202L
Commercial assay or kitSingle-Cell RNA Purification KitNorgenCat:#51800
Commercial assay or kitOvation RNA-Seq System V2NuGENCat:#7102–32
Commercial assay or kitKAPA LTP Library Preparation KitRocheCat:#KK8232
Commercial assay or kitQiagen DNeasy Blood and Tissue kitQiagenCat:#69504
Commercial assay or kitEZ methylation Direct kitZymoCat:#D5021
Commercial assay or kitAccel-NGS 2S Plus DNA Library KitSwift BiosciencesCat:#21024
Chemical compound, drugProtease inhibitorSigma-AldrichCat:#11836153001
Chemical compound, drugPhosphatase inhibitorPierceCat:#A32957
Chemical compound, drugKapa HiFi Hotstart ReadymixKapa Biosystems07958935001-KK2602
Software, algorithmMinianalysisMinianalysis (https://www.synaptosoft.com/MiniAnalysis/)RRID:SCR_002184
Software, algorithmpCLAMPMolecular Devices (https://www.moleculardevices.com/products/software/pclamp.html)RRID:SCR_011323
Software, algorithmReconstructSynapse Web (https://synapses.clm.utexas.edu/tools/reconstruct/reconstruct.stm)RRID:SCR_002716
Software, algorithmActivity MonitorMED Associates (https://www.med-associates.com/product/activity-monitor/)RRID:SCR_014296
Software, algorithmANY-mazeSan Diego Instruments (https://sandiegoinstruments.com/product/any-maze/)RRID:SCR_014289
Software, algorithmFastQCBabraham Bioinformatics (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)RRID:SCR_014583v0.11.8
Software, algorithmTrim GaloreBabraham Bioinformatics (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/)RRID:SCR_011847v0.5.0
Software, algorithmCutadaptDOI:10.14806/ej.17.1.200RRID:SCR_011841v1.16
Software, algorithmSTARPMID:23104886RRID:SCR_004463v2.5.1b
Software, algorithmRSEMPMID:21816040RRID:SCR_013027v1.2.30
Software, algorithmedgeRPMID:19910308RRID:SCR_012802v3.24.1
Software, algorithmSQuIREPMID:30624635
Software, algorithmclusterProfilerPMID:22455463RRID:SCR_016884v3.10.0
Software, algorithmGATKPMID:21478889RRID:SCR_001876v3.7
Software, algorithmBWAPMID:19451168RRID:SCR_010910v0.7.13-r1126
Software, algorithmPicardBroad (https://broadinstitute.github.io/picard/)RRID:SCR_006525V2.10.6
Software, algorithmBEDToolsPMID:20110278RRID:SCR_006646v2.27.1
Software, algorithmmethylpymethylpy (https://github.com/yupenghe/methylpy)v1.3.2
Software, algorithmbowtie2PMID:22388286RRID:SCR_016368v2.3.0
Software, algorithmMatchItDOI:10.18637/jss.v042.i08v3.0.2
Software, algorithmBARTPMID:29608647
Software, algorithmHOMERPMID:20513432RRID:SCR_010881v4.8.3
Software, algorithmENCODE Transcription Factor and Histone ChIP-Seq processing pipelinePMID:22955991RRID:SCR_021323v1.1.6
Software, algorithmsamtoolsPMID:19505943RRID:SCR_002105v1.2
Software, algorithmSICERPMID:19505939RRID:SCR_010843
Software, algorithmDiffBindPMID:22217937RRID:SCR_012918v2.10.0
Software, algorithmDESeq2PMID:25516281RRID:SCR_015687v1.22.1
Software, algorithmGREATPMID:20436461RRID:SCR_005807v3.0.0
Software, algorithmMethylSeekRPMID:23828043RRID:SCR_006513v1.22.0
Software, algorithmregioneRPMID:26424858V1.14.0
Software, algorithmRR Project for Statistical Computing (https://www.r-project.org/)RRID:SCR_001905v3.5.0
Software, algorithmMATLABMathWorks (https://www.mathworks.com/products/matlab.html)RRID:SCR_001622v2017a
Software, algorithmPythonPython Programming Language (https://www.python.org/)RRID:SCR_008394v3.x

Data availability

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/.

The following data sets were generated
The following previously published data sets were used

References

Decision letter

  1. Anne E West
    Reviewing Editor; Duke University, United States
  2. Lu Chen
    Senior Editor; Stanford University, United States
  3. Alvaro Rada-Iglesias
    Reviewer; University of Cologne, Germany

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Dnmt3a knockout impairs synapse maturation and is partly compensated by repressive histone modification H3K27me3" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by a Reviewing Editor and Lu Chen as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Alvaro Rada-Iglesias (Reviewer #3).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1) The reviewers all felt that claims regarding the functional compensation by H3K27me3 for the loss of Dnmt3a were too strongly stated. The reviewers discussed this point extensively and decided that both the title and the sections of the paper that describe and discuss these data should be changed to reflect that the replacement of the marks has been observed but that the functional relevance of this change has not been directly tested. The reviewers appreciate that functional testing of this hypothesis is well beyond the scope of the current manuscript – only text revisions are required.

2) The reviewers strongly suggest (but do not absolutely require) that the authors perform bioinformatic analysis of the H3K27me3 compensation against gene expression in the Dnmt3a cKO in order to provide experimental evidence for the claim of functional compensation. One of the reviewers also offered a very thorough list of other possible analyses (especially regarding studies of enhancer marks) that would strengthen the study, however these are entirely at the discretion of the authors to consider.

3) The reviewers have described below specific examples where addition of further technical details would strengthen the manuscript and we ask the authors to expand on these experimental points in the text to help the readers evaluate the study. In particular all of the reviewers expressed some concern about the use of only two replicates, and we had a long discussion of whether a minimum of three should be required. It was not entirely clear in every case whether each replicate was obtained from only a single mouse or whether the samples might be pooled (which would strengthen the argument for using only 2). In either case, in the end we felt the manuscript was strong enough that we did not feel it was necessary for the authors to add more replicates, but we would like the authors to address in the text any limitations of the study that readers may need to keep in mind with the use of duplicates as opposed to larger replicate numbers.

Reviewer #1 (Recommendations for the authors):

1) Because of the way the authors present the examples and the comparisons, it is unclear if the compensatory upregulation of H3K27me3 occurs de novo, or whether regions that are previously methylated with H3K27me3 are the regions where this compensation occurs. The latter could suggest this might be a consequence of having more substrate for histone methylation – and may not actually serve any function. This is particularly important to delineate as there are several other DMV clusters that are resistant to changes in gene expression that do not have this compensatory upregulation of H3K27me3 – like clusters 4, 5 and 6 in Figure 8.

2) It would be more narratively and biologically significant if the authors could connect their data analysis and observations about compensatory H3K27me3 upregulation and changes to gene expression in the context of the developing cortical neurons, or the behavioral phenotypes observed in the introductory figures.

3) The inferences made in Figure 8 are a little inconsistent with those described earlier. For example, in Figure 5A the authors describe the overrepresentation of predicted PRC2-TF binding motifs at upregulated and downregulated genes during Dnmt3a-cko, while later discussions move toward describing PRC2-mediated regulation as the reason for resistance to changes in gene expression due to the cko. Additionally, H3K27me3 peaks in Figure 5B do not appear to be all that different between the control and cko mouse at Mab21l2 – diluting the punchline of the figure.

4) Figure 8 is difficult to read.

Reviewer #2 (Recommendations for the authors):

As described in the public comments, this study provides interesting new insights into organismal and cellular consequences of DNMT3A conditional knockout in excitatory neurons and presents novel effects on H3K27me3 in these mice. I have several concerns I think would be important to address before publication in eLife however. Here I elaborate in specific detail:

1. The analysis and core conclusions about changes in histone modifications in the study are based on n=2 experiments. While it is potentially beyond the scope of a revision to repeat additional replicates for all analyses that are carried out, it would substantially strengthen the study to perform additional replicates of H3K27me3 ChIP in order to rigorously show that the core finding regarding H3K27me3 changes is robust to sample-to-sample technical variability. This is particularly important given the subtle magnitude of changes that are described.

2. Though the authors provide strong characterization of developmental changes of H3K27me3 and show changes in this mark in response to loss of DNMT3A, their analyses fall short of definitively showing that this is a true functionally compensatory effect. Ideally, an experimental analysis would be done to test this by blocking H3K27me3 build up in DNMT3A mutants and determining if de-compensated gene expression changes occur. Given that this is outside the scope of a revision however, additional computational approaches could be employed to better test the hypothesis:

a. Figures 6E and S10A present the analysis that is most relevant to whether accumulation of K27me3 in the cKO impacts DEG and therefore may play a compensatory role. However, additional analysis could more decisively test the hypothesis. For example, do non-dysregulated genes containing DMRs show less K27me3 accumulation than down-regulated genes, but more K27me3 accumulation than upregulated genes? The authors could include a boxplot of a well-matched set of non-DE genes with overlapping DMRs in panel 6E to test this. A more powerful analysis, however, might be to examine a large set of genes from the genome that are selected to contain similar numbers of DMRs (in order to control for relative potential effects of demethylation) and ask if changes in K27me3 correspond with up-regulation (low K27me3 accumulation), no dysregulation (intermediate accumulation) or downregulation (high K27me3 accumulation).

b. On p19 the authors state that "P39 cKO DMRs were significantly enriched in peaks and DM regions of H3K27me3 (28.1% overlapped H3K27me3 peaks in cKO, p<1e-100) (Figure 6B)." However, 6B only presents the overlap with peaks in each genotype as an odds ratio and doesn't make it clear how the enrichment of DMRs in DM regions was assessed. This appears to be assessed to some degree in S8F (note the S8 figure legend appears to be wrong for this panel), but no analysis is presented to control for the general tendency of DMRs and H3K27me3 peaks to overlap. It seems imperative, given that DMRs already overlap with H3K27me3 peaks, that a resampling approach is used comparing DM peaks to non-DM K27me3 peaks that are matched to have similar wild-type H3K27me3 signal and size as DM peaks. Likewise, what would the changes in H3K27me3 signal presented in 6C look like at non-DMR regions that contain similar baseline wild-type H3K27me3 signal? Without these analyses it is not clear that DMRs are particularly unique in showing increased K27me3 signal.

c. The increase in H3K27me3 in the cKO at the mab21l2 locus at P0 (figure 5B), before methylation has or has not been established, suggests that H3K27me3 may change before methylation differences, rather than in reaction to changes in methylation. Is this a common effect seen across other genes? How does this impact the idea that compensation occurs due to lack of methylation? Related to this issue, on p.19 no specific plot is referenced to support the statement that "There was no corresponding increase of H3K27me3 at these DMRs in newborn (P0) or fetal (E14) neurons." A quantitative analysis as shown in 6C but for the E14 and P0 time-points could assess this directly.

d. Why were developmental DMs called between E14 and P39 while DMRs were called between P0 and P39? Developmental comparisons would be more justified if both used P0 and P39.

e. In the absence of additional data supporting H3K27me3 accumulation as a functional compensation mechanism for loss of DNA methylation, the title of Figure 6 should not so strongly claim that H3K27me3 changes compensate for mC loss.

3. While initially assessed by the authors, mCH changes are not integrated into the analysis of H3K27me3 changes. Previous studies by this group and others have suggested the importance of mCH in gene regulation within specific neuron populations and, given the interest to the field, more thorough analysis is warranted regarding this neuronal methylation and H3K27me3.

a. Given the association, even if limited, between differential mCH loss and gene dysregulation in the cKO, it is relevant to examine how changes in K27me3 relate to, or compensate for, loss of mCH. For example, do genes that are not dysregulated in the cKO, and normally have high mCH (which is lost in the cKO) show higher increases in K27me3 compared to non-dysregulated genes that have low mCH (and therefore lose less mCH)? Reciprocally, one could ask if a population of genes selected to have similar, high mCH loss in the cKO, show differential gene dysregulation as a function of differential K27me3 accumulation.

b. What are the changes in mCH at DM H3K27me3 sites? Since mCG and mCH loss are correlated, one might expect that these are sites that normally also contain high mCH. It might not be feasible to call mCH DMRs given the technical limitations, but quantification of this signal can be done at DMs and DMRs, and the potential for it to contribute to these effects could be discussed.

4. The authors present evidence of PRC2 binding at DEGs (Figure 5A) as a motivation to assess H3K27me3 in their mice. However, they then observe changes in H3K27me3 that they contend do not occur under normal conditions and propose that these changes are an underlying mechanism that blocks gene dysregulation in these mice. How are we to reconcile this original rationale of a normal function for PRC2 at dysregulated genes with the model of de novo compensatory H3K27me3 accumulation? Some explanation by the authors of how these separate notions may relate to one another would be helpful.

5. Some effects on H3K27ac are observed at DMRs, but the authors ignore the prospect that these effects could drive changes in gene expression. Previous groups have shown significant changes in H3K27ac in DNMT3A models, and some of the data the authors present suggest that this could also explain some of the observations in their model.

a. While DM K27ac peaks are not identified in the cKO, in figure 6B there is a similar enrichment of DMRs to associate with K27ac peaks in the cKO vs control as there is for K27me3. In addition, there is a similar increase in odds ratio of enrichment of K27ac peaks to the of K27me3 at DMRs for the cKO in figure 6B. The analysis in S8F-G that is referenced when referring to K27ac is either missing (S8G) or is not showing data for K27ac (S8F). The authors should at least comment in the Results section that these signals exist and may be significant for gene regulation.

b. Are two biological replicates sensitive enough to call differences in H3K27ac? In the context of a small number of replicates and the subtle magnitude of changes in gene expression that are observed, negative results for differential peaks should be interpreted with caution. In contrast, the larger peaks observed in K27me3 may have led to increased ability to call differential peaks.

6. Though the characterization of DNMT3A cKO in excitatory neurons is novel, if the authors are claiming that these neurons specifically are driving behavioral phenotypes then further contextualization and characterization of the extent of conditional deletion is needed.

a. It would be useful to justify the rationale for investigating excitatory neurons, and how using this class of neurons can uncover novel findings. Existing studies of DNMT3A cKO have shown behavioral and epigenomic effects in pan-neuronal and inhibitory neuron conditional knockouts. Comparisons between these models and discussion of differences could offer strong insights into what unique circuits may be most disrupted by DNMT3A cKO in excitatory neurons.

b. If the authors are claiming that behaviors and phenotypes are the result of specifically knocking out DNMT3A in all excitatory neurons, further evidence should be provided that the Cre-line is being expressed sufficiently and specifically in excitatory neurons, and that their genomic data is faithfully capturing the majority of excitatory cells while excluding non-excitatory cells.

7. The manuscript would be further improved by clarification of some methods and statistics, and modifications in data presentation:

a. The authors should provide information on the statistics performed on the behavior data. Are there significant differences between sexes, or significant sex by genotype interactions for these behavioral tasks? A two-way ANOVA using genotype and sex would be useful in understanding the impact of these variables. There appears to be many effects that are sex-specific, yet this appears to be overlooked in the text. Furthermore, if there are robust sex-specific effects then the subsequent molecular and genomic analyses should be done using both sexes, and comparisons could be done to explain the molecular underpinnings of these differences.

b. The authors might also consider showing preference index when analyzing social approach data when comparing genotypes.

c. What is the background set of genes for 2B-C? Synaptic genes are enriched in both up- and down-regulated gene lists. How is gene expression normalized? Is there a discovery bias for more highly expressed genes?

d. What value was used to identify equally expressed unchanged genes in figure 4B, TPM? Identifying equally expressed genes can be complicated for nuclear RNA where unprocessed RNAs contribute differentially to total counts based on the length of the gene. In addition, previous analysis of mCH enrichment at dysregulated genes in inhibitory neuron-specific KO of DNMT3A also found DEGs to be differentially methylated based on direction of change but showed that all other genes were more lowly methylated than these two groups (Lavery et al. 2020). Here the authors see the same relationship between DEG groups but find other genes to be higher in mCH rather than lower. Do the authors have any perspective on why these differences are observed?

e. On p22 the authors state "We found that developmental peaks (DevLoss and DevGain) were relatively unaffected by the cKO(H3K27me3 = 0.11, -0.20 respectively, in units of log10(CPM+1)), whereas Group cKO had ten-fold larger mean effect (H3K27me3 = 1.12 respectively)." Given the y scales in figures 7E and S10H it is unclear how they derived a ∆ value of 1.12 on a log 10 scale for Group cKO (e.g. in S10H the median values look much less 1 unit different on a log10 scale). Please clarify the actual magnitude of effect sizes.

f. The authors describe the association between gene body mCH and gene expression changes as only accounting for 1% of the variance. It would be useful to have a point of reference to understand if this is a low value or a high value in comparison to other potential predictors of gene expression changes. For example, how much does the proposed overcompensation by H3K27me3 accumulation contribute to the genome-wide variance in gene expression? In addition, is it the case that variance in the measured fold-changes will be mitigated by addition of more replicates to more accurately estimate true fold-changes? If this is the case, it should be noted.

g. Using different color schemes would be helpful in 8A and 8B to avoid confusion.

h. How was clustering analysis done in 8D? What do the values for K27me3 Diff. Mod and DEG enrichment indicate (i.e. what does the scale from 0 to 0.02 represent)?

Reviewer #3 (Recommendations for the authors):

– The authors delete exon 19 of Dnmt3a. I think they should provide a more thorough description of what this deletion is expected to cause: loss of function due to premature stop codon, loss of enzymatic activity, etc. Moreover, it would be desirable to show by immunostaining whether DNMT3A protein is specifically lost in excitatory neurons.

– Regarding the DMRs identified in DNMT3A KO neurons, the authors should use the generated profiles for H3K27ac in P0 and P39 WT neurons to determine which DMRs overlap enhancers that are active in any of these two neuronal stages. This is more relevant than overlaps with enhancers that could be active in any given cell type. Most importantly, the authors should evaluate the impact the loss of DNA methylation can have on enhancer activity by measuring H3K27ac levels and eRNAs in DNMT3A KO neurons. Measuring eRNA levels can be particularly relevant, as recent evidences suggest that is a better marker of enhancer activity than H3K27ac.

– The authors should aim at performing experiments that address the proposed compensatory role for PRC2 and H3K27me3. I understand that generating a KO for PRC2 could be rather complicated, but perhaps the authors could consider using new drugs that act as very specific inhibitors against the enzymatic activity of PRC2.

https://doi.org/10.7554/eLife.66909.sa1

Author response

Essential revisions:

1) The reviewers all felt that claims regarding the functional compensation by H3K27me3 for the loss of Dnmt3a were too strongly stated. The reviewers discussed this point extensively and decided that both the title and the sections of the paper that describe and discuss these data should be changed to reflect that the replacement of the marks has been observed but that the functional relevance of this change has not been directly tested. The reviewers appreciate that functional testing of this hypothesis is well beyond the scope of the current manuscript – only text revisions are required.

We agree with the Reviewers that it is important to clearly define the conclusions we can draw from our data, as well as their limitations. Indeed, our data demonstrate that loss of Dnmt3a-dependent gain of methylation during the perinatal period leads to specific changes in H3K27me3, but they do not fully clarify the functional significance of these processes or show that Polycomb associated marks fully compensate for reduced DNA methylation.

To clarify this, we have revised the title of the paper: Dnmt3a knockout in excitatory neurons impairs postnatal synapse maturation and increases the repressive histone modification, H3K27me3

We have further removed the following phrase from the Abstract: “partially compensating for the loss of DNA methylation.” We have limited to the Discussion section our comments about a potential compensatory role for the gain of H3K27me3.

2) The reviewers strongly suggest (but do not absolutely require) that the authors perform bioinformatic analysis of the H3K27me3 compensation against gene expression in the Dnmt3a cKO in order to provide experimental evidence for the claim of functional compensation. One of the reviewers also offered a very thorough list of other possible analyses (especially regarding studies of enhancer marks) that would strengthen the study, however these are entirely at the discretion of the authors to consider.

We appreciate the suggestions from the reviewers. We have added new analyses to further support our claim, including comparisons of the changes in H3K27me3 in different groups of genes (shown in Figure 5E and Figure 5 — figure supplement 2A). These are fully described below, in the response to Reviewer #2 (comment 2a).

3) The reviewers have described below specific examples where addition of further technical details would strengthen the manuscript and we ask the authors to expand on these experimental points in the text to help the readers evaluate the study. In particular all of the reviewers expressed some concern about the use of only two replicates, and we had a long discussion of whether a minimum of three should be required. It was not entirely clear in every case whether each replicate was obtained from only a single mouse or whether the samples might be pooled (which would strengthen the argument for using only 2). In either case, in the end we felt the manuscript was strong enough that we did not feel it was necessary for the authors to add more replicates, but we would like the authors to address in the text any limitations of the study that readers may need to keep in mind with the use of duplicates as opposed to larger replicate numbers.

We appreciate this concern, and we agree that it is important to show that findings are replicable across multiple independent samples.

When tracing back the animal IDs to create a table of the number of animals included in each sample, we discovered an error in the dataset. Although we intended to use independent pools of mice for each biological replicate, we found that due to an experimental tracking error some of the replicates contained tissue from overlapping pools of animals. In particular, the P39 Dnmt3a cKO group (pooled from 3 mice) was split into two technical replicates but mislabeled as two biological replicates and used as such in the RNA-seq and MethylC-seq analyses. For this reason, for this resubmission we decided to repeat our experiments and to generate new RNA-seq data and MethylC-seq data for both P39 Dnmt3a cKO and control mice. Each condition now has two independent, newly-generated biological replicates, each pooled from two mice. We have added the information about the number of animals contributing tissue to each sample (including P0 MethylC-seq data and all ChIP-seq data which are not re-sequenced) in the “num_pooled_animals” column in Supplementary File 1. Most of our sequencing data come from tissue samples pooled from two mice, with the exception of P0 MethylC-seq data. For that time point, we used 6 control and 2 cKO samples that each came from a single mouse.

To show that the new dataset largely agrees with the previously collected data, we plotted the correlation across replicates from the two batches computed using mCG and mCH in 10kb genomic bins for MethylC-seq data, and gene expression for RNA-seq data (see Author response image 1). The molecular profiles are largely consistent across the two batches (Spearman correlation > 0.80; note that mCH in Dnmt3a cKO samples are essentially absent, so those samples were excluded from the mCH analysis). Moreover, the consistency across biological replicates in the new batch is also high (r > 0.93). We have also added plots in the manuscript to demonstrate the consistency across the replicates for RNA-seq (Figure 3 — figure supplement 1A), MethylC-seq (Figure 3 — figure supplement 3A), and ChIP-seq (Figure 6 — figure supplement 2).

Author response image 1
Heatmap to show the Spearman correlations across the control and Dnmt3a cKO samples from the two batches of data (new vs old). The correlations coefficients were computed using CG methylation levels in 10Kb genomic bins (left), CH methylation levels in 10Kb genomic bins (middle), and gene expression (log10TPM, right).

As is often the case for RNA-seq data, we found differences between the new and old batches (i.e. batch effects) that cannot be ignored (see Author response image 2A). Due to the mislabeling of technical replicates as biological replicates in our previous analysis, our original submission likely had inflated statistical significance for detecting differentially expressed genes. Using the new batch of RNA-seq data that includes true biological replicates, we detected fewer significant DE genes under the same FDR < 0.05 criteria (70 genes compared to the previous 1720 genes). However, the fold-changes of these DE genes are in general concordant across the two batches (see Author response image 2B). This new dataset and analysis replicate and validate the disruption of transcriptome we observed earlier. The smaller number of differentially expressed genes we now report is consistent with the fact that statistical significance was inflated in our previous analysis by the non-independence of the replicate samples.

Author response image 2
(A) PCA plot using log10TPM gene expression in the two batches of RNA-seq samples. (B) Scatter to show the consistency of gene expression fold-changes (Dnmt3a cKO vs. control) across the two batches of RNA-seq samples using significant DE genes detected in the new batch (left) and significant DE genes detected in the old batch (right).

Confidence in our differential expression findings is also supported by the strong correspondence we now find between genes that lose mCH and those that increase mRNA expression (Figure 3D-E). This gives additional confidence that the differential expression analysis identifies reliable and consistent effects in the Dnmt3a cKO.

We apologize for the confusion caused by the replacement of these new sequencing datasets. We hope that we have demonstrated that the new data largely agrees with the old, and that most of our conclusions are still supported by the new data albeit with reduced statistical power. In the revised manuscript, all results generated using P39 RNA-seq and MethylC-seq data were updated with the new datasets. Nevertheless, some analysis is no longer feasible due to fewer detected DE genes, and we have modified the manuscript accordingly. In the following responses, we address each reviewers’ specific concerns.

Reviewer #1 (Recommendations for the authors):

1) Because of the way the authors present the examples and the comparisons, it is unclear if the compensatory upregulation of H3K27me3 occurs de novo, or whether regions that are previously methylated with H3K27me3 are the regions where this compensation occurs. The latter could suggest this might be a consequence of having more substrate for histone methylation – and may not actually serve any function. This is particularly important to delineate as there are several other DMV clusters that are resistant to changes in gene expression that do not have this compensatory upregulation of H3K27me3 – like clusters 4, 5 and 6 in Figure 8.

We appreciate the reviewer’s suggestion and apologize for the confusion. To address this question, we have added the following analysis (Figure 4 — figure supplement 3A) to compare the baseline H3K27me3 levels in the differentially modified regions with their shuffle controls:

Here two types of random shuffle controls were used. In the first shuffle control (the middle two violins), we randomly shuffled the H3K27me3 DM regions in the whole genome, keeping the peak sizes and the chromosomes where the peaks reside in same as the original peaks. In the second set (the rightmost two violins) we required the same constraints as the first set but further restricted the shuffles must be within the H3K27me3 peaks (the union of control and cKO peaks). We can see that the regions with up-regulation of H3K27me3 (the leftmost two violins), had an intermediate level of H3K27me3 (median ~8 RPKM) in the control samples to begin with, which is higher than the level in random shuffles in the genome but lower than the level in random shuffles in H3K27me3 peaks. Therefore, these regions with upregulation of H3K27me3 after Dnmt3a cKO include regions that were previously methylated with H3K27me3 in the control samples, and we presume such a medium but non-zero level of H3K27me3 is fine-tunable and hence can act as a “failsafe” triggered after the original repression mechanism, DNA methylation, was disrupted. Of course, it is important to note that our ChIP-seq data come from a heterogeneous pool of pyramidal neurons, and it is possible that the quantitative increase in the H3K27me3 level corresponds to the genuine gain of de novo histone methylation in a particular subpopulation of excitatory cells.

We can also reach a similar conclusion in our updated DMV analysis (now in Figure 7C). C1 clusters are DMVs that highly overlap with the H3K27me3 DM regions, and they have an intermediate level of H3K27me3 (~8 RPKM) that is not as high as those in cluster 6 (~20 RPKM) but still higher than those in clusters 5, 7 and 8 (~0 RPKM). As a comparison, regions in cluster 6 were under tight regulation with high levels of H3K27me3 and a low level of DNA methylation (~10% mCG); genes within these DMVs were mostly not expressed. Therefore, losing the small amount of DNA methylation in these regions did not alter gene expression. As for the regions in clusters 5, 7 and 8, they are active with both low DNA methylation and H3K27me3, and hence they are not affected by the cKO.

We have added these results to the Result section of the manuscript:

“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 Dnmt3a cKO.”

2) It would be more narratively and biologically significant if the authors could connect their data analysis and observations about compensatory H3K27me3 upregulation and changes to gene expression in the context of the developing cortical neurons, or the behavioral phenotypes observed in the introductory figures.

We agree that linking the effects of the cKO on adult neurons with developmental patterns of epigenetic regulation would be valuable. Our study includes developmental ChIP-seq profiles of H3K27me3, including three time points (E14, P0 and P39; see Figure 5B-C). These data show that the majority of the changes in H3K27me3 emerge in the adult neurons (P39) and are not present at earlier developmental time points (Figure 4 — figure supplement 1D). Moreover, we related the patterns of regulations at DMVs with their developmental history of H3K27me3 (Figure 7).

We have now added a new analysis to address the relationship between upregulated H3K27me3 signal and gene expression changes at P39 (see Figure 5 — figure supplement 2C and also our replies to reviewer #2, comment 3a).

Unfortunately, we did not profile gene expression after cKO in development so we can’t examine these observations in the context of development. Although we agree that this would be of interest, we believe that this question would be better addressed using single-cell transcriptomics (scRNA-seq) to disentangle the contribution from diverse subtypes of excitatory pyramidal neurons. We, therefore, plan to apply single-cell methods in the future to address this question.

3) The inferences made in Figure 8 are a little inconsistent with those described earlier. For example, in Figure 5A the authors describe the overrepresentation of predicted PRC2-TF binding motifs at upregulated and downregulated genes during Dnmt3a-cko, while later discussions move toward describing PRC2-mediated regulation as the reason for resistance to changes in gene expression due to the cko.

The BART results for up-regulated and down-regulated DE genes (now in Figure 4A) are bioinformatics predictions that motivated us to look into Polycomb repression via H3K27me3 profiling. With this rationale, we then performed the ChIP-seq analysis of H3K27me3 in the control and cKO neurons. We agree with the reviewer that there is not a clear explanation for the enrichment of PRC2 related sequence motifs at both up- and down-regulated genes. In light of this, and in response to other reviewer comments, we have significantly softened our interpretation of the direct link (“compensation”) between loss of DNA methylation and increased H3K27me3. Instead, we now report that increased H3K27me3 occurs preferentially in regions that lose DNA methylation, but we do not draw a causal link or interpret these data in terms of direct functional compensation.

Additionally, H3K27me3 peaks in Figure 5B do not appear to be all that different between the control and cko mouse at Mab21l2 – diluting the punchline of the figure.

Regarding the size of the difference in H3K27me3 at the Mab21l2 locus, we appreciate that it can be hard to judge the magnitude of the differential enrichment of the chromatin modification (due to the normalized scale to make it comparable across E14, P0, and P39) in the browser view (now in Figure 4B). The quantitative analysis of differential binding at this locus showed that it was highly significant (tested with DEseq2, fold-change = 2.41, FDR=1.33e-4). We have added a shaded box in Figure 4B to highlight the H3K27me3 signal differences between the P39 control and cKO mouse at the Mab21l2 locus, and also a new plot in Figure 4C to show the quantification of the increase in H3K27me3 ChIP-seq signal in each replicate at this locus.

We have also added analysis to show the increases of H3K27me3 were also enriched in non-DE genes with the largest magnitude of mCH (see Figure 5 — figure supplement 2B and also our replies to reviewer #2, comment 3a).

4) Figure 8 is difficult to read.

We have re-organized the figure (now Figure 7) to make it easier to read.

Reviewer #2 (Recommendations for the authors):

As described in the public comments, this study provides interesting new insights into organismal and cellular consequences of DNMT3A conditional knockout in excitatory neurons and presents novel effects on H3K27me3 in these mice. I have several concerns I think would be important to address before publication in eLife however. Here I elaborate in specific detail:

1. The analysis and core conclusions about changes in histone modifications in the study are based on n=2 experiments. While it is potentially beyond the scope of a revision to repeat additional replicates for all analyses that are carried out, it would substantially strengthen the study to perform additional replicates of H3K27me3 ChIP in order to rigorously show that the core finding regarding H3K27me3 changes is robust to sample-to-sample technical variability. This is particularly important given the subtle magnitude of changes that are described.

We understand the reviewer’s concern regarding the limited number of biological replicates. To reduce variability due to individual differences, each of our ChIP-seq samples included pooled tissue from 2 mice. This information is now included in the

“num_pooled_animals” column in Supplementary File 1 and explained in the main text:

“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.”

We have added additional analyses showing the consistency of the ChIP-seq signal in genomic bins and peak regions across biological replicates. These results are now shown in Figure 6 — figure supplement 2 and we have added descriptions of these results in the manuscript as follows:

“We also examined differentially modified (DM) regions at P0 vs. E14 and P0 vs. P39 (Figure 6 — figure supplement 1A). However, due to greater biological variability at the perinatal time point, the two replicate ChIP-seq datasets at P0 were less consistent than the E14 and P39 samples (Figure 6 — figure supplement 2). As a result, our data were not well powered to detect changes at P0 and we chose to focus on the E14 vs. P39 DM regions as sites of developmental chromatin remodeling.”

2. Though the authors provide strong characterization of developmental changes of H3K27me3 and show changes in this mark in response to loss of DNMT3A, their analyses fall short of definitively showing that this is a true functionally compensatory effect. Ideally, an experimental analysis would be done to test this by blocking H3K27me3 build up in DNMT3A mutants and determining if de-compensated gene expression changes occur. Given that this is outside the scope of a revision however, additional computational approaches could be employed to better test the hypothesis:

We agree that additional experiments that block H3K27me3 in Dnmt3a cKO animals are required to establish a causal role for the compensating effect. However, as noted by the Reviewer, such experiments are beyond the scope of this study. Instead, we have performed additional bioinformatic analyses to better test the hypothesis according to the reviewer’s suggestion.

As suggested by the reviewer, we used box plots to compare the gain of H3K27me3 in different groups of genes. These analyses (now in Figure 5E and Figure 5 — figure supplement 2A) compare up- and down-regulated genes with the control gene set chosen to have matching expression levels in the control samples. Moreover, we separately analyzed genes that overlap DMRs and those which do not overlap DMRs; the corresponding control genes were likewise grouped by whether or not they overlap DMRs. (See Methods and Figure 3 — figure supplement 4B on how we selected the expression-matched non-DE genes).

This analysis shows that down-regulated DE genes with overlapping DMRs (the first box in Figure 5E, upper panel) and non-DE genes selected with matched expressions of up-regulated DE genes (the fourth box in Figure 5E, upper panel) show a small but significant median increase of H3K27me3 in the cKO (Wilcoxon rank-sum test; asterisks in the middle of each boxplot denotes p-value: *, p < 0.05). When we considered a larger set of DE genes with a more relaxed threshold (FDR<0.2, Figure 5 — figure supplement 2A), we observed that down-regulated DE genes containing DMRs accumulate significantly more H3K27me3 than non-DE genes containing DMRs. By contrast, up-regulated genes had no significant accumulation of H3K27me3 and were not significantly different from the control genes (Figure 5 — figure supplement 2A). No such differences were observed between DE and non-DE genes without overlapping DMRs (Figure 5E lower panel and Figure 5 — figure supplement 2A). These results suggest that the effect of H3K27me3 is likely specific to genes containing DMRs and the effect is stronger in the down-regulated DE genes.

We have added descriptions of these results in the manuscript as follows:

“We further examined whether the increased H3K27me3 could account for the reduced expression of some genes in the Dnmt3a cKO neurons. We found that down-regulated genes that overlap DMRs and the non-DE genes selected with matched expressions of up-regulated DE genes that overlap DMRs showed a small but significant median increase of H3K27me3 in the cKO (Wilcoxon rank-sum test p-value < 0.05; Figure 5E upper panel). When we considered a larger set of DE genes with a more relaxed threshold (FDR<0.2), we observed that down-regulated DE genes containing DMRs accumulate significantly more H3K27me3 than non-DE genes containing DMRs and up-regulated genes containing DMRs (Wilcoxon rank sum test p = 0.0029, Figure 5 — figure supplement 2A). By contrast, up-regulated genes had no significant accumulation of H3K27me3 and were not significantly different from the control genes (Figure 5 — figure supplement 2A). No such differences were observed between DE and non-DE genes without overlapping DMRs (Figure 5E lower panel and Figure 5 — figure supplement 2A). These results suggest that the effect of H3K27me3 is likely specific to genes containing DMRs and the effect is stronger in the downregulated DE genes, which may partially explain the fact that 24 genes were significantly down-regulated after the loss of repressive DNA methylation in the Dnmt3a cKO (Figure 2A-B).”

a. Figures 6E and S10A present the analysis that is most relevant to whether accumulation of K27me3 in the cKO impacts DEG and therefore may play a compensatory role. However, additional analysis could more decisively test the hypothesis. For example, do non-dysregulated genes containing DMRs show less K27me3 accumulation than down-regulated genes, but more K27me3 accumulation than upregulated genes? The authors could include a boxplot of a well-matched set of non-DE genes with overlapping DMRs in panel 6E to test this. A more powerful analysis, however, might be to examine a large set of genes from the genome that are selected to contain similar numbers of DMRs (in order to control for relative potential effects of demethylation) and ask if changes in K27me3 correspond with up-regulation (low K27me3 accumulation), no dysregulation (intermediate accumulation) or downregulation (high K27me3 accumulation).

We apologize for the confusion. We have added a panel using the same method used in Figure 5B (see “Enrichment test of DMRs and other genomic regions” in the Materials and methods section) to show the odds ratio of DMRs overlapping the H3K27me3 DM regions (now in Figure 4 — figure supplement 3D, right panel).

In addition, in the previous version of our manuscript, the enrichment of DMRs in DM regions was also assessed in the Venn diagram in Figure 4 — figure supplement 3D, and we have corrected the figure reference in the main text and the figure legend for Figure 4 — figure supplement 3. We tested this using the `fisher` sub-commands in the BEDTools Suite, where we checked if the amount of overlap between the 2 sets of regions is more than expected given their coverage and the size of the genome. This is in most cases analytically the same as shuffling the genome and checking the simulated (shuffled) versus the observed.

As the reviewer suggested, we have selected a more stringent set of regions as control. For each DM H3K27me3 peak, one non-DM H3K27me3 peak with a similar peak size and similar H3K27me3 signal in the control (wild-type) mice was selected as the matched control region using the greedy nearest neighbor matching in the R package `MatchIt`. Figure 4 — figure supplement 3C shows the distribution of H3K27me3 signal in the DM and non-DM peak sets:

We then counted the number of overlaps with DMRs in these matched non-DM peaks and compared those with the observed DM peaks (Figure 4 — figure supplement 3D, left panel). We found that the DM-peaks have significantly more overlaps with DMRs compared to the non-DM peaks (Fisher exact test p < 2.2e-16, odds ratio OR=2.70):

We have added these additional figures in Figure 4 — figure supplement 3 and described these results as follows in the main text:

“Likewise, P39 cKO DMRs were significantly enriched in peaks (22.4% overlapped H3K27me3 peaks in cKO, p<1e-100, Figure 5B) and DM regions of H3K27me3 (Figure 4 — figure supplement 3D). The DM regions of H3K27me3 had significantly more overlaps with DMRs compared to the non-DM control regions (Fisher exact test p < 2.2e-16, OR: 2.70, Figure 4 — figure supplement 3C-D).”

The methods used here are clarified in the Material and Methods section as follows:

“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).”

b. On p19 the authors state that "P39 cKO DMRs were significantly enriched in peaks and DM regions of H3K27me3 (28.1% overlapped H3K27me3 peaks in cKO, p<1e-100) (Figure 6B)." However, 6B only presents the overlap with peaks in each genotype as an odds ratio and doesn't make it clear how the enrichment of DMRs in DM regions was assessed. This appears to be assessed to some degree in S8F (note the S8 figure legend appears to be wrong for this panel), but no analysis is presented to control for the general tendency of DMRs and H3K27me3 peaks to overlap. It seems imperative, given that DMRs already overlap with H3K27me3 peaks, that a resampling approach is used comparing DM peaks to non-DM K27me3 peaks that are matched to have similar wild-type H3K27me3 signal and size as DM peaks. Likewise, what would the changes in H3K27me3 signal presented in 6C look like at non-DMR regions that contain similar baseline wild-type H3K27me3 signal? Without these analyses it is not clear that DMRs are particularly unique in showing increased K27me3 signal.

As for the reviewer’s suggestion of looking at non-DMR regions that contain similar baseline wild-type H3K27me3, it will be computationally intense to search for such matching regions as non-DMR regions are essentially the complement of DMR regions from the whole genome. And it will take time and luck to shuffle the DMR regions, compute the baseline H3K27me3 signal and get a match to the observed DMR regions. Instead, we addressed this comment with the following two approaches.

First, in Figure 5 — figure supplement 1A we analyze the ChIP-seq signal enrichment for two types of shuffle control regions. The first group of shuffle regions (dotted lines) is generated by shuffling the DMRs randomly across the same chromosome (excluding blacklist regions). These regions represent random background noise. Indeed, we find that the ChIP-seq signal is flat for these control regions in every sample and every histone mark. The second group of shuffle regions (dashed lines) is selected to meet the same criteria as the first group 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). These regions should at least have some P39 H3K27me3 signal, and as seen in the figure when averaged they formed a line that peaks at the center in the P39 H3K27me3 samples. Even though on average they have a higher enrichment of H3K27me3 compared to the observed DMRs, the differences between the cKO and Control are much smaller than those in the observed DMRs.

Secondly, we tiled the genome with 1kb bins and grouped them into bins with overlapping DMRs and bins without overlapping DMRs. For these two groups of bins, we plotted the mean changes of H3K27me3 in cKO vs. Control as a function of the baseline H3K27me3 signal in the control animals, stratified by the deciles of the baseline H3K27me3 signal (Figure 5 — figure supplement 1B). To avoid double-dipping (using the control signal in both the x-axis and y-axis), we used the baseline signal from the Control replicate 1 on the x-axis and the baseline signal from the Control replicate 2 on the y axis, and vice versa. In both cases, we observed bigger increases in H3K27me3 in bins with overlapping DMRs (Wilcoxon rank-sum test p < 0.001 in each bin), even after controlling for the baseline H3K27me3 level:

These new analyses indicate that DMRs are especially enriched in regions of increased H3K27me3. We have added these figures in Figure 5 — figure supplement 1A-B and described these results in the main text as follows:

“Moreover, H3K27me3 was more abundant at the center of DMRs in cKO compared to control neurons at P39 (∆ = 0.15 in the unit of fold enrichment vs. IgG, 5.04% increases, Figure 5C). There were much smaller changes of H3K27me3 at these DMRs in newborn (P0, ∆ = 0.0075, 0.27% changes) or fetal (E14, ∆ = -0.028, -1.02% changes) neurons (Figure 5C), and the increases were not seen in randomly shuffled regions (Figure 5 — figure supplement 1A). The signal of H3K4me3 and H3K27ac at the DMRs was also elevated in the cKO, to a lesser extent (H3K4me3: ∆ = 0.068, 2.92% changes; H3K27ac: ∆ = 0.13, 4.64% changes). … When accessing the H3K27me3 changes as a function of baseline H3K27me3 levels in the control across the genome tiled in 1kb bins, we observed bigger H3K27me3 increases in bins with overlapping DMRs compared to bins without overlapping DMRs (Wilcoxon rank-sum test p < 0.001 in each bin, Figure 5 — figure supplement 1B). These results indicate that DMRs are particularly unique in showing increased H3K37me3 signal.”

c. The increase in H3K27me3 in the cKO at the mab21l2 locus at P0 (figure 5B), before methylation has or has not been established, suggests that H3K27me3 may change before methylation differences, rather than in reaction to changes in methylation. Is this a common effect seen across other genes? How does this impact the idea that compensation occurs due to lack of methylation? Related to this issue, on p.19 no specific plot is referenced to support the statement that "There was no corresponding increase of H3K27me3 at these DMRs in newborn (P0) or fetal (E14) neurons." A quantitative analysis as shown in 6C but for the E14 and P0 time-points could assess this directly.

We apologize for the confusion. At the Mab21l2 locus (now in Figure 4B) there is no significant increase of H3K27me3 at P0 or E14 (i.e. there is no horizontal red bar under the P0 and E14 tracks; and our DiffBind analysis used a threshold of FDR<0.05). We have added a shaded box to highlight the region with significant increases of H3K27me3 at P39 to better illustrate our finding here (Also see our replies to reviewer #1, comment 3 for the quantification of the H3K27me3 levels at this locus).

We agree with the reviewer that quantitative analysis is needed to support our statement regarding the lack of increased H3K27me3 at DMRs at P0 and E14. Indeed, this analysis was included (now in Figure 5C) and we have updated the text to directly reference this figure panel as follows:

“Moreover, H3K27me3 was more abundant at the center of DMRs in cKO compared to control neurons at P39 (∆ = 0.15 in the unit of fold enrichment vs. IgG, 5.04% increases, Figure 5C). There were much smaller changes of H3K27me3 at these DMRs in newborn (P0, ∆ = 0.0075, 0.27% changes) or fetal (E14, ∆ = -0.028, -1.02% changes) neurons (Figure 5C), and the increases were not seen in randomly shuffled regions (Figure 5 — figure supplement 1A).”

Finally, we have also shown the H3K27me3 changes did not happen in earlier timepoint as we found almost no significant differentially modified H3K27me3 regions in E14 or P0 (Figure 4 — figure supplement 1D). Also, if we focus on the 4,040 significant up-regulated H3K27me3 regions found in P39, the increased signal changes in cKO vs. Control in P39 are not correlated with those in E14 or P0 (Figure 4 — figure supplement 2B and D).

d. Why were developmental DMs called between E14 and P39 while DMRs were called between P0 and P39? Developmental comparisons would be more justified if both used P0 and P39.

Our choice to analyze the DMs between E14 and P39 was motivated by several characteristics of our dataset. Unfortunately, we do not have methylation data for E14 excitatory neurons, therefore developmental DMRs could only be called between P0 and P39. Moreover, we found much more significant developmental DM H3K27me3 regions in P39 vs. E14 compared to P39 vs. P0 (Figure 6 — figure supplement 1A). We believe that one of the reasons for this is that the E14 ChIPseq data are of better quality than the perinatal data from P0 (Figure 6 — figure supplement 2). For this reason, in our original submission, we chose to use E14 and P39 to call developmental H3K27me DMs. We now provide a new figure for the reviewers (Figure R2) to show corresponding results for Figure 6B-E and Figure 6 — figure supplement 1C-G using the set of developmental H3K27me3 DMs between P39 and P0. In these figures, we demonstrate that the conclusions drawn in the previous version are verified in the P39 vs. P0 comparison.

e. In the absence of additional data supporting H3K27me3 accumulation as a functional compensation mechanism for loss of DNA methylation, the title of Figure 6 should not so strongly claim that H3K27me3 changes compensate for mC loss.

We have revised the title of this figure (now Figure 5): Increased H3K27me3 correlates with the loss of postnatal DNA methylation.

3. While initially assessed by the authors, mCH changes are not integrated into the analysis of H3K27me3 changes. Previous studies by this group and others have suggested the importance of mCH in gene regulation within specific neuron populations and, given the interest to the field, more thorough analysis is warranted regarding this neuronal methylation and H3K27me3.

a. Given the association, even if limited, between differential mCH loss and gene dysregulation in the cKO, it is relevant to examine how changes in K27me3 relate to, or compensate for, loss of mCH. For example, do genes that are not dysregulated in the cKO, and normally have high mCH (which is lost in the cKO) show higher increases in K27me3 compared to non-dysregulated genes that have low mCH (and therefore lose less mCH)? Reciprocally, one could ask if a population of genes selected to have similar, high mCH loss in the cKO, show differential gene dysregulation as a function of differential K27me3 accumulation.

We appreciate the reviewer’s insightful suggestions for analysis of the relationship between the loss of mCH and the changes in the H3K27me3 signal. We agree that since the loss of mCH and the loss of mCG are highly correlated (Figure 3 — figure supplement 3D), we would expect to see to a certain extent a similar compensation of H3K27me3 in genes/regions with loss of mCH. We have added a figure to Figure 5 — figure supplement 2B to show the changes in the H3K27me3 signal as a function of the loss of gene body mCH in non-DE genes (stratified by the deciles of δ mCH):

We indeed observed a larger increase of H3K27me3 signal in the gene body of genes that lost most mCH in the cKO (genes that have high mCH in the control; the leftmost decile in the plot), compared to genes that lost less mCH (genes that normally have low mCH).

In addition, we grouped the genes by the quantiles of loss of mCH in cKO and plotted the gene expression fold-changes as a function of the changes of H3K27me3 in gene body for each of these gene sets (Figure 5 — figure supplement 2C):

This shows that genes that lost the most mCH in the cKO (genes that have high mCH in the control; red line) tend to have increased expression in the cKO, compared to genes that lost the least mCH (genes that have low mCH to begin with; blue line), which is consistent with our results in Figure 3E. Moreover, we can see a negative trend between the changes of gene body H3K27me3 signal and gene expression fold-change in the gene set with the highest baseline mCH, in which genes that showed the biggest increases of H3K27me3 were inclined to decrease in expression whereas genes that had small increases or decreases of H3K27me3 signal were generally up-regulated in gene expression. Such correlation was not observed in genes with the lowest baseline mCH.

We have added these plots in Figure 5 — figure supplement 2B-C and described the results in the main text as follows:

“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.”

b. What are the changes in mCH at DM H3K27me3 sites? Since mCG and mCH loss are correlated, one might expect that these are sites that normally also contain high mCH. It might not be feasible to call mCH DMRs given the technical limitations, but quantification of this signal can be done at DMs and DMRs, and the potential for it to contribute to these effects could be discussed.

As suggested by the reviewer, we have estimated the changes of mCG and mCH in the DM H3K27me3 regions and also compared them with the changes in matched non-DM regions (see our response to comment 2b for the methods we used to select these matched non-DM regions). We found that the DM regions showed significantly larger decreases in both mCG and mCH when compared to non-DM regions (Wilcoxon rank-sum test p < 0.0001).

We have added these plots in Figure 4 — figure supplement 3E and described the results in the main text as follows:

The DM regions also have larger decreases in both mCG and mCH when compared to non-DM regions (Wilcoxon rank-sum test p < 0.0001, Figure 4 — figure supplement 3E).

4. The authors present evidence of PRC2 binding at DEGs (Figure 5A) as a motivation to assess H3K27me3 in their mice. However, they then observe changes in H3K27me3 that they contend do not occur under normal conditions and propose that these changes are an underlying mechanism that blocks gene dysregulation in these mice. How are we to reconcile this original rationale of a normal function for PRC2 at dysregulated genes with the model of de novo compensatory H3K27me3 accumulation? Some explanation by the authors of how these separate notions may relate to one another would be helpful.

The BART results for up-regulated and down-regulated DE genes (now in Figure 4A) were bioinformatics predictions that drove us to look into the Polycomb repression via H3K27me3 profiling. BART predictions are based on ChIP-seq data in various cell lines and tissues from the ENCODE project. This means that the potential PRC2 binding at our DEGs may normally happen in systems other than the brain or pyramidal neurons, or at other time points during development. Our results showed that these DEGs were normally regulated by DNA methylation in the pyramidal neurons but they can still be regulated by PRC2 in the current system if the DNA methylation was abolished.

We have added these explanations in the Discussion section of the manuscript as follows:

“Overall, our results suggest that when DNA methylation is disrupted, H3K27me3 might partially compensate for the loss of mCG and/or mCH and act as an alternative mode of epigenetic repression. Nevertheless, we did not find differential expression in any of the four core components of PRC2 (Ezh2, Suz12, Eed and Rbbp4) in adult Dnmt3a cKO animals. It is possible that the increased H3K27me3 was mediated by transient expression of PRC2 components during development in the cKO. Furthermore, the predictions from BART (Figure 4A) were derived from various cell lines and tissues from the ENCODE project (Davis et al., 2018; ENCODE Project Consortium, 2012), suggesting that the potential PRC2 binding at our DEGs may normally happen in systems other than the brain or pyramidal neurons, or at other time points during development. Additional experiments which directly manipulate components of the PRC2 system are required to further test the potential compensation mechanism.”

5. Some effects on H3K27ac are observed at DMRs, but the authors ignore the prospect that these effects could drive changes in gene expression. Previous groups have shown significant changes in H3K27ac in DNMT3A models, and some of the data the authors present suggest that this could also explain some of the observations in their model.

a. While DM K27ac peaks are not identified in the cKO, in figure 6B there is a similar enrichment of DMRs to associate with K27ac peaks in the cKO vs control as there is for K27me3. In addition, there is a similar increase in odds ratio of enrichment of K27ac peaks to the of K27me3 at DMRs for the cKO in figure 6B. The analysis in S8F-G that is referenced when referring to K27ac is either missing (S8G) or is not showing data for K27ac (S8F). The authors should at least comment in the Results section that these signals exist and may be significant for gene regulation.

In the new set of DMRs from the new batch of DNA methylation data (see responses to the editor on why we introduced a new batch of sequencing data), we didn’t see the enrichment of overlapping DMRs in H3K27ac (now in Figure 5B). But indeed, we still see the enrichment of the P39 H3K27ac signal in the center of DMRs (Figure 5C). We decided to focus on the changes of H3K27me3 since it’s more evident, but we agree with the reviewer that H3K27ac may also have an impact on regulating the expression of the DEGs. We apologize for the confusion in referring to Figure S8F-G, and we have corrected the text and added some comments regarding H3K27ac as follows:

“The signal of H3K4me3 and H3K27ac at the DMRs was also elevated in the cKO, but to a lesser extent (H3K4me3: ∆ = 0.068, 2.92% changes; H3K27ac: ∆ = 0.13, 4.64% changes).”

b. Are two biological replicates sensitive enough to call differences in H3K27ac? In the context of a small number of replicates and the subtle magnitude of changes in gene expression that are observed, negative results for differential peaks should be interpreted with caution. In contrast, the larger peaks observed in K27me3 may have led to increased ability to call differential peaks.

We agree with the reviewer that our power for discovering changes in H3K27ac is somewhat limited by the sample size of two replicates per condition. We have added some comments to remind the readers of this potential issue in the Discussion section as follows:

“Indeed, we found that the Dnmt3a cKO DMRs overlapped with regions that gain methylation during normal postnatal development (Figure 3G-H). These DMRs had increased H3K27ac in the Dnmt3a cKO (Figure 5C), consistent with observations in adult animals lacking Dnmt3a specifically in GABAergic Sst- or Vip-expressing interneurons (Stroud et al., 2020). In those experiments, embryonic gene-regulatory elements had lower cytosine methylation and increased H3K27ac and H3K4me1 (Stroud et al., 2020). This suggests an essential role for Dnmt3a and DNA methylation in shaping the transcriptome during development in part via the inactivation of embryonic enhancers, with potentially long-lasting effects on the gene expression pattern of mature neurons. Although we did not detect individual peaks with a statistically significant difference in H3K27ac between cKO and control (Figure 4 — figure supplement 1C), this could be due to a small number of replicates and the subtle magnitude of changes in gene expression relative to the large number of tested peaks. Future experiments using more replicates or orthogonal techniques like CUT&RUN and CUT&Tag could help to further investigate the relationship between the Dnmt3a-dependent DNA methylation and the activation of enhancers.”

6. Though the characterization of DNMT3A cKO in excitatory neurons is novel, if the authors are claiming that these neurons specifically are driving behavioral phenotypes then further contextualization and characterization of the extent of conditional deletion is needed.

a. It would be useful to justify the rationale for investigating excitatory neurons, and how using this class of neurons can uncover novel findings. Existing studies of DNMT3A cKO have shown behavioral and epigenomic effects in pan-neuronal and inhibitory neuron conditional knockouts. Comparisons between these models and discussion of differences could offer strong insights into what unique circuits may be most disrupted by DNMT3A cKO in excitatory neurons.

Our previous studies (e.g. Lister et al., 2013; Mo et al. 2015; Luo et al. 2017) have shown that CG and non-CG DNA methylation across all neuron types has a highly cell type specific distribution that is globally remapped during postnatal brain development. To understand the functional significance of this epigenetic regulation, we chose to focus on the major neuronal populations with established roles in regulating cognitive behavior. Our previous study (Lavery et al., 2020) analyzed the loss of Dnmt3a in GABAergic neurons. Here, we chose to focus on the other major neuronal population, the excitatory (glutamatergic) neurons in the frontal cortex.

Note that both pan-neuronal and pan-GABAergic cKOs (e.g. (Feng et al., 2010; Morris et al., 2014; Nguyen et al., 2007)) affect many more regions of the brain than what is affected in our study, thus leading to confounding compensatory/deleterious phenotypes. Previous studies using those cKOs had profound phenotypes and died before reaching early adulthood. Because we wanted to analyze the effect of blunting the increase in Dnmt3a during the perinatal period in a cell-specific manner, without killing the mice, we believe our approach has allowed for a much more detailed study of the regulation of neuronal DNA methylation. Although the phenotypes observed in the present cKO are mild, both at the functional and transcriptional level, we still found significant changes that are more reminiscent of those found in neurodevelopmental disorders. Dnmt3a de novo mutations have been linked to autism spectrum and other developmental disorders. Our present results may point to the period of gain of methylation in postmitotic excitatory neurons as a critical period in the development of these disorders.

b. If the authors are claiming that behaviors and phenotypes are the result of specifically knocking out DNMT3A in all excitatory neurons, further evidence should be provided that the Cre-line is being expressed sufficiently and specifically in excitatory neurons, and that their genomic data is faithfully capturing the majority of excitatory cells while excluding non-excitatory cells.

We thank the reviewer for raising this question. The Neurod6-Cre mouse line (also referred to as NEX-Cre) used in our study is a knock-in mouse line that expresses Cre recombinase under the control of Neurod6 regulatory sequences (Goebbels et al., 2006). Neurod6 is a transcription factor exclusively expressed in the central nervous system – most prominently in excitatory neurons of the neocortex and hippocampus (e.g. see Yao et al., Cell 2021), and its expression parallels excitatory neuronal differentiation and synaptogenesis during brain development (see Figure 1 — figure supplement 1B).

Consistent with the expression of NeuroD6, Goebbels et al. (2006) observed the most prominent Cre activity in the neocortex and hippocampus, starting from around embryonic day 11.5. Moreover, Cre-mediated recombination marked pyramidal neurons and dentate gyrus mossy and granule cells and was absent from proliferating neural precursors of the ventricular zone, interneurons, oligodendrocytes, and astrocytes (Goebbels et al., 2006). We confirmed that Neurod6-dependent Cre recombination occurred only in excitatory neurons using the INTACT mouse (Mo et al., 2005), in which Cre-mediated excision of the transcription stop signals activates the expression of the nuclear membrane tag (Sun1 -sfGFP-myc) in the cell type of interest. Since recombination only occurs in excitatory neurons, we expect excision of exon 19 from Dnmt3a will only occur in those neurons expressing Cre recombinase. Figure 1 — figure supplement 2A-C depicts the staining of Sun1-GFP across the brain in non-inhibitory Neurod6+ neurons, confirming previous results (Goebbles et al., 2006). Note that the sfGFP-tag decorates the nucleus of cells where recombination occurs.

We also added the following Figure 1 — figure supplement 2D to show the expression of pan-neuronal, excitatory neuron and inhibitory neuron marker genes in our RNA-seq data, in which we observed high expression of excitatory neuron markers and very low expression of inhibitory neuron markers, further demonstrating that the neurons carrying the GFP label belong to the excitatory neuron population.

7. The manuscript would be further improved by clarification of some methods and statistics, and modifications in data presentation:

a. The authors should provide information on the statistics performed on the behavior data. Are there significant differences between sexes, or significant sex by genotype interactions for these behavioral tasks? A two-way ANOVA using genotype and sex would be useful in understanding the impact of these variables. There appears to be many effects that are sex-specific, yet this appears to be overlooked in the text. Furthermore, if there are robust sex-specific effects then the subsequent molecular and genomic analyses should be done using both sexes, and comparisons could be done to explain the molecular underpinnings of these differences.

We thank the reviewer for raising this point, for which we would like to provide clarification. The experiments in this section (Figures 1, Figure 1 — figure supplement 4 and 5) were first and foremost designed to investigate whether the loss of Dnmt3a-dependent DNA methylation during the maturation of pyramidal cells might lead to endophenotypic markers of neurodevelopmental illnesses. Male mice were used in generating all genomic datasets. We performed behavioral experiments in cohorts of male mice, which led to the discovery of impairments in social behavior, PPI, and working memory. We subsequently investigated if those same behaviors were affected in a cohort of female mice, finding that the loss of Dnmt3a caused similar impairments in working memory, but not in PPI or social behavior. Because the male and female cohorts were tested in completely independent experiments, it would not be appropriate to perform a two-way ANOVA using sex as a variable, as this could result in an inaccurate depiction of results. Yet, the wider range of behavioral impairments identified in the male cohort is noteworthy, as it correlates well with human neurodevelopmental disorders, such as schizophrenia and autism, in which males are disproportionately more affected than females. We have edited the Discussion to clarify that even though it was outside the scope of the present work, our data provide support for designing a comprehensive follow up study aiming to dissect sex-specific effects of Dnmt3a and their underlying causes.

b. The authors might also consider showing preference index when analyzing social approach data when comparing genotypes.

We thank the reviewer for this suggestion, which we have carefully considered. We decided to keep our results as raw "time spent on side" to take into account the differing time spent in the empty space, which showed high significance and would be lost if the results were reported as preference index. We edited the results text in order to mention this finding clearly:

“Moreover, when tested in a three-chamber box in which one of the sides contained a novel mouse and the opposite a novel object, male Dnmt3a cKO animals spent less time in the former, opting, instead, to remain significantly longer in the center (empty compartment), which is suggestive of reduced exploration and social interest (Figure 1D, left panel; p = 0.01048).”

c. What is the background set of genes for 2B-C? Synaptic genes are enriched in both up- and down-regulated gene lists. How is gene expression normalized? Is there a discovery bias for more highly expressed genes?

In the previous version of our manuscript, we used expressed genes in P39 as the background set in the GO analysis for the DE genes. Expressed genes were defined as genes with CPM (counts per million) greater than 2 in at least two samples. In the new batch of RNA-seq data (see responses to the editor on why we introduced a new batch of sequencing data), we identified fewer significant DE genes and hence found no significant enrichments in any GO terms. Therefore, we have removed all the figures regarding the GO analysis for the DE genes.

d. What value was used to identify equally expressed unchanged genes in figure 4B, TPM? Identifying equally expressed genes can be complicated for nuclear RNA where unprocessed RNAs contribute differentially to total counts based on the length of the gene. In addition, previous analysis of mCH enrichment at dysregulated genes in inhibitory neuron-specific KO of DNMT3A also found DEGs to be differentially methylated based on direction of change but showed that all other genes were more lowly methylated than these two groups (Lavery et al. 2020). Here the authors see the same relationship between DEG groups but find other genes to be higher in mCH rather than lower. Do the authors have any perspective on why these differences are observed?

We used the R package `MatchIt` to select non-DE genes with matched expression levels in the unit of TPM. We have added Figure 3 — figure supplement 4B to show these non-DE genes have similar expression levels as the DE genes in the P39 Control samples, and their expressions were not changed in the P39 cKO samples.

With the new batch of RNA-seq and MethylC-seq data with bona fide biological replicates (see responses to the editor on why we introduced a new batch of sequencing data), we updated the mCH profile in the DEGs and corresponding non-DE control genes (now in Figure 3D):

We now see that the up-regulated genes have higher mCH than the downregulated genes, and non-DE genes have lower mCH than the up-regulated genes but also indistinguishable levels from the down-regulated genes. This in general agrees with the findings in Lavery et al. We presumed the discrepancy between the findings in Lavery et al. and our previous version of the manuscript might be due to the differences in the power of detecting DE genes, the choice of thresholding significant DEGs, the method of selecting non-DE control genes, and probably some unknown distinctions of biological mechanism in regulating gene expression after Dnmt3a cKO in inhibitory neurons vs. in excitatory neurons.

e. On p22 the authors state "We found that developmental peaks (DevLoss and DevGain) were relatively unaffected by the cKO(H3K27me3 = 0.11, -0.20 respectively, in units of log10(CPM+1)), whereas Group cKO had ten-fold larger mean effect (H3K27me3 = 1.12 respectively)." Given the y scales in figures 7E and S10H it is unclear how they derived a ∆ value of 1.12 on a log 10 scale for Group cKO (e.g. in S10H the median values look much less 1 unit different on a log10 scale). Please clarify the actual magnitude of effect sizes.

We are sorry for the confusion. In the main figure (now in Figure 6C-E) we normalized the H3K27me3 signal to the means of the two P39 Control samples and plotted the fold changes. In the said supplemental figures (now in Figure 6 — figure supplement 1E-G) we plotted the H3K27me3 signal in units of log10(CPM+1) for each sample. In the previous version of our manuscript, the δ values were achieved by running a paired t-test on the raw CPM scale. To avoid confusion, we have rerun the paired t-test on the log10(CPM+1) scale and updated the δ values in Figure 6 — figure supplement 1E-G. We have also edited the relevant statements in the main text as follows:

“To further stratify the joint distribution of developmental and Dnmt3a cKOdependent changes in both H3K27me3 and DNA methylation, we assigned peaks to three groups (Figure 6C-E and Figure 6 — figure supplement 1EG). “Group DevLoss” and “Group DevGain” peaks lose or gain H3K27me3 during development, respectively (Figure 6C-D and Figure 6 — figure supplement 1E-F). “Group cKO” peaks have higher H3K27me3 in the Dnmt3a cKO compared to control at P39 (Figure 6E and Figure 6 — figure supplement 1G). We found that developmental peaks (DevLoss and DevGain) were relatively unaffected by the cKO (H3K27me3 = 0.02, -0.03 respectively, in units of log10(CPM+1)), whereas Group cKO had 4.5-fold larger mean effect (H3K27me3 = 0.09 respectively). Group cKO peaks also experienced greater loss of mCG (mCG = -23.5%) than Group DevLoss (13.0%) or DevGain (18.6%) (Figure 6C-E and Figure 6 — figure supplement 1E-G, middle and right panels). These results suggest that regions prone to alteration of H3K27me3 by Dnmt3a cKO are distinct from the regions affected by developmentally dynamic H3K27me3.”

f. The authors describe the association between gene body mCH and gene expression changes as only accounting for 1% of the variance. It would be useful to have a point of reference to understand if this is a low value or a high value in comparison to other potential predictors of gene expression changes. For example, how much does the proposed overcompensation by H3K27me3 accumulation contribute to the genome-wide variance in gene expression? In addition, is it the case that variance in the measured fold-changes will be mitigated by addition of more replicates to more accurately estimate true fold-changes? If this is the case, it should be noted.

Thank you for the valuable suggestion to provide a point of reference for the magnitude of the correlation between mCH and gene expression. The variance explained (R-squared) between the gene body mCH and RNA log fold-change was 0.46%. To give a quantitative context for this, we estimated the total explainable variance as the R-squared between RNA log fold-change from independent biological replicates. That is, we selected one of the cKO and one control sample to calculate the first replicate logFC value; we then used the second cKO and second control sample to get a second, independent replicate logFC estimate. We found that the correlation between these was 1.30%. This shows that, although the mCH R-squared is small, it is nevertheless a significant fraction of the total explainable variance given the noise in our RNA-seq data. As the reviewer suggests, it is likely that adding more data (more biological replicates) would reduce the noise in both the RNA logFC and the mCH measurements, which might increase the explained variance. We have revised the Results as follows:

“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.”

Furthermore, we used a linear model to test how much the δ H3K27me3 signal may contribute to the variance of RNA logFC between cKO and control in 14,755 expressed genes. We found that the R-squared value is 0.0599%, which is smaller than the R-squared for mCH (0.456%, see Figure 3 — figure supplement 4D).

g. Using different color schemes would be helpful in 8A and 8B to avoid confusion.

We have modified the said figure (now in Figure 7A-B) as suggested to avoid confusion.

h. How was clustering analysis done in 8D? What do the values for K27me3 Diff. Mod and DEG enrichment indicate (i.e. what does the scale from 0 to 0.02 represent)?

The clustering of DMVs (now in Figure 7C) was originally performed using k-means. In the resubmission, we have updated this analysis to make the results more evident, and as a result we modified the clustering procedure. The new procedure sorts 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).

We now show directly the mean mRNA logFC for all genes contained within each DMV.

We have added the description of this analysis in the method section as follows:

“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.”

Reviewer #3 (Recommendations for the authors):

– The authors delete exon 19 of Dnmt3a. I think they should provide a more thorough description of what this deletion is expected to cause: loss of function due to premature stop codon, loss of enzymatic activity, etc. Moreover, it would be desirable to show by immunostaining whether DNMT3A protein is specifically lost in excitatory neurons.

We thank the reviewer for raising this question and allowing us to clarify an important point. Deletion of exon 19 in the mouse Dnmt3a gene is expected to produce a deletion of 50 amino acids (aa 623-673) in the methyltransferase domain (encompassing the S-adenosyl-L-methionine binding site), thus leading to disruption of the catalytic activity of the enzyme (Lyko 2018). In addition, this approach resulted in a marked downregulation of Dnmt3a protein levels, as quantified by western blot experiments (Figure 1 — figure supplement 3B). Such a downregulation in Dnmt3a expression was matched by an almost complete loss of mCH and loss of mCG remethylation in excitatory neurons confirming the lack of activity of Dnmt3a in excitatory neurons from frontal cortex (Figure 3C). Please also see our response to Reviewer #2 comment 6b where we performed triple labeling of brain slices, showing that the Cre-dependent recombination, and thus the deletion of Dnmt3a, occurs in excitatory, but not inhibitory, neurons.

– Regarding the DMRs identified in DNMT3A KO neurons, the authors should use the generated profiles for H3K27ac in P0 and P39 WT neurons to determine which DMRs overlap enhancers that are active in any of these two neuronal stages. This is more relevant than overlaps with enhancers that could be active in any given cell type. Most importantly, the authors should evaluate the impact the loss of DNA methylation can have on enhancer activity by measuring H3K27ac levels and eRNAs in DNMT3A KO neurons. Measuring eRNA levels can be particularly relevant, as recent evidences suggest that is a better marker of enhancer activity than H3K27ac.

We agree that it is important to analyze the impact of the loss of DNA methylation on enhancer activity. For P39 samples, we analyzed the overlap of H3K27ac peaks in pyramidal neurons with DMRs (Figure 5A-B and Figure 5 — figure supplement 1A). Moreover, in Figure 5C we plotted the mean level of H3K27ac in the region around DMRs. These analyses show that enhancer activity is indeed increased in regions that lose DNA methylation in the cKO. Please also see our responses to Reviewer #1 comment 3 and Reviewer #2 comment 5 regarding the analyses on H3K27ac.

We would like to clarify that we did not perform H3K27ac ChIP-seq for P0, so we were unable to specifically examine the association of enhancer activity at the perinatal timepoint.

We agree that it would be interesting to further examine the enhancers using eRNAs, however, we found that the depth of coverage of our RNA-seq data was not sufficient to address this.

– The authors should aim at performing experiments that address the proposed compensatory role for PRC2 and H3K27me3. I understand that generating a KO for PRC2 could be rather complicated, but perhaps the authors could consider using new drugs that act as very specific inhibitors against the enzymatic activity of PRC2.

We appreciate the reviewer’s suggestion, and we agree that it will need further experiments to directly test the hypothesis of the compensatory role for PRC2 and H3K27me3. But such experiments are out of the scope of this manuscript. We have modified the text to weaken our claims regarding the compensation effect of H3K27me3.

https://doi.org/10.7554/eLife.66909.sa2

Article and author information

Author details

  1. Junhao Li

    Department of Cognitive Science, University of California, San Diego, La Jolla, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review and editing
    Contributed equally with
    Antonio Pinto-Duarte
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6784-3780
  2. Antonio Pinto-Duarte

    Computational Neurobiology Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    Present address
    Ionis Pharmaceuticals, Carlsbad, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Validation, Visualization, Writing – original draft, Writing – review and editing
    Contributed equally with
    Junhao Li
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2215-7653
  3. Mark Zander

    Genomic Analysis Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    Present address
    Waksman Institute of Microbiology, Rutgers, The State University of New Jersey, Piscataway, United States
    Contribution
    Data curation, Investigation, Methodology, Validation, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8643-1407
  4. Michael S Cuoco

    Bioinformatics and Systems Biology Graduate Program, University of California, San Diego, La Jolla, United States
    Contribution
    Formal analysis
    Competing interests
    No competing interests declared
  5. Chi-Yu Lai

    Computational Neurobiology Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    Contribution
    Data curation, Investigation, Methodology, Validation
    Competing interests
    No competing interests declared
  6. Julia Osteen

    Computational Neurobiology Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    Contribution
    Data curation, Investigation, Methodology, Validation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7058-3297
  7. Linjing Fang

    Waitt Advanced Biophotonics Core, Salk Institute for Biological Studies, La Jolla, United States
    Present address
    Division of Biology and Biological Engineering, California Institute of Technology, Pasadena, United States
    Contribution
    Data curation, Investigation, Methodology, Resources, Validation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2232-2601
  8. Chongyuan Luo

    1. Genomic Analysis Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    2. Howard Hughes Medical Institute, Salk Institute for Biological Studies, La Jolla, United States
    Present address
    Department of Human Genetics, University of California, Los Angeles, Los Angeles, United States
    Contribution
    Data curation, Formal analysis, Investigation, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8541-0695
  9. Jacinta D Lucero

    Computational Neurobiology Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    Contribution
    Investigation, Methodology, Resources
    Competing interests
    No competing interests declared
  10. Rosa Gomez-Castanon

    Genomic Analysis Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
  11. Joseph R Nery

    1. Genomic Analysis Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    2. Howard Hughes Medical Institute, Salk Institute for Biological Studies, La Jolla, United States
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
  12. Isai Silva-Garcia

    Computational Neurobiology Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    Contribution
    Investigation, Methodology, Validation
    Competing interests
    No competing interests declared
  13. Yan Pang

    Computational Neurobiology Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    Contribution
    Investigation, Methodology, Resources
    Competing interests
    No competing interests declared
  14. Terrence J Sejnowski

    Computational Neurobiology Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    Contribution
    Resources, Supervision
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0622-7391
  15. Susan B Powell

    Department of Psychiatry, University of California, San Diego, La Jolla, United States
    Contribution
    Funding acquisition, Investigation, Methodology, Resources, Supervision, Writing – review and editing
    Competing interests
    No competing interests declared
  16. Joseph R Ecker

    Genomic Analysis Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    Contribution
    Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Writing – review and editing
    For correspondence
    ecker@salk.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5799-5895
  17. Eran A Mukamel

    Department of Cognitive Science, University of California, San Diego, La Jolla, United States
    Contribution
    Conceptualization, Formal analysis, Funding acquisition, Project administration, Resources, Supervision, Visualization, Writing – original draft, Writing – review and editing
    For correspondence
    emukamel@ucsd.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3203-9535
  18. M Margarita Behrens

    1. Computational Neurobiology Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    2. Department of Psychiatry, University of California, San Diego, La Jolla, United States
    Contribution
    Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Writing – original draft, Writing – review and editing
    For correspondence
    mbehrens@salk.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7168-8186

Funding

National Institute of Mental Health (R01MH112763)

  • Joseph R Ecker
  • Eran A Mukamel
  • M Margarita Behrens

Kavli Foundation

  • Antonio Pinto-Duarte
  • Susan B Powell
  • M Margarita Behrens

Howard Hughes Medical Institute

  • Joseph R Ecker

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was supported by R01MH112763 to MMB and JRE, and a Kavli Foundation award to MMB, APD, and SBP. JRE is an Investigator of the Howard Hughes Medical Institute. We acknowledge stimulating discussions with Drs. Huda Zoghbi and Laura Lavery. We thank the members of the Salk Biophotonics Core, Dr Uri Manor, Sammy Weiser Novak, and Dr Tong Zhang for their insightful suggestions. We also thank Joseph Chambers and Caitlin Chambers for technical assistance in animal handling, Colleen Heller for technical assistance in behavioral experiments, and Faith Zhang for her initial involvement in the morphometric analysis of dendritic spines. The Waitt Advanced Biophotonics Core Facility at the Salk Institute receives funding from NIH-NCI CCSG: P30 014195 and the Waitt Foundation. The Flow Cytometry Core Facility of the Salk Institute receives funding from NIH-NCI CCSG: P30 014195. The authors have no conflict of interest in relation to the work described here.

Ethics

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).

Senior Editor

  1. Lu Chen, Stanford University, United States

Reviewing Editor

  1. Anne E West, Duke University, United States

Reviewer

  1. Alvaro Rada-Iglesias, University of Cologne, Germany

Publication history

  1. Preprint posted: December 20, 2020 (view preprint)
  2. Received: January 28, 2021
  3. Accepted: May 22, 2022
  4. Accepted Manuscript published: May 23, 2022 (version 1)
  5. Version of Record published: June 6, 2022 (version 2)

Copyright

© 2022, Li, Pinto-Duarte et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,070
    Page views
  • 335
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Junhao Li
  2. Antonio Pinto-Duarte
  3. Mark Zander
  4. Michael S Cuoco
  5. Chi-Yu Lai
  6. Julia Osteen
  7. Linjing Fang
  8. Chongyuan Luo
  9. Jacinta D Lucero
  10. Rosa Gomez-Castanon
  11. Joseph R Nery
  12. Isai Silva-Garcia
  13. Yan Pang
  14. Terrence J Sejnowski
  15. Susan B Powell
  16. Joseph R Ecker
  17. Eran A Mukamel
  18. M Margarita Behrens
(2022)
Dnmt3a knockout in excitatory neurons impairs postnatal synapse maturation and increases the repressive histone modification H3K27me3
eLife 11:e66909.
https://doi.org/10.7554/eLife.66909

Further reading

    1. Cell Biology
    2. Genetics and Genomics
    Heyun Guo et al.
    Research Article

    In the first meiotic cell division, proper segregation of chromosomes in most organisms depends on chiasmata, exchanges of continuity between homologous chromosomes that originate from the repair of programmed double-strand breaks (DSBs) catalyzed by the Spo11 endonuclease. Since DSBs can lead to irreparable damage in germ cells, while chromosomes lacking DSBs also lack chiasmata, the number of DSBs must be carefully regulated to be neither too high nor too low. Here, we show that in Caenorhabditis elegans, meiotic DSB levels are controlled by the phosphoregulation of DSB-1, a homolog of the yeast Spo11 cofactor Rec114, by the opposing activities of PP4PPH-4.1 phosphatase and ATRATL-1 kinase. Increased DSB-1 phosphorylation in pph-4.1 mutants correlates with reduction in DSB formation, while prevention of DSB-1 phosphorylation drastically increases the number of meiotic DSBs both in pph-4.1 mutants as well as in the wild type background. C. elegans and its close relatives also possess a diverged paralog of DSB-1, called DSB-2, and loss of dsb-2 is known to reduce DSB formation in oocytes with increasing age. We show that the proportion of the phosphorylated, and thus inactivated, form of DSB-1 increases with age and upon loss of DSB-2, while non-phosphorylatable DSB-1 rescues the age-dependent decrease in DSBs in dsb-2 mutants. These results suggest that DSB-2 evolved in part to compensate for the inactivation of DSB-1 through phosphorylation, to maintain levels of DSBs in older animals. Our work shows that PP4PPH-4.1, ATRATL-1, and DSB-2 act in concert with DSB-1 to promote optimal DSB levels throughout the reproductive lifespan.

    1. Genetics and Genomics
    2. Microbiology and Infectious Disease
    Liselot Dewachter et al.
    Research Article

    Antibiotic resistance in the important opportunistic human pathogen Streptococcus pneumoniae is on the rise. This is particularly problematic in the case of the β-lactam antibiotic amoxicillin, which is the first-line therapy. It is therefore crucial to uncover targets that would kill or resensitize amoxicillin-resistant pneumococci. To do so, we developed a genome-wide, single-cell based, gene silencing screen using CRISPR interference called sCRilecs-seq (subsets of CRISPR interference libraries extracted by fluorescence activated cell sorting coupled to next generation sequencing). Since amoxicillin affects growth and division, sCRilecs-seq was used to identify targets that are responsible for maintaining proper cell size. Our screen revealed that downregulation of the mevalonate pathway leads to extensive cell elongation. Further investigation into this phenotype indicates that it is caused by a reduced availability of cell wall precursors at the site of cell wall synthesis due to a limitation in the production of undecaprenyl phosphate (Und-P), the lipid carrier that is responsible for transporting these precursors across the cell membrane. The data suggest that, whereas peptidoglycan synthesis continues even with reduced Und-P levels, cell constriction is specifically halted. We successfully exploited this knowledge to create a combination treatment strategy where the FDA-approved drug clomiphene, an inhibitor of Und-P synthesis, is paired up with amoxicillin. Our results show that clomiphene potentiates the antimicrobial activity of amoxicillin and that combination therapy resensitizes amoxicillin-resistant S. pneumoniae. These findings could provide a starting point to develop a solution for the increasing amount of hard-to-treat amoxicillin-resistant pneumococcal infections.