Loss of Ezh2 in the medial ganglionic eminence alters interneuron fate, cell morphology and gene expression profiles

1 Introduction

Inhibitory GABAergic interneurons are a heterogeneous cell population with dozens of subtypes displaying distinct morphologies, connectivity, electrophysiology properties, neurochemical markers and gene expression profiles (Bandler et al., 2017; Hu et al., 2017b; Lim et al., 2018; Williams and Riedemann, 2021). Perturbation of interneuron development and inhibition is associated with a range of disorders including epilepsy, schizophrenia and autism (Valiente and Marin, 2010; Bozzi et al., 2012; Inan et al., 2013; Contractor et al., 2021), and many disease-associated genes are enriched in prenatal immature interneurons and affect their development (Schork et al., 2019; Trevino et al., 2020; Paulsen et al., 2022). Forebrain interneurons originate from two transient structures in the embryonic ventral forebrain, the medial and caudal ganglionic eminences (MGE and CGE, respectively), and mature over the course of embryonic and postnatal development (Bandler et al., 2017; Hu et al., 2017b; Miyoshi, 2019; Williams and Riedemann, 2021). The MGE gives rise to distinct, non-overlapping interneuron subtypes, parvalbumin- (PV+) and somatostatin-expressing (SST+) interneurons (fast-spiking (FS) and non-fast spiking (NFS) interneurons, respectively).

Several factors regulate initial interneuron fate decisions within the MGE, including gradients of diffusible cues, spatial location of progenitors, temporal birthdates and the mode of neurogenesis (Flames et al., 2007; Glickstein et al., 2007; Wonders et al., 2008; Inan et al., 2012; Petros et al., 2015; Bandler et al., 2017, 2021; Mi et al., 2018; Allaway et al., 2021). The advent of single cell sequencing technologies over the last decade has generated a transcriptional and epigenetic ‘ground truth’ in the ganglionic eminences in mice (Mayer et al., 2018; Mi et al., 2018; Allaway et al., 2021; Bandler et al., 2021; Lee et al., 2022a; Rhodes et al., 2022), and more recently, in primates and humans (Nowakowski et al., 2017; Eze et al., 2021; Shi et al., 2021; Yu et al., 2021; Braun et al., 2022; Schmitz et al., 2022; Velmeshev et al., 2022; Zhao et al., 2022). With this baseline in place, researchers can better characterize how genetic and epigenetic perturbations affect the fate and maturation of GABAergic interneurons.

Epigenetic mechanisms play critical roles in gene expression during neurogenesis, and modifications of the chromatin landscape regulate cell state changes during neurodevelopment (Hirabayashi and Gotoh, 2010; Yao et al., 2016; Podobinska et al., 2017; Albert and Huttner, 2018). Alterations in epigenetic regulation can be associated with numerous neurodevelopmental disorders (Sokpor et al., 2017; Dall’Aglio et al., 2018; Wang et al., 2018). Enhancer of Zeste Homolog 2 (Ezh2) is the primary methyltransferase component of the Polycomb Repressive Complex 2 (PRC2) that is critical for trimethylation of histone 3 at lysine 27 (H3K27me3) resulting in gene repression (Muller, 1995; Cao et al., 2002; Aranda et al., 2015). Ezh2 is an evolutionary conserved gene that is aberrantly overexpressed in several forms of cancerous tumors (Sun et al., 2022; Zeng et al., 2022). Ezh2 expression is enriched in neural precursor cells where it represses target genes crucial to cell fate decisions and, in concert with other epigenetic marks, generates a transcriptional memory of specific gene expression patterns through cell divisions (Rice and Allis, 2001; Zaidi et al., 2011; Aranda et al., 2015). EZH2 variants can lead to Weaver Syndrome, a complex disease with variable degrees of intellectual disability (Tatton-Brown et al., 2013; Cohen et al., 2016), and dysregulation of H3K27me3 may be the primary driver in ataxia-telangiectasia (Li et al., 2013). Loss of Ezh2 can lead to ectopic exiting of the cell cycle and premature neuronal differentiation (Pereira et al., 2010; Zemke et al., 2015; Kim et al., 2018; Buontempo et al., 2022; Wu et al., 2022), neuronal migration defects (Di Meglio et al., 2013; Zhao et al., 2015), altered neuronal fate (Zemke et al., 2015; Feng et al., 2016; Wever et al., 2019) and changes in neuronal morphology and cognitive defects (von Schimmelmann et al., 2016; Zhang et al., 2019). Thus, Ezh2 is an important player in epigenetic regulation of neuronal fate and maturation, but a role for Ezh2 in forebrain GABAergic interneurons has not been explored.

We generated conditional Ezh2 knockout (cKO) mice to remove Ezh2 from the MGE and observed an increase in SST+ and decrease in PV+ interneurons across multiple brain regions. These fate changes were due to Ezh2 loss in cycling neural progenitors, as removing Ezh2 in post-mitotic cells did not alter interneuron fate. While the intrinsic physiology of MGE-derived interneurons in Ezh2 cKO mice was normal, fast-spiking cells displayed an increase in axonal length and branching. Fewer cortical MGE-derived interneurons were observed during the first postnatal week, which likely indicates decreased neurogenesis compared to WT mice. Single nuclei transcriptome analysis revealed an increase in SST expression and a decrease in genes predictive of PV-fated cells in the MGE of Ezh2 cKO mice, consistent with the fate changes observed in the adult. Lastly, while a global downregulation of H3K27me3 was observed in the Ezh2 cKO MGE, we observed that specific genomic loci were more susceptible to loss of Ezh2 while other loci were more resistant, indicating that global loss of Ezh2 had a differential effect of H3K27me3 at specific loci. In sum, we demonstrate that loss of Ezh2 disrupts H3K27me3, alters gene expression and cell proliferation in the MGE, which in turn disrupts the normal balance of SST+ and PV+ interneurons in the forebrain.

2 Materials and methods

2.1 Animals

All experimental procedures were conducted in accordance with the National Institutes of Health guidelines and were approved by the NICHD Animal Care and Use Committee (protocol #20-047). The following mouse lines were used in this study: Nkx2-1-Cre (Jax# 008661) (Xu et al., 2008), Ezh2^F/F (Jax# 022616) (Shen et al., 2008), Dlx5/6-Cre (Jax# 008199) (Monory et al., 2006) and Ai9 (Jax# 007909) (Madisen et al., 2010). Ezh2 genotyping was performed as previously described (Shen et al., 2008). For timed matings, noon on the day a vaginal plug was observed was denoted E0.5. Both embryonic and adult male and female embryonic mice were used without bias for all experiments.

2.2 Harvesting and fixing brain tissue

MGE dissections: E12.5 and E15.5 embryos were removed and placed in ice-cold carbogenated ACSF. Tails were clipped for PCR genotyping. During genotyping, embryonic brains were harvested and MGEs were dissected from each embryo and stored in ice-cold carbogenated artificial cerebral spinal fluid (ACSF, in mM: 87 NaCl, 26 NaHCO₃, 2.5 KCl, 1.25 NaH₂PO₄, 0.5 CaCl₂, 7 MgCl₂, 10 glucose, 75 sucrose, saturated with 95% O₂, 5% CO₂, pH 7.4). For E12.5, the entire MGE was removed as previously described (Lee et al., 2022b). For E15.5, the MGE was dissected under a fluorescent dissecting microscope to ensure collection of only Tom+ MGE cells and minimize collection of post-mitotic MGE-derived cells in the striatum anlage and other ventral forebrain structures. Upon obtaining genotyping results (∼90 min), MGEs from embryos of the same genotype were combined (when applicable) to generate single nuclei dissociations.

Postnatal brain fixations: All mice ≥ P5 were terminally anesthetized with an i.p. injection of Euthasol (270 mg/kg, 50 μl injection per 30 g mouse) and perfused with 4% paraformaldehyde (PFA). In some cases, tail snips were collected for genotyping prior to perfusion. Brains were removed and post-fixed in 4% PFA O/N at 4°C, then cryoprotected in 30% sucrose in PBS O/N at 4°C before embedding in OCT. Brains were sectioned at 30 μm on a CryoStar™ NX50 cryostat and stored as floating sections in antifreeze solution (30% ethylene glycol, 30% glycerol, 40% PBS) at −20°C in 96-well plates.

Embryonic brain fixations: Pregnant dams were terminally anesthetized with an i.p. injection of Euthasol (270 mg/kg, 50 μl injection per 30 g mouse). E12.5-E15.5 embryos were removed and placed in ice cold carbogenated ACSF. Embryonic brains were removed and incubated in 4% PFA O/N at 4°C. In some cases, tail snips were collected for genotyping. Brains were washed in PBS, transferred to 30% sucrose in PBS and embedded in OCT upon sinking. Brains were sectioned at 14–16 μm, mounted directly onto Permafrost slides, and stored at −80°C.

2.3 Immunohistochemistry and fluorescent in situ hybridizations (FISH)

30 μm free floating sections from P30-40 brains were washed in PBS to remove antifreeze solution and incubated in blocking solution (10% Normal Donkey Serum in PBS + 0.3% Triton X-100) for 1–2 h. Sections were incubated with primary antibodies in blocking solution for 48 h at 4°C, then washed in PBS for 2–4 h at RT. Sections were incubated with secondary antibodies with DAPI in blocking solution O/N at 4°C, washed in PBS, mounted and imaged.

Cryosectioned E12.5-E15.5 brains sections were incubated with blocking buffer (10% Normal Donkey Serum in PBS + 0.1% Triton X-100) for 1 h at RT, then incubated in primary antibody solution in blocking buffer O/N at 4°C, washed in PBS for 1–2 h, then incubated with secondary antibodies with DAPI in blocking solution for 2 h at RT or O/N at 4°C. Slides were washed in PBS and imaged. The following antibodies were used in this study: rabbit-anti H3K27me3 (1:100, Cell Signaling 9733T), rat-anti SST (1:300, Millipore MAB354), goat-anti PV (1:1000, Swant PVG213), rabbit-anti PV (1:1000, Swant PV27), rabbit-anti nNos (1:500, Millipore MAB5380), rabbit anti-Olig2 (1;500, Millipore AB9610). Species-specific fluorescent secondary antibodies used were conjugated to AlexaFluor^® 488, 647 and 790, and all used at 1:500.

RNAscope ISH assays™ (Advanced Cell Diagnostics) were performed on E12.5-E13.5 brain sections according to the manufacturer’s instructions. The following probes were used in this study: Nkx2-1 (434721), Ezh2 (802751-C3), Ezh1 (895231) and tdTomato (317041-C2).

All images were captured at 20X on a Zeiss Axioimager.M2 with ApoTom.2 (with Zen Blue software) or an Olympus VS200 Slide Scanner (VS200 ASW). Image post-processing was performed with Adobe Photoshop and ImageJ software.

2.4 Western blots and analysis

Core histone proteins from E13.5 MGE samples (pooled from 2 to 4 brains per genotype) were extracted using EpiQuik Total Histone Extraction kit (Epigentek# OP-0006). We obtained 40–80 μg of histone proteins per extraction and loaded 20 μg of protein onto a 4–12% Bis-Tris Plus mini gel (Invitrogen# NW04120BOX). Gels were run for ∼20 min at 200V using the Invitrogen mini gel tank with Blot MES SDS Running buffer. Gels were transferred to iBlot2 polyvinylidene difluoride (PVDF) membranes (Invitrogen# IB24001) in the iblot2 at 20V for 7 min. Blots were incubated with the primary antibody, mouse anti-H3 (Cell Signaling Technology# 3638S; 1:1000) and Rabbit anti-H3K27me3 (Cell Signaling Technology# 9733S; 1:500) overnight at 4°C, and then incubated with secondary anti-mouse-Starbright Blue-520 (Bio-Rad# 64456855; 1:2000) and anti-Rabbit-Starbright Blue-700 (Bio-Rad# 64484700; 1:2000) for 1 h at room temperature. Blots were imaged on a Bio-Rad ChemiDoC MP imaging system.

For analysis, a box was drawn around the ROI in each lane, with the same sized box used for both H3 and H3K27me3 signals for the WT, Het and cKO lanes in each blot. The average (mean) gray value was calculated in each box, then a lane-specific background signal taken just below each ROI was subtracted from each value. For normalization, each H3K27me3 value was divided by the corresponding H3 value for each lane (e.g., WT H3K27me3/WT H3). Then these Het and cKO values were divided by the WT value to determine the % of H3K27me3 signal compared to WT.

2.5 Cell counting

Adult brains: All cell counts were performed by hand and blind to genotypes. Total brains counted for Nkx2-1-Cre;Ezh2;Ai9 mice: WT = 5, Het = 5, cKO = 6, from 4 different litters. Total brains counted for Dlx5/6-Cre;Ezh2;Ai9 mice: WT = 4, Het = 3, cKO = 5, from 3 different litters. Counted cells consisted of either Tom+, Tom+/SST+, or Tom+/PV+ (and in the hippocampus, Tom+/nNos+); any SST+, PV+ or nNos+ cells that were Tom- were excluded from counts since they likely did not recombine at the Ezh2 (or Ai9) locus. For all sections, area was calculated using “Measurement” function in Photoshop, and all average areas described below include WT, het and cKO brains. Cortex: Counts were performed in somatosensory cortex on 3 non-consecutive sections per brain, then averaged together. Individual cortical images were divided into superficial and deep sections using DAPI staining to define the layer III-IV boundary. Average cortical area/section = 0.85 mm². Striatum: Counts were performed on 3 sections per brain, one section each through the anterior, middle and posterior striatum, then averaged together. Average striatal area/section = 3.15 mm². Hippocampus: Counts were performed on 8 non-consecutive sections per brain, then averaged together. More hippocampal sections were counted per brain due to the comparatively low number and section-to-section variability of Tom+ cells in the hippocampus. Sections were restricted to the anterior and middle hippocampus; the posterior hippocampus was excluded due to greater variability in interneuron density in this region. Hippocampal sections were divided into CA1, CA2/3 and DG regions using DAPI staining. Small Tom+ cell bodies (identified as oligodendrocytes) in CA2/3 were excluded from interneuron counts and counted as separate group. Average CA1 area/section = 0.98 mm², average CA2/3 area/section = 0.49 mm², average DG area/section = 0.60 mm².

P5 brains: All cell counts were performed by hand and blind to genotypes. Total brains counted for Nkx2-1-Cre;Ezh2;Ai9 mice: WT = 4 and cKO = 4, from 4 different litters. Counts were performed in somatosensory cortex on 3 non-consecutive sections per brain, then averaged together. Individual cortical images were divided into superficial and deep sections using DAPI staining. Average cortical area/section = 0.92 mm².

2.6 In vitro electrophysiology

Slice preparation: P30-P40 Nkx2-1-Cre;Ezh2 WT and cKO mice were anesthetized with isoflurane [5% isoflurane (vol/vol) in 100% oxygen], perfused transcardially with an ice-cold sucrose solution containing (in mM) 75 sucrose, 87 NaCl, 2.5 KCl, 26 NaHCO₃, 1.25 NaH₂PO₄, 10 glucose, 0.5 CaCl₂, and 2 MgSO₄, saturated with 95% O₂ and 5% CO₂ and decapitated. Brain was rapidly removed from the skull and transferred to a bath of ice-cold sucrose solution. Coronal slices of 300 μm were made using a vibratome (Leica Biosystems) and were stored in the same solution at 35°C for 30 min and at room temperature (RT) for an additional 30–45 min before recording.

Electrophysiology: Whole-cell patch clamp recordings on tdTomato+ cells in cortical layers V/VI cells were performed in oxygenated ACSF containing (in mM) 125 NaCl, 2.5 KCl, 26 NaHCO₃, 1.25 NaH₂PO₄, 10 glucose, 2 CaCl₂ and 1 MgCl₂. The ACSF was equilibrated with 95% O₂ and 5% CO₂ throughout an entire recording session which typically lasted between 30 min to 1 h to ensure sufficient permeation of neurobiotin. Recordings were performed at 30°C–33°C. Electrodes (3–7 MΩ) were pulled from borosilicate glass capillary (1.5 mm OD). The pipette intracellular solution contained (in mM) 130 potassium gluconate, 6.3 KCl, 0.5 EGTA, 10 HEPES, 5 sodium phosphocreatine, 4 Mg-ATP, 0.3 Na-GTP and 0.3% neurobiotin (pH 7.4 with KOH, 280–290 mOsm). Membrane potentials were not corrected for the liquid junction potential. During patching, cell-attached seal resistances were > 1 GΩ. Once whole-cell configuration was achieved, uncompensated series resistance was usually 5–30 MΩ and only cells with stable series resistance (<20% change throughout the recording) were used for analysis. Data were collected using a Multiclamp 700B amplifier (Molecular Devices), low-pass filtered at 10 kHz and digitally sampled at 20 kHz, and analyzed with pClamp10 (Molecular Devices). To characterize the intrinsic membrane properties of neurons, hyperpolarizing and depolarizing current steps were injected at 0.1 Hz under current-clamp configuration.

Data analysis: All intrinsic properties were measured in current-clamp configuration and calculated from 800 millisecond-long current injections unless noted otherwise. The resting membrane potential (in mV) was measured with 0 pA current injection a few minutes after entering whole-cell configuration. All other properties were measured holding the cell at −70 mV. Input resistance (in MΩ) was calculated using Ohm’s law from averaged traces of 100 ms long negative current injections of −20 pA. Action potential (AP) threshold was calculated as the potential when voltage change over time was 10 mV/ms using the first observed spike. AP amplitude (in mV) was calculated as the time difference in potential from the spike peak to spike threshold. AP/spike half-width (in ms) was calculated as the difference in time between the ascending and descending phases of a putative spike at the voltage midpoint between the peak of spike and spike threshold. Adaptation ratio was calculated as the ratio of the number of APs in the number of spikes in the last 200 ms over the number of APs in the first 200 ms of a positive current injection that elicited approximately 20 Hz firing. Afterhyperpolarization (AHP) amplitude was calculated as the difference between AP threshold and the most negative membrane potential after the AP, measured on the response to the smallest current step evoking an AP (Rheobase). Membrane time constant (in ms) was determined from a monoexponential curve best fitting the falling phase of the response to a small hyperpolarizing current step.

2.7 CUBIC clearing and streptavidin staining

After performing electrophysiological recordings, brain slices were fixed in 4% PFA in 0.1M PB and kept overnight at 4°C and then kept in 20% sucrose (in PB). The brain slices were processed for CUBIC (Clear, Unobstructed Brain/Body Imaging Cocktails and Computational analysis) clearing (Susaki et al., 2014). Slices were first washed with 0.1M PB (3 times for 10 min) at RT, followed by immersion in CUBIC reagent 1 for 2 days at 4°C. After 2 days of incubation, slices were washed with 0.1M PB (4 times for 30 min) at RT to ensure complete removal of CUBIC reagent 1. Slices were then incubated in fluorophore-conjugated streptavidin (1:500; ThermoFisher Scientific) in 0.1M PB (0.5% TritonX-100) overnight at 4°C. Slices were subsequently washed with 0.1M PB (4 times, 30 min) at RT and mounted with CUBIC reagent 1. Filled neurons were imaged with a Nikon A1R confocal microscope. Z-stacked images (each stack 1 μm) were acquired with a 40X oil-immersion objective.

2.8 Sholl analysis

Neurobiotin-filled neurons were processed and reconstructed using the Neurolucida 360 software (NL360) (MBF Biosciences). Briefly, image stack files were converted to JPEG 2000 file format with a MicroFile+ software, and the converted images were loaded into the NL360 software package. Neurites of PV were then traced manually in a 2D environment. Among the traced neurites, dendrites were easily distinguishable from axons, whose extensive ramifications maintained a constant diameter and had varicosities. Branched structure and Sholl analysis were performed using the built-in functions of the Neurolucida explorer, in which a series of concentric spheres (10 μm interval between radii) were created around the middle point (soma of the traced neuron).

2.9 Generating single nuclei suspensions for CUT&Tag and multiome experiments

Single nuclei suspensions were prepared as previously described (Lee et al., 2022a,b) with slight modifications. CUT&Tag: MGEs were transferred to a 1 mL Dounce homogenizer containing DNA Lysis Buffer (10 mM Tris–HCl (Ph.7.4), 10 mM NaCl, 3 mM CaCl₂, 0.1% Tween-20, 1.5% BSA and 0.1% IGEPAL CA-630 in nuclease-free water, 1 mL per sample). Cells were dounced with pestle A and pestle B, ∼10–15 times each, and pipetted through a 40 μm filter onto a pre-chilled 50 mL conical tube on ice and wet with 1 mL of DNA Nuclei Wash Buffer [in mM: 10 Tris–HCl (Ph.7.4), 10 NaCl, 3 CaCl₂, with 0.1% Tween-20 and 1.5% BSA in nuclease-free water, 5 mL per sample]. Lysed nuclei suspensions were then passed through a pre-wetted filter, and dounce was rinsed with 1 mL Nuclei Wash Buffer and passed through filter. After filtering, nuclei suspension was divided into 2 pre-chilled 2 mL tubes and spun at 500 g for 5 min at 4°C. After removing supernatant, nuclei pellet was washed once with 1 mL Nuclei Wash Buffer and then with 1 mL of 1X CUT&Tag Wash Buffer (from Active Motif CUT&Tag IT Assay kit), spinning with same conditions above. CUT&Tag Wash Buffer was removed, leaving ∼30–50 μL in each tube, and solutions were triturated and combined. Nuclei concentration was determined on Countess II FL Automated Cell Counter. 100,000–125,000 nuclei were used for each CUT&Tag reaction using the CUT&Tag IT Assay Kit (Active Motif, #53610) per manufacturer’s instructions. Primary antibody was rabbit-anti H3K27me3 (Cell Signaling, #9733T, 1:50) or rabbit-anti H3K27ac (Cell Signaling, #8173; 1:50), and secondary antibody was guinea pig anti-rabbit (Active Motif, 105465 from CUT&Tag IT kit, 1:100).

10x Genomics Multiome kit: MGEs were transferred to a 1 mL Dounce homogenizer containing Multiome Lysis Buffer [in mM: 10 Tris–HCl (Ph.7.4), 10 NaCl, 3 CaCl₂, 1 DTT, with 0.1% Tween-20, 1.5% BSA and 0.1% IGEPAL CA-630 in nuclease-free water, 1 mL per sample]. Cells were lysed by douncing with pestle A and pestle B, ∼10–15 times each. A 40 μm filter was placed onto a pre-chilled 50 mL conical tube on ice and wet with 1 mL of Multiome Nuclei Wash Buffer [in mM: 10 Tris–HCl (Ph.7.4), 10 NaCl, 3 CaCl₂, 1 DTT, with 0.1% Tween-20, 1.5% BSA and 1 U/μL in nuclease-free water, 5 mL per sample]. Nuclei suspension was then passed through pre-wetted filter, and dounce was rinsed with 1 mL Multiome Wash Buffer and passed through filter. After filtering, nuclei suspension was divided into 2 pre-chilled 2 mL tubes and spin at 500 g for 5 min at 4°C. After removing supernatant, nuclei pellet was washed twice with 1 mL Multiome Wash Buffer and spun as above. Multiome Wash Buffer was removed, leaving ∼20–30 μL of solution in each tube. Solution was then triturated to dissociate pellets and nuclei suspensions combined into 1 tube. Nuclei concentration was determined on Countess II FL Automated Cell Counter. Nuclei suspensions were diluted to ∼3,000–4,000 nuclei/μL, with 5 μL being used for 10x Genomics Multiome kit per manufacturer’s instructions. E15.5 data: MGE from 4 WT, Het and cKO embryos from 2 different E15.5 litters were combined to generate 1 biological rep. E12.5 data: MGE from 1 WT, Het and cKO mouse from a single E12.5 litter were used to generate 1 biological rep. Total cell numbers that passed QC: E12.5 WT = 6,391; E12.5 Het = 6,608; E12.5 cKO = 8,546; E15.5 WT = 11,477; E15.5 Het = 10,027; E15.5 cKO = 8,607.

2.10 CUT&Tag sequencing and analysis

Following library amplification, DNA quantity was determined with a Qubit and library quality characterized with an Agilent Tapestation. Libraries were balanced for DNA content and pooled before performing a final SPRI select bead 1x left size selection and paired-end sequenced (50 x 50 bp) on an Illumina NovaSeq. Sequencing read depths per library are detailed in Supplementary Table 5.

Paired-end run Illumina NovaSeq produced 2 compressed fastq.gz files for each replicate. For each age (E12.5 and E15.5) and each genotype, a total of 3 different biological reps from 3 different experiments were combined. Adaptor sequences were removed using cutadapt (Martin, 2011) v3.4 with the following parameters: -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT -q 20 –minimum-length 25. Trimmed reads were mapped to mouse reference genome (GRCm38/mm10) using bowtie2 (Langmead and Salzberg, 2012) v2.4.2 with the following parameters: –no-unal -N 1 –no-mixed –no-discordant –very-sensitive-local –local –phred33 -I 10. Aligned reads in sam files were further processed to remove multimappers if MAPQ was less than 10 (-q 10) and then sorted using samtools (Danecek et al., 2021) v1.12. Aligned reads that intersected blacklist regions (Amemiya et al., 2019) were removed and saved to bam files using bedtools (Quinlan and Hall, 2010) v2.30.0. PCR duplicates were removed from the bam files using picard MarkDuplicates v2.25.2.¹ Bigwig files were created from bam files using deepTools (Ramirez et al., 2016) (v3.5.1) bamCoverage with the following parameters: –binSize 5 –normalizeUsing RPKM. Peaks/domains were called using Epic2 (Stovner and Saetrom, 2019) on each bam file with the parameters –genome mm10 –guess-bampe and saved to bed files. Peak widths ranged from 200 bp to 261 kb, with a mean of 5.1 kb and median of 3.2 kb. Peaks were visualized using the Integrative Genomics Viewer (IGV) (Robinson et al., 2011) v2.16.1 by importing bigwig and bed files, respectively.

2.10.1 Differential analysis

Differentially enriched peaks were analyzed using edgeR (Robinson et al., 2010) v3.32.1 implemented in DiffBind (Ross-Innes et al., 2012) v3.0.15 in R v4.0.3.² Similarity in raw peak profiling across the samples was analyzed through PCA and a sample correlation heatmap using DiffBind::plotPCA and DiffBind::plot, respectively. For differential testing, CUT&Tag reads were counted in consensus peaks with default width in DiffBind::dba.count. The counts were subsequently normalized using the TMM (Robinson and Oshlack, 2010) method by setting the normalize argument to DBA_NORM_TMM in DiffBind::dba.normalize. In the TMM method, raw CUT&Tag counts, which were aggregated per consensus peak, were normalized using normalization factors that represented library size for moderately counted peaks. To determine normalization factors, peak-wise proportions of peak abundance in log2 scale were calculated between cKO and WT for individual consensus peaks. These proportions are defined as M values. M values calculated from peaks with extreme abundance were subsequently removed. A normalization factor for a sample is the weighted mean of the trimmed M values, where the weight is the inverse of the consensus peak’s approximate variance. The false discovery rate (FDR) was determined using the Benjamini–Hochberg (BH) (Benjamini and Hochberg, 1995) method for multiple hypothesis testing. Peaks with an FDR < 0.1 were considered statistically significant. Given the range of peak sizes, multiple settings for the “summit width” of DiffBind were used: (1) the default of using only reads in the 400 bp around the highest detected position within a peak; (2) extending this width to 2kb; and (3) disabling any summit to use reads from the whole peak/domain. Differential binding results differed slightly but were largely concordant across parameter choices, and the default of 400 bp was selected for further analysis because it best matched biological expectations.

Peaks were annotated to the nearest genes using ChIPseeker package (Yu et al., 2015) in R v4.0.3. For the TxDb and annoDb arguments in ChIPseeker::annotatePeak, we used TxDb.Mmusculus.UCSC.mm10.knownGene ( doi: 10.18129/B9.bioc.TxDb.Mmusculus.UCSC.mm10.knownGene) and org.Mm.eg.db ( doi: 10.18129/B9.bioc.org.Mm.eg.db), respectively. The transcript start site (TSS) region was defined as ranging from −5kb to +5kb. Statistically significant loci were visualized in MA plots (also known as Bland-Altman plots) using ggplot2 v3.3.3 in R.

2.11 Multiome sequencing and analysis

Joint libraries for snRNA-seq and snATAC-seq were created using 10x Genomics Single Cell Multiome ATAC + Gene Expression kit (1000285) by following manufacturer’s protocols. Sequencing was conducted with paired-end (50 x 50 bp) using an Illumina HiSeq 2500 (E12.5 snATAC-seq) or NovaSeq 6000 (E12.5 snRNA-seq, E15.5 snATAC-seq, E15.5 snRNA-seq), with detailed sequencing depth data in Supplementary Table 6.

The Cell Ranger ARC (v2.0.0) pipeline was used to process sequenced libraries with default parameters unless otherwise noted. Demultiplexed FASTQ files were generated by cellranger-arc mkfastq from BCL files. Reads were aligned to custom-built mouse (GRCm38/mm10) reference genome modified to include tdTomato using cellranger-arc count. Reads with de-duplicated and valid cell barcodes were used to build gene-by-barcode (snRNA-seq) and peak-by-barcode (snATAC-seq) matrices by cellranger-arc count per genotype. Individual matrices were aggregated to a single feature-barcode matrix file containing every genotype using cellranger-arc aggr without depth normalization (–normalize = none).

2.11.1 snRNA-seq data analysis

Seurat: An aggregated feature-barcode matrix was used as input to Seurat (Satija et al., 2015) (v4.0.5³) in R (v4.1.1, see text footnote 2). After imputing missing values to zero in metadata, outlier removal was performed on the number of counts per gene and percent reads mapping to mitochondrial genome (mitochondrial percentage). Lower limits for the number of counts per gene and mitochondrial percentage were set to 100 counts per gene and three standard deviations (SD) below the mean, respectively. Upper limits were set to three SD above the mean for both metrics. Negative datapoints created by subtraction of three SD from the mean were reset to 1, while datapoints that exceeded the upper limits were reset to the maximum datapoint. Cells were removed if they were more extreme then the upper/lower limits, or if they were eliminated from the snATAC-seq dataset during QC. The numbers of cells passing QC were: WT E12.5: 6,391; WT E15.5: 11,477; Het E12.5: 6,608; Het E15.5 10,027; cKO E12.5: 8,546; cKO E15.5: 8,607. Remaining cells were proceeded to the normalization workflow using Seurat::SCTransform using default parameters. For integration of snRNA-seq datasets from E12.5 and E15.5, 3,000 variable features were found using Seurat::SelectIntegrationFeatures on SCTransformed data. Prior to integration, anchors were identified using Seurat::FindIntegrationAnchors with the parameters dims, anchors.features, and normalization.method set to 1:30, the 30,000 variable features, and SCT, respectively. The integration was performed using Seurat::IntegrateData with identical dims and normalization.method to those from Seurat::FindIntegrationAnchors, along with the computed anchors. Dimensionality reduction was performed using Seurat::RunPCA and Seurat::RunUMAP with the parameters dims and umap.method set to 1:25 and uwot, respectively, on the integrated data.

2.11.2 snATAC-seq data analysis

Signac: An aggregated peak-by-barcode matrix was used as input to Signac (Stuart et al., 2021) (v1.4.0⁴) pipeline in R (v4.1.1). After imputing missing values in metadata as described above, outlier removal was performed on the number of chromatin accessibility peaks, transcription start site (TSS) enrichment score, and mitochondrial percentage. Lower limits were set to 1,000 counts per feature for the number of peaks and 2 for TSS enrichment score. Lower limit for mitochondrial percentage and all the upper limits were determined as described in the snRNA-seq analysis. Cells were removed if they were more extreme than the lower/upper limits in individual metrics, or if they were eliminated in corresponding snRNA-seq dataset during QC. Peaks were normalized via Term Frequency–Inverse Document Frequency (TF-IDF) method using Signac::RunTFIDF with default parameters. Highly variable peaks were found by Signac::FindTopFeatures with the min.cutoff set to 10. Dimensionality reduction was performed using Signac::RunSVD with default parameters. For integration of snATAC-seq datasets from E12.5 and E15.5, anchors were identified using Seurat::FindIntegrationAnchors with the parameters anchor.features and dims set to all features and 2:30, respectively. Integration of two snATAC-seq datasets was conducted using Seurat::IntegrateEmbeddings taking advantage of reciprocal LSI projection (RLSI), as instructed by Signac standard workflow. The Seurat::IntegrateEmbeddings ran using the previously computed anchors and 1 to 30 dimensions.

2.11.3 Multimodal analysis

Seurat: Multimodal data integration was started by finding Weighted Nearest Neighbor (WNN) (Hao et al., 2021) using Seurat::FindMultiModalNeighbors on snRNA-seq and snATAC-seq datasets with or without age-specific integration along with reduction lists set to Principal Component Analysis (PCA) (1–25 dimensions) for snRNA-seq and RLSI (1–20 dimensions) for snATAC-seq. Subsequent dimensionality reduction was performed using Seurat::RunUMAP on weighted.nn with default parameters. Cells were clustered on the WNN graph with Leiden algorithm and resolution 0.8 using Seurat::FindClusters. To visualize the DEGs between populations subsetted using Seurat, we employed the EnhancedVolcano tool. For significant determination, parameters were set to pCutoff = 10e^–6 and FCcutoff = 0.2. By default, the Uniform Manifold Approximation and Projection (UMAP) from the age-integrated multimodal analysis was utilized throughout the study. The expression of Sst and Pde1a genes in was examined on a multimodal UMAP at E15.5 without age integration.

For QC, the TSS enrichment score and the approximate ratio of mononucleosomal to nucleosome-free fragments (nucleosome signal) were computed using the functions Signac::TSSEnrichment and Signac::NucleosomeSignal with default parameters, respectively. The mitochondrial percentage was computed using the function Seurat::PercentageFeatureSet, which matches the pattern of gene names to “^[Mm][Tt]-”.

2.11.4 Differential analysis via DA-seq

Single cell DA analysis was conducted using DA-seq (Zhao et al., 2021) (v1.0.0⁵) in R (v4.1.1). The DA-seq used a PCA matrix as input, which was computed by Seurat after the multimodal integration of age-integrated snRNA-seq and snATAC-seq datasets. Cells were projected onto 2D space based on UMAP coordinates, which were computed by Seurat after multimodal integration of age-integrated snRNA-seq and snATAC-seq datasets. DA cells were determined by running DAseq::getDAcells, with the k.vector parameter set to every 50 between 50 and 500, and DAseq::updateDAcells with the pred.thres parameter set to +/-0.7. DA-seq conducted a random permutation test on abundance scores, using a threshold of +/- 0.7, to identify cells with an abundance score greater than 0.7 or less than −0.7.

2.11.5 Visualization of Multiome data

Uniform Manifold Approximation and Projection coordinates and WNN clustering, computed by Seurat on multimodal integrated datasets with or without age-specific integration, were imported to the Loupe Browser (v6.0.0, 10x Genomics). The expression of genes of interest and multimodal clustering were visualized using the Loupe Browser. DA analysis was visualized on a UMAP created using the functions DAseq::getDAcells and DAseq::updateDAcells in R, which take advantage of ggplot2 (v3.3.5). UMAPs visualizing the genotype and age in the DA analysis were generated using the Seurat::DimPlot. For QC, the following metrics were extracted from the metadata of the Seurat object: the number of ATAC fragments (nCount_ATAC), nucleosome signal (nucleosome_signal), the TSS enrichment score (TSS.enrichment), the number of RNA reads (nCount_RNA) and mitochondrial percentage (percent.mt). These metrics were then visualized using violin plots with the ggplot2 package in R.

To infer the developmental trajectories of E12.5 and E15.5 MGE cells across all conditions, Slingshot (Street et al., 2018) was used to generate trajectories onto UMAP projections. The Slingshot function was applied to the WNN matrix as described above, enabling the identification of neural developmental progression along the inferred trajectories.

2.12 Statistics and reproducibility

Cell counts: All cell counts were performed by hand and blind to genotype. Number of brain sections per mice and mice per genotype are described above and in figure legends. One-way ANOVA was used to compare WT, Het and cKO for all brain regions, followed by Tukey’s Multiple Comparison Test to identify significant differences between conditions. All statistical analysis was performed on GraphPad Prism (version 9.4.1). All raw cell counts, ANOVA F- and P-values, and results of Tukey’s multiple comparison tests related to Figures 2–4, 6, Supplementary Figures 2–4, 6 can be found in the Source Data file.

Electrophysiology: Collection of data and analysis were not performed blind and were non-randomized. Data from electrophysiological recordings are presented throughout as mean ± s.e.m unless otherwise noted. Unless indicated, all statistical comparisons were non-parametric. Number of recorded neurons (n) and number of animals (N) used were reported for each figure. All data were analyzed using pClamp10, GraphPad Prism and Neurolucida 360 software.

Sample size determination: No statistical method was used to predetermine sample size. For Western Blots, we desired n ≥ 2 biological replicates for each genotype to determine the percent of decreased H3K27me3 signal in the Ezh2 cKO. For cell counts with Nkx2-1-Cre mice, we wanted n ≥ 5 mice per genotype from at least 3 different litters. For cell counts with Dlx5/6-Cre mice, we wanted n ≥ 4 WT and cKO mice per genotype from at least 3 different litters. For electrophysiological recordings, we wanted to record from ≥ 10 neurons for each condition. For the single nuclei Multiome experiments, we wanted a minimum of 5,000 high-quality sequenced nuclei per condition (age and genotype), which should be sufficient to identify significant differences between conditions. This goal required ∼15,000 nuclei input for each sample (with the expectation of recovering ∼30–70% of nuclei/cells for each reaction based on 10x Genomics recommendations and our previous experience). Viable nuclei that passed QC ranged from 6,391 to 11,477 nuclei per condition (see above). Per standard single nuclei sequencing protocols, nuclei that did not pass stringent QC measurements (nCount_RNA and % mitochondria reads for snRNA; nCount_ATAC, nucleosome_signal and TSS enrichment for snATAC) in the Multiome datasets were excluded from analysis (as detailed in Supplementary Figure 7). For the CUT&Tag experiments, we strove for 100,000 nuclei for each reaction (actual range from 90,000 to 120,000 per reaction), with n = 3 biological replicates for each condition (age and genotype). All computational and statistical analysis are discussed in detail above and/or in the legends of the relevant figures and tables. All attempts at replication were successful.

3 Results

3.1 Loss of Ezh2 alters MGE-derived interneuron fate in the cortex and striatum

To characterize the function of Ezh2 in MGE-derived interneurons, we crossed Nkx2-1-Cre^C/C;Ezh2^F/+ males with Ai9^F/F;Ezh2^F/+ females to generate Nkx2-1-Cre^C/+;Ezh2;Ai9^F/+ wildtype, heterozygous and conditional knockout mice (hereafter WT, Het and cKO mice, respectively). In these mice, Ezh2 perturbation is restricted to MGE-derived cells in the telencephalon, and these cells also express tdTomato (Tom+). We first confirmed that Ezh2 expression is strongly downregulated in the MGE of cKO mice (Figure 1A). Since Ezh2 is critical for tri-methylation at histone H3K27, we also verified a significant reduction in H3K27me3 in the MGE in cKO mice (Figure 1B). To quantify this decrease, we performed Western Blots on MGE tissue from WT, Het and cKO mice. We observed a 15% decrease of H3K27me3 signal in Hets and a 47% decrease in cKO MGE compared to WT levels (Figure 1C). Some of this remaining H3K27me3 in the cKO MGE likely arises from the dorsal MGE where Cre expression is lacking in Nkx2-1-Cre mice (Xu et al., 2008), and from non-Nkx2-1 lineage cells within the MGE (endothelial cells, LGE- and CGE-derived cells migrating through this region, etc.). Additionally, the methyltransferase Ezh1 is ubiquitously expressed and can partially compensate for loss of Ezh2 in some conditions (Margueron et al., 2008; Shen et al., 2008). We do not observe significant changes in Ezh1 expression in the MGE of Ezh2 cKO mice (Supplementary Figure 1).

To determine if loss of Ezh2 affects interneuron fate, we counted the number of SST+ and PV+ interneurons in the somatosensory cortex, hippocampus and striatum from WT, Het and cKO mice. We observed a moderate but significant decrease in the density of Tom+ MGE-derived cortical interneurons in the cortex of cKO mice. More striking, there was a significant increase in the density of SST+ interneurons and a corresponding decrease in the density of PV+ interneurons in the cortex of cKO mice (Figures 2A, B). This shift in cell fate is also apparent when examining the proportion of Tom+ cells that express either SST or PV (Figure 2C).

To examine if these differences were consistent throughout cortical layers, we divided the somatosensory cortex into superficial (I-III) and deep (IV-VI) layers based on DAPI staining (Figure 2A). There was no significant decrease in the density of Tom+ cells between WT and cKO mice in either the superficial or deep layers (Supplementary Figure 2). The strongest effect in SST and PV fate changes was found in deep cortical layers, whereas superficial layers displayed a more moderate increase in SST+ and decrease in PV+ cells (Supplementary Figure 2). Thus, loss of Ezh2 in the MGE results in a slight reduction in total MGE-derived cortical neurons, with a significant increase in SST+ and decrease in PV+ cells, most notably in the deeper cortical layers.

MGE-derived SST+ and PV+ interneurons also populate the adult striatum. We observed a significant decrease in the density and proportion of PV+ interneurons in the striatum of cKO mice, but no change in the total number of Tom+ or SST+ cells (Supplementary Figure 3), indicating this decrease in the density of PV+ interneurons is observed in multiple brain regions.

3.2 Alteration in both hippocampal interneurons and oligodendrocytes in Ezh2 cKO

The hippocampus contains a population of MGE-derived, neuronal nitric oxide synthase nNos-expressing (nNos+) neurogliaform and ivy cells that are not found in the cortex (Tricoire et al., 2010, 2011; Jaglin et al., 2012). SST+, PV+ and nNos+ interneurons each make up ∼1/3 of MGE-derived hippocampal interneurons (Quattrocolo et al., 2017). We characterized the densities and percentages of SST+, PV+ and nNos+ interneurons throughout the hippocampus (Figures 3A, B). No significant differences were found in the densities of Tom+, SST+, PV+ or nNos+ interneurons in the whole hippocampus (Figure 3C). However, there was a significant increase in the percentage of Tom+ cells that expressed SST or nNos in cKO mice compared to WT (Figure 3D).

Since the prevalence of interneurons differs between the CA1, CA2/3 and dentate gyrus (DG) regions of the hippocampus, we divided hippocampal sections into these three regions based on DAPI staining (Figures 3A, B). There were almost no differences in densities or the proportion of subtypes in the CA1, only a slight but significant increase in the percentage of Tom+/nNos+ cells in the cKO (Supplementary Figure 4A). This increase in the percentage of nNos+ cells was also detected in the DG. More striking was a strong reduction of both PV+ cell densities and the percentage of Tom+/PV+ in the DG (Supplementary Figure 4B).

The CA2/3 region displayed the greatest differences between WT and cKO mice, many of which mimicked cell fate changes in the cortex. First, there was a decrease in the total density of CA2/3 Tom+ cells in the cKO compared to WT (Figure 3E). Second, there was a significant decrease in PV+ cell densities in the CA2/3 region of cKO mice, and a corresponding increase in the percentage of Tom+/SST+ (Figures 3E, F). Third, there was a significant increase in the density and percentage of nNos+ cells in CA2/3 region of cKO mice compared to WT (Figures 3E, F). Thus, loss of Ezh2 in the MGE had both broad and region-specific effects on interneuron fate in the hippocampus: an increase in the proportion of MGE-derived nNos+ cells in all three hippocampal regions, an increase in the percentage of Tom+/SST+ cells specifically in the CA2/3 region (which resulted in a significant increase in the entire hippocampus), and a decrease in PV+ cell density in both CA1 and CA2/3. This decrease in PV cells was observed in the cortex, striatum and hippocampus.

While performing cell counts in the CA2/3 region, we observed Tom+ cell bodies that were too small to be interneurons, and we did not observe these cells in other brain regions (Figure 4A). Counting these cells separately, we found a very strong reduction of these CA2/3-specific cells in the cKO hippocampus (Figures 4B, C). We stained WT hippocampal sections with various glia and microglia markers and found that many of these small Tom+ cell bodies were Olig2+, indicating that they are likely oligodendrocytes (Figure 4B). This decrease in oligodendrocytes in Ezh2 cKO mice is consistent with findings demonstrating that loss of Ezh2 can block or delay gliogenesis (Hirabayashi et al., 2009; Wang et al., 2020).

3.3 Normal intrinsic properties but altered morphology of Ezh2 cKO interneurons

To characterize the intrinsic physiology of MGE-derived interneurons in cKO mice, we performed patch clamp recording of layer V/VI Tom+ cortical cells in acute brain slices. Cells were classified as FS or NFS based on their intrinsic electrophysiological properties characterized under current clamp-recording. NFS cells had larger half-width, input resistance and membrane time constant/Tau, but smaller rheobase compared to FS cells. We analyzed action potential shapes, resting membrane potential, spike adaptation ratio, afterhyperpolarization (AHP) amplitude, input resistance and rheobase. There were no differences in intrinsic properties of FS or NFS cortical interneurons between WT and cKO mice (Supplementary Figure 5).

However, reconstructions of recorded cells did reveal morphological changes in FS cells. The axonal arbor of FS cells from cKO mice were larger and more complex compared to WT cells (Figure 5A). Sholl analysis revealed a significant increase in axon intersections, axon length and axon volume in FS cKO cells, while there were no changes in dendritic arbors (Figure 5B). This increased axonal arbor is similar to what was observed when trkB signaling was blocked in PV cells (Guyon et al., 2021). Thus, while intrinsic properties of Tom+ MGE-derived cortical interneurons were normal in Ezh2 cKO mice, FS cells displayed greater complexity in their axonal arbors compared to FS cells from WT mice.

3.4 Loss of Ezh2 in cycling progenitors is required for cell fate changes

We next wanted to determine at what stage of development loss of Ezh2 results in cell fate changes. Ezh2 is enriched in cycling progenitors throughout the embryonic brain (Figure 1A), but Ezh2 may play a critical at other stages as well. To investigate this possibility, we generated Dlx5/6-Cre;Ezh2^F/F;Ai9 conditional cKO mice in which loss of Ezh2 is restricted to post-mitotic neurons arising from the ganglionic eminences. We verified that Ezh2 is still expressed in MGE ventricular and subventricular zone cycling progenitors in Dlx5/6-Cre cKO mice (Supplementary Figure 6A). There were no differences in the densities or percent of SST+ and PV+ cells in the cortices of these Dlx5/6-Cre cKO mice (Supplementary Figure 6B), indicating that Ezh2 is required in cycling MGE progenitors for proper interneuron fate and maturation.

A wave of programmed apoptosis occurs between the first and second postnatal weeks that eliminates ∼20–40% of cortical interneurons (Southwell et al., 2012; Denaxa et al., 2018; Priya et al., 2018; Duan et al., 2020). To determine if there were changes in the overall production of MGE-derived interneurons during embryogenesis, we counted the number of Tom+ cells in the cortex at P5, prior to programmed apoptosis. We found a significant decrease in the number of Tom+ cells in the cKO cortex compared to WT at P5 (Figures 6A, B). This finding supports the hypothesis that loss of Ezh2 in cycling MGE progenitors decreases the overall production of MGE-derived cortical interneurons, which would lead to more prominent loss of PV+ interneurons due to their bias production at later embryonic timepoints compared to SST+ cells (Inan et al., 2012).

3.5 Changes in gene expression and chromatin accessibility in the MGE of Ezh2 cKO mice

We then investigated whether transcriptomic and epigenetic changes in the Ezh2 cKO MGE are predictive of these cell fate changes in the mature forebrain. We generated single nuclei suspensions from E12.5 and E15.5 MGE and used the 10x Genomics Multiome kit to define the gene expression profile and chromatin accessibility within individual cells. We obtained a total of 51,656 nuclei that passed QC (Supplementary Figure 7A), with a range of ∼6,300–11,500 nuclei per condition. These nuclei can be clustered based on their transcriptome, chromatin accessibility or integrated RNA and ATAC using the Weighted Nearest Neighbor (WNN) function in Seurat (Hao et al., 2021; Figures 7A, B; Supplementary Figure 7B). Analysis of several genes revealed the expected progression from radial glia cells (Nes) to immature neuronal progenitors (Ascl1) to post-mitotic neural precursors (Dcx and Rbfox3), which was verified with pseudotime analysis (Figure 7C).

Since there was a general increase in SST+ interneurons in the Ezh2 cKO mouse, we examined SST expression in this single nuclei dataset. While no obvious differences in SST expression was apparent at E12.5, we did observe an increase in SST expression in the MGE of E15.5 cKO compared to WT (Figure 8A). While PV is not expressed in the embryonic mouse brain, two genes that are enriched in PV-fated interneurons and critical for their development are the transcription factors Mef2c and Maf (Mayer et al., 2018; Pai et al., 2020). We found that both genes are strongly reduced in E15.5 cKO MGEs compared to WT (Figure 8B). Complementing this gene expression analysis, the motifs for these transcription factors are enriched in accessible regions of the E15.5 WT MGE compared to cKO mice (Figure 8C). Thus, our gene expression analysis reveals an apparent increase in SST expression and decrease in Maf and Mef2c expression in the MGE of E15.5 cKO mice, which is predictive for the increase in SST+ and decrease of PV+ interneurons in Ezh2 cKO brains.

Additionally, WT and cKO cells appeared to display differential abundance (DA) in specific regions of the UMAP plot, most notably in the E15.5 dataset (Figure 7B; Supplementary Figure 7B). To confirm this observation, we performed DA analysis using DA-seq (Zhao et al., 2021). DA-seq determines a DA score for each cell, whereby a cell that is surrounded by cKO cells in a k-nearest neighbor (KNN) graph has a score closer to +1, and a cell surrounded by WT cells has a score closer to −1. This DA score does not rely upon previously identified clusters, and it does not require similar cell numbers between different conditions. Our DA-seq analysis revealed that most cells with a DA score above +0.7 or below −0.7 were from the E15.5 MGE (Figure 9A), whereas cells from the E12.5 MGE displayed little differential abundance. Most cells with a DA score below −0.7 (blue, WT bias) were in the clusters enriched for Maf and Mef2c, which are putative PV+ interneurons, whereas the cluster containing SST+ cells contained many cells with DA score above 0.7 (red, cKO bias) (Figure 9A). Thus, we observe a differential abundance of WT and cKO cells specifically in the E15.5 dataset, with an increase of cKO cells in clusters expressing SST and an increase of WT cells in clusters where Maf and Mef2c are strongly enriched.

To specifically focus on these clusters enriched for SST- and PV-fated interneurons, we extracted the four clusters containing these cells from the dataset (clusters 3, 4, 5 and 7) (Figures 7A, 9B). Using thresholds of a Log₂ fold change (FC) > ± 0.2 and a false discovery rate (FDR) of -Log₁₀ Adjusted P-value (-Log₁₀Adj P) < 10^–6, we identified 66 differentially expressed genes at E12.5 (52 downregulated and 14 upregulated in cKO) and 174 differentially expressed genes at E15.5 (46 downregulated and 128 upregulated in cKO) (Figure 9C; Supplementary Tables 1–2). Notably, SST was significantly upregulated in the E15.5 cKO MGE whereas both Maf and Mef2c were significantly downregulated in the E15.5 cKO MGE, consistent with the observations above (Figure 9C). Additionally, the MGE-specific transcription factors Nkx2-1 and Lhx6 were also upregulated in the E15.5 cKO MGE.

Since differentially expressed genes of interest were restricted to E15.5, we reexamined the integrated E15.5 RNA+ATAC dataset alone (Supplementary Figure 7B). In this dataset, SST was enriched in cluster 10 (Figure 9D). The top gene enriched in this cluster was Phosphodiesterase 1A, Pde1a, which is also significantly upregulated in the E15.5 cKO MGE (Figure 9C). Based on the Allen Brain Institute’s single cell transcriptomic adult mouse brain dataset, Pde1a is enriched in many SST+ interneuron subtypes while it’s excluded from PV+ interneurons (Figure 9D), providing additional evidence that there is an increase in SST-fated cells in the MGE of Ezh2 cKO mouse. In sum, this analysis indicates that the shift in MGE-derived interneuron subtypes is most evident in the E15.5 MGE, with an increase in both SST expression levels and the number of SST-expressing cells in the MGE of Ezh2 cKO mice, with a corresponding decrease in PV-fated cells and PV-predictive genes.

3.6 Changes of H3K27me3 at distinct genomic loci in MGE of Ezh2 cKO mice

H3K27me3 levels are strongly downregulated in the MGE of Ezh2 cKO mice (Figure 1). This downregulation could occur homogenously throughout the genome, or alternatively, there could be differential changes where some loci are more resistant to loss of Ezh2 while other sites are more susceptible. For example, the promoter-enhancer interaction at the Sox2 locus is extremely resistant to epigenetic perturbations, possibly due to its critical function in in the epiblast stage of pre-implantation embryos (Chakraborty et al., 2023). To explore differential H3K27me3 changes at specific genomic loci, we performed bulk CUT&Tag (Kaya-Okur et al., 2019) with a H3K27me3 antibody in the MGE of WT and cKO mice. We conducted 3 biological replicates for each age and genotype, using Epic2 (Stovner and Saetrom, 2019) for peak calls and the DiffBind package (Stark and Brown, 2023) with edgeR’s trimmed mean of M values (TMM) (Robinson et al., 2010) for comparative analysis.

We hypothesized that H3K27me3 would not be uniformly lost throughout the genome in the MGE of Ezh2 cKO mice and that some loci might preferentially retain (or lose) H3K27me3. CUT&Tag experiments with global differences between conditions are challenging to interpret. Spike-ins can help, but the amount of added material must be known with high accuracy and precision to avoid propagating measurement error in spike-in concentration to the normalized CUT&Tag data. Without spike-ins, sequencing and library size normalization effectively equalizes any global differences between samples, and the resulting normalized data can be interpreted as relative differences (in this case H3K27me3 levels at specific genomic loci between WT and cKO mice). This interpretation enables us to examine the redistribution of H3K27me3 throughout the genome. We interpret normalized, equalized H3K27me3 signal that is lower in cKO compared to WT to indicate H3K27me3 was preferentially lost from that locus (i.e., a locus more susceptible to loss of Ezh2). Conversely, we interpret H3K27me3 signal that is higher in cKO compared to WT as preferentially retained at that locus (i.e., a locus more resistant to loss of Ezh2).

Biological replicates were grouped together in both principal component analysis (PCA) and unbiased correlation heatmaps, with age being a greater differentiation factor compared to genotype (Figures 10A, B). Comparative analyses revealed significant changes in the distribution of H3K27me3 at 82 loci at E12.5 (59 “decreased” and 23 “increased” loci in cKO) and 51 loci at E15.5 (14 “decreased” and 37 “increased” loci in cKO) (Figure 10C; Supplementary Tables 3–4). One of the most prominent changes was at the Foxp4 locus, which is very susceptible to Ezh2 loss with a strong ‘decrease’ in H3K27me3 at both E12.5 and E15.5 (Figures 10C, D). Despite this significant change, we do not observe a concomitant increase in Foxp4 gene expression in the cKO MGE (Figure 10D). While Foxp4 is strongly enriched in the LGE (Takahashi et al., 2008), its function in forebrain development has not been explored. Several genes critical for development of MGE-derived interneurons displayed a relative ‘increase’ in H3K27me3 (and thus increased resistance to loss of Ezh2) in the cKO MGE: Nkx2-1 and the Dlx1/2 locus were increased at E12.5, and Lhx6 and the Dlx5/6 locus were increased at E15.5 (Figures 10C, E–G). We did not observe a corresponding decrease in expression levels of these genes in the cKO MGE. In fact, several of these genes were upregulated in the E15.5 cKO MGE (Figure 9C), which is likely due to the absolute downregulation of H3K27me3 at these loci.

We identified 240 differentially expressed genes in the Multiome dataset and 133 peaks with significantly different relative H3K27me3 levels in the CUT&Tag dataset. Integrating these two datasets is challenging because (1) the CUT&Tag data represents ‘relative’ changes in H3K27me3 levels rather than global changes, and (2) we cannot definitively associate H3K27me3 peaks with specific genes. With these caveats, we identified 7 genes/loci that were significantly different between WT and Ezh2 cKO MGE in both datasets. Three genes displayed increased relative H3K27me3 levels and increased gene expression in the cKO MGE (Lhx6, Myt1l and Nkx2-1). It is intriguing that two of the most important genes for MGE development (Nkx2-1 and Lhx6) display resistance to changes in H3K27me3 levels and increased gene expression in Ezh2 cKO mice. One gene displayed decreased relative H3K27me3 levels and decreased gene expression in the cKO MGE (Zeb2). ZEB2 is co-expressed with MAF+ cells in the human MGE (Yu et al., 2021), so decreased of Zeb2 expression in the cKO MGE is another predictor of decreased PV+ cells. Three genes displayed increased relative H3K27me3 levels and decreased gene expression in the cKO MGE (Kirrel3, Meis1 and Zfhx3), which represents a logical correlation between H3K27me3 levels and gene expression. Variants in KIRREL3 and ZFHX3 have recently been associated with ASD and neurodevelopmental disorders, respectively (Taylor et al., 2020; Del Rocio Perez Baca et al., 2023), so it will be interesting to determine how these genes regulate fate of MGE-derived cells in the future. There were zero genes that displayed decreased relative K27me3 levels and increased gene expression in the cKO MGE.

Even though there is a strong global decrease of H3K27me3 in the Ezh2 cKO MGE (Figure 1), we identified shifts in the relative level of H3K27me3 at various loci in the genome. This indicates that specific genomic loci are more susceptible (e.g., Foxp4) or resistant (e.g., Nkx2-1) to H3K27me3 loss in the absence of Ezh2. The mechanism by which loss of Ezh2 generates these differential effects at genomic loci, and how these changes in H3K27me3 levels relate to gene expression, require further investigation.

4 Discussion

There is growing evidence that dysregulation of epigenetic mechanisms can lead to a variety of human diseases and neurodevelopmental disorders (Kalish et al., 2014; Sokpor et al., 2017; Dall’Aglio et al., 2018; Qureshi and Mehler, 2018; Wang et al., 2018). For example, postmortem tissue from schizophrenic patients displays alterations in genome organization and other epigenomic characteristics (Kozlenkov et al., 2018; Rajarajan et al., 2018). Additionally, many genes associated with neurological and psychiatric diseases are enriched in immature interneurons during embryonic development (Cross-Disorder Group of the Psychiatric Genomics Consortium, 2019; Grove et al., 2019; Schork et al., 2019; Trevino et al., 2020). Thus, advancing our knowledge of gene regulation mechanisms during interneuron development is critical for understanding both normal development and disease etiologies.

In this study, we find that loss of Ezh2 in the MGE decreases the density and proportion of PV+ cells, often with a corresponding increase in SST+ cells. This shift in interneuron fate was most prominent in the cortex and the CA2/3 region of the hippocampus, with an overall decrease in total MGE-derived interneurons also observed in the cortex. A decrease in PV+ cells was observed in the CA1, DG and striatum without a significant increase in SST+ cells (Figures 2–3; Supplementary Figures 2–4). In the hippocampus, we also observe an increase in MGE-derived nNos+ cells in the cKO (Figure 3). Unlike PV+ and SST+ interneurons, the spatial and temporal origin of hippocampal nNos+ in the MGE is not well characterized, so it’s unclear how loss of Ezh2 increased nNos+ cells. These phenotypes were due to Ezh2 function in cycling MGE progenitors, as no changes were observed when Ezh2 was removed in post-mitotic MGE cells (Supplementary Figure 6). This differential severity of interneuron fate changes in specific brain regions highlights the importance of examining interneuron fate in multiple brain regions.

In the hippocampus, CA2/3 displayed the strongest shift in interneuron fate compared to other hippocampal regions, more closely matching changes in the cortex. CA1 and CA2/3 have similar densities and percentages of MGE-derived interneurons (Figure 3; Supplementary Figure 4), so why these regions display such different phenotypes in the Ezh2 cKO remains unknown. The severity of CA2/3-specific changes could arise from differences in interneuron migration, cell survivability or specific circuit interactions within this region. Additionally, the high density of MGE-derived oligodendrocytes in the WT CA2/3 was quite striking, as we did not observe these cells in the cortex or other hippocampal regions. MGE-derived oligodendrocytes do migrate into the cortex, but they are almost entirely eliminated during the second postnatal week (Kessaris et al., 2006). We are unaware of any reports describing MGE-derived oligodendrocytes in the hippocampus, and specifically in the CA2/3 region. Further studies are needed to understand why a population of MGE-derived oligodendrocytes perdures in CA2/3 into adulthood. Assuming these MGE-derived oligodendrocytes are generated toward the end of the cell cycle, then loss of this population is consistent with preemptive cell cycle exit in the Ezh2 cKO mouse and the overall reduced production of MGE-derived cortical interneurons observed at P5 (Figure 6).

Previous studies revealed that the spatial origin and temporal birthdates of MGE-derived progenitors play a role in cell fate decisions. Cells born in the dorsal-posterior MGE preferentially give rise to SST+ interneurons whereas PV+ cells predominantly arise from the ventral-anterior MGE (Flames et al., 2007; Fogarty et al., 2007; Wonders et al., 2008; Inan et al., 2012; Taniguchi et al., 2013; McKenzie et al., 2019). Additionally, SST+ cells are preferentially born early (E11.5-E15.5) from apical progenitors whereas PV+ interneurons arise later (E12.5-E17.5) predominantly from basal progenitors (Inan et al., 2012; Taniguchi et al., 2013; Petros et al., 2015). As loss of Ezh2 has been shown to induce cell cycle exit and neurogenesis (Pereira et al., 2010; Zemke et al., 2015; Kim et al., 2018; Buontempo et al., 2022; Wu et al., 2022), this could in part explain the fate switch in Ezh2 cKO mice. Premature neurogenesis would generate an excess of neurons at early stages when SST+ interneurons are born, while simultaneously deplete the progenitor pool so that there are fewer neurogenic divisions at later stages when PV+ cells arise. This could also result in an overall decrease in MGE-derived interneurons, as was apparent at P5 (Figure 6). Manipulation of cell cycle dynamics has previously been shown to shift the percentage of SST+ or PV+ cells born at a particular timepoint (Inan et al., 2012; Taniguchi et al., 2013; Petros et al., 2015).

Despite these changes in cell fate, mature PV+ and SST+ cortical interneurons displayed normal intrinsic physiological properties. However, cortical PV+ interneurons displayed significantly greater axon length and complexity in the cKO mouse. This could be a form of compensation, as one way to increase inhibition with fewer PV+ cells is for the surviving PV+ interneurons to have increased synaptic contacts. Alternatively, enhanced activation of PV+ interneurons can result in more elaborate axon morphologies (Stedehouder et al., 2018): if individual PV+ interneurons are receiving more glutamatergic inputs because of decreased number of PV+ cells in Ezh2 cKO mice, than this could in part explain the morphological changes. Additionally, PV+ interneurons normally form exuberant axon projections and synaptic connections during the first two postnatal weeks which are then pruned and retracted during the 3rd and 4th weeks (Micheva et al., 2021). Thus, it’s possible that this normal developmental retraction is disrupted in the cKO mouse.

Single nuclei sequencing revealed that these cell fate changes are already apparent in the developing MGE. First, the majority of transcriptome differences were observed in the E15.5 MGE while the E12.5 MGE from WT and cKO mice were quite similar. This could imply that the loss of Ezh2 takes time to manifest (due to perdurance of H3K27me3 marks, rate of histone turnover, etc.), with the strongest affects becoming more evident during the second half of MGE neurogenesis. Second, there was an increase in SST expression and a decrease in expression of PV-predictive genes Mef2c and Maf in the Ezh2 cKO mouse at E15.5 (Figures 8, 9C). Third, differential abundance analysis reveals that there is an enrichment of cells from the cKO MGE in the SST-expressing region and an enrichment for WT MGE cells in the Mef2c/Maf-expressing region (Figure 9A). Fourth, the top gene expressed in the SST-enriched cluster at E15.5 is Pde1a, which is expressed by many mature SST+ interneuron subtypes but excluded from PV+ interneurons (Figure 9D). In sum, these transcriptional and cellular differences in the MGE drive the shifts in SST+ and PV+ interneurons in the adult brain of Ezh2 cKO mice.

Despite the global downregulation of H3K27me3, we found that specific genomic loci were either more resistant (“increased” H3K27me3) or susceptible (“decreased” H3K27me3) to loss of Ezh2. One of the most susceptible regions at both E12.5 and E15.5 with a strong “decrease” in H3K27me3 was the Foxp4 locus (Figures 10C, D). Foxp4 is enriched in the LGE (Takahashi et al., 2008), but it has not been well-studied in neurodevelopment. In heterologous cell lines, FOXP4 can directly interact with the transcription factors SATB1, NR2F1 and NR2F2, all of which are critical for development of MGE-derived interneurons (Close et al., 2012; Hu et al., 2017a; Estruch et al., 2018). Whether these interactions occur in the developing brain is unclear. A study on medulloblastoma found that Foxp4 and Ezh2 are both targets of the microRNA miR-101-3p (Xue et al., 2022), indicating the function of these genes might be linked in some scenarios. Why the Foxp4 locus is extremely sensitive to Ezh2 loss requires further study.

In the Ehz2 cKO MGE, we observed that several genes critical for fate determination of MGE-derived interneurons (Nkx2-1, Lxh6, Dlx1/2 and Dlx5/6) were more resistant to loss of Ezh2 (Figures 10C, E–G). This raises the possibility that some genes playing critical roles in fate determination may be more resistant to epigenetic changes. As Nkx2-1 is a ‘master regulator’ of MGE fate and the Nkx2-1 locus displays unique chromatin organization specific to the MGE (Rhodes et al., 2022), this locus could be more resistant to epigenetic modifications. The interaction between Sox2, a transcription factor essential in the epiblast of pre-implantation embryos, and a critical enhancer downstream is maintained even when artificial boundaries are introduced between these regions (Chakraborty et al., 2023). It will be interesting to explore if this theme of specific genomic loci being resistant to epigenetic modifications holds true at critical fate-determining genes in other tissues and/or stages of development, both in terms of normal development and with respect to neurodevelopmental and psychiatric diseases.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary material.

Ethics statement

The animal study was approved by the NICHD Animal Use and Care Committee (ACUC). The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

CR: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Visualization, Writing—original draft, Writing—review and editing. DA: Data curation, Formal analysis, Investigation, Methodology, Writing—review and editing. MS: Data curation, Formal analysis, Investigation, Methodology, Software, Visualization, Writing—review and editing. SN: Data curation, Formal analysis, Investigation, Methodology, Visualization, Writing—review and editing. LE: Data curation, Formal analysis, Investigation, Writing—review and editing. PS: Formal analysis, Investigation, Visualization, Writing—review and editing. DL: Formal analysis, Methodology, Visualization, Writing—original draft, Writing—review and editing. YZ: Investigation, Visualization, Writing—review and editing. PR: Funding acquisition, Methodology, Supervision, Writing—review and editing. RD: Data curation, Formal analysis, Funding acquisition, Methodology, Supervision, Writing—review and editing. SL: Methodology, Supervision, Writing—review and editing. TP: Conceptualization, Funding acquisition, Resources, Supervision, Visualization, Writing—original draft, Writing—review and editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by Eunice Kennedy Shriver NICHD Intramural Awards to TP, PR, and RD; NIMH Intramural Award to SL; NICHD Scientific Director’s Award to TP and PR.

Acknowledgments

We thank S. Coon, F. Faucz, J. Iben, T. Li, and other members of the NICHD Molecular Genomic Core; members of the Petros Lab for helpful discussions and feedback on the manuscript. Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, TP (tim.petros@nih.gov).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncel.2024.1334244/full#supplementary-material

Footnotes

1. ^ https://broadinstitute.github.io/picard/
2. ^ https://cran.r-project.org/
3. ^ https://satijalab.org/seurat/
4. ^ https://stuartlab.org/signac
5. ^ https://klugerlab.github.io/DAseq

References