Abstract
Epigenetic mechanisms are essential to establish and safeguard cellular identities in mammals. They dynamically regulate the expression of genes, transposable elements and higher-order chromatin structures. Consequently, these chromatin marks are indispensable for mammalian development and alterations often lead to disease, such as cancer. Bivalent promoters are especially important during differentiation and development. Here we used a genetic screen to identify new regulators of a bivalent repressed gene. We identify BEND3 as a regulator of hundreds of bivalent promoters, some of which it represses, and some of which it activates. We show that BEND3 is recruited to a CpG-containg consensus site that is present in multiple copies in many bivalent promoters. Besides having direct effect on the promoters it binds, the loss of BEND3 leads to genome-wide gains of DNA methylation, which are especially marked at regions normally protected by the TET enzymes. DNA hydroxymethylation is reduced in Bend3 mutant cells, possibly as consequence of altered gene expression leading to diminished alpha-ketoglutarate production, thus lowering TET activity. Our results clarify the direct and indirect roles of an important chromatin regulator, BEND3, and, more broadly, they shed light on the regulation of bivalent promoters.
Graphical Abstract
INTRODUCTION
The development and life of organisms depend on the execution of precise transcriptional programs, which allow the different cell types to respond to stimuli and fulfill their biological function. Part of the information underpinning these programs is hardwired in the genome: promoters and enhancers drive the context-dependent expression of genes, in response to transcriptional regulators that recognize specific sequences (1). In addition to this sequence-based information, chromatin constitutes an additional regulatory layer that is vital for proper gene expression. Nucleosomes are omnipresent in the genome, and the modifications borne by their histone tails profoundly influence gene expression. This regulation occurs through modulating chromatin compaction and directly recruiting or repelling nucleosome remodelers, as well as transcriptional activators and repressors (2). These chromatin marks can be remodeled in response to certain cues and can be (at least partially) carried through mitosis, thereby constituting an epigenetic ‘memory’ (3).
Chromatin marks do not just turn genes ‘On’ or ‘Off’, at least one other state exists, the ‘Poised’ situation. Poised genes simultaneously harbor activating and repressive histone marks on their promoter, rendering it ‘bivalent’ (4). Bivalency is especially prevalent on key developmental genes, where it may facilitate their rapid induction during developmental transitions when they resolve into active or repressed states. Bivalency exists in different cell types, yet a particularly useful model to study bivalency is mouse ES cells, in which poised promoters were in fact initially discovered (5–7). Work with this system over the past 15 years has led to several key insights, yet the molecular foundations are not well understood. One important question that remains to be elucidated is: which factors maintain bivalent promoters in a poised state?
This paper aims to identify regulators of bivalent genes using minimally biased discovery tools (8). To accomplish this, we built a reporter cell line in mouse ES cells and carried out a primary CRISPR KO screen, followed by a custom secondary screen using Fluidigm screening. Through this approach, we discover that BEND3, a transcription factor previously known to act at ribosomal DNA (9), G4 quadruplexes (10) and pericentric heterochromatin (11,12), binds and regulates hundreds of bivalent promoters, via a consensus site that contains a CpG and is only bound when unmethylated. This confirms a recent report, which arrived at similar conclusions by independent approaches (13). Beyond confirming these data, we describe a new role of BEND3, as we show that Bend3 loss leads to a decrease of DNA hydroxymethylation and a global increase of DNA methylation, which indirectly affect hundreds of non-bivalent genes.
MATERIALS AND METHODS
Cell culture
Mouse ES cells were grown on gelatin-coated dishes in serum/LIF medium containing DMEM/GlutaMAX supplemented with 15% fetal bovine serum (FBS), non-essential amino acids (NEAA), penicillin/streptomycin and 1000 U/ml leukemia inhibitory factor (LIF). For certain experiments, ESCs were adapted to 2i/Vitamin C/LIF medium containing serum-free DMEM-F12 and Neurobasal medium supplemented with 1% N2, 2% B27, 100 μg/ml ascorbic acid, 1 μM of PD0325901 and 3 μM of CHIR99021. The cells were incubated in a humidified atmosphere at 37°C under 5% CO2.
Cloning of sgRNA, transfection and transduction in mES cells
Single guide RNAs (sgRNAs) were designed using Benchling/CRISPOR software. The sgRNAs were cloned into either vector PX459 (Addgene #62988) or vector lentiCRISPRv2 (Addgene #52961). For the transfection of ESCs, we used an Amaxa 4D nucleofector (Lonza). Lentiviral particles were produced by calcium phosphate transfection of HEK293T with psPAX2 and pMD2.G plasmids in a BSL3 tissue culture system. Oligonucleotide sequences are listed in Supplementary Table S1.
Generation of DASH (dazl-mScarlet-HygromycinR) reporter cell line
A P2A-mScarlet-T2A-HygroR reporter cassette (GenScript) was inserted in-frame into exon 6 of the mouse Dazl gene. The cassette was flanked by Dazl homology arms corresponding to the endogenous intron 5-exon 6 and intron 6 sequences, respectively. To prevent re-cutting by Cas9 after cassette insertion by homologous recombination, the Protospacer Adjacent Motif (PAM) sites of the two sgRNAs targeting Dazl exon 6 were mutated within the homology arms. The synthesized cassette was cloned into pUC57-Simple. Two sgRNAs targeting Dazl exon 6 were cloned into the pSpCas9(BB)-2A-GFP backbone (Addgene #48138). Homologous incorporation of the reporter cassette into one of the Dazl alleles was confirmed by PCR and sequencing. Oligonucleotide sequences are listed in Supplementary Table S1.
Uhrf1 KO: validation of DASH mES cells
Two sgRNAs targeting Uhrf1 exon 2 were designed, cloned in the PX459 vector and transfected into DASH cells to test the reporter, followed by a short pulse of 1 μg/ml puromycin for 24 h. Knockout efficiency was assessed by western blotting. Oligonucleotide sequences are listed in Supplementary Table S1.
CRISPR KO screen and generation of clonal bend3 KO lines
The screen and the CRISPR KOs were performed as detailed in Gupta et al. (14). Screen results (top 100 hits) are listed in Supplementary Table S2.
Microfluidic RT-qPCR (fluidigm): cDNA synthesis, RT-qPCR and analysis
Total RNA was extracted and quantified as described above. RNA quality and integrity were analyzed by capillary electrophoresis with the Fragment Analyzer (Agilent Technologies) to calculate the RNA quality number (RQN) for each sample. Defined on a scale ranging from 1 to 10, the mean RQN of the 52 samples was 9.9, indicating very good RNA quality. The cDNA synthesis was performed using Reverse Transcription Master Mix (Fluidigm) according to the manufacturer's protocol in a final volume of 5 μl containing 200 ng total RNA. Reverse transcription was performed using a Nexus Thermocycler (Eppendorf) with the following parameters: 5 min at 25°C, 30 min at 42°C followed by heat-inactivation of reverse transcriptase for 5 min at 85°C and immediately cooled to 4°C. cDNA samples were diluted 5× by adding 20 μl of low TE buffer (10 mM Tris, 0.1 mM EDTA, pH 8.0) (TEKNOVA) and stored at –20°C before specific target pre-amplification. Diluted cDNA was used for multiplex pre-amplification in a total volume of 5 μl containing 1 μl of 5× PreAmp Master Mix (Fluidigm), 1.25 μl of cDNA, 1.25 μl of pooled assay (Exiqon) with an original concentration of each assay of 10 μM and 1.5 μl of nuclease-free water. The cDNA samples were subjected to pre-amplification with the following conditions: 95°C for 2 min, followed by 10 cycles at 95°C for 15 s and 60°C for 4 min. The pre-amplified cDNA was diluted 5× by adding 20 μl of low TE buffer and stored at –20°C before RT-qPCR. The expression of 48 target genes was quantified in 52 samples by quantitative PCR using the high-throughput platform BioMark HD System on two 48 × 48 GE Dynamic Arrays (Fluidigm). Pre-amplified cDNA was used in a total volume of 15 μl containing 7.5 μl of 2× Gene Expression PCR Master Mix (Thermo Fisher Scientific), 3 μl of pre-amplified cDNA, 0.75 μl of 20× EvaGreen (Biotium), 1.5 μl of each primer (Exiqon) and 2.25 μl nuclease-free water. Thermal conditions for qPCR were: 50°C for 2 min before initial heat-activation of DNA polymerase at 95°C for 10 min, followed by 40 cycles at 95°C for 15 s and 60°C for 1 min. Melt curve analysis was performed with an increase of 0.2°C/s, starting from 65°C for 5 s and ending at 95°C. The parameters of the thermocycler were set with ROX as a passive reference and a single probe EvaGreen as a fluorescent detector. To determine the quantification cycle Cq, data were processed by an automatic threshold for each assay, with linear derivative baseline correction using BioMark Real-Time PCR Analysis Software 4.0.1 (Fluidigm). The quality threshold was set at the default setting of 0.65. Second data analysis was performed using R software to measure the relative gene expression using the comparative Cq method with the efficiency corrected method of Pfaffl, after normalization with reference genes (Cq mean of Actinb, Ppia and Rplp0). Fold change between experimental and control groups was calculated for each sample as the difference of Cq between reference genes and the gene of interest (GOI) in control and experimental conditions.
Flow cytometry
The number of cells expressing mScarlet was determined in the ImagoSeine Core Facility (Institut Jacques Monod) using the Influx or FACSAria Fusion cell sorter (BD Biosciences) with yellow laser (561 nm). The percentage of positive cells was determined by flow cytometry using the mScarlet signal intensity threshold (minus background fluorescence from wild-type mESCs). Data were analyzed using FlowJo software (BD Biosciences).
Isolation of genomic DNA
Genomic DNA was isolated from cells by a 200 μg/ml proteinase K treatment at 55°C overnight, followed by 20 μg/ml RNase A treatment at 37°C for 1 h, and extraction using the standard phenol/chloroform/alcohol method. Genomic DNA was suspended in water and quantified with the Qubit dsDNA BR Assay kit on a Qubit 2.0 Fluorometer (Thermo Fisher Scientific) system (Agilent), and only samples with DNA Integrity Number >9 were used for subsequent analysis.
Identification of CRISPR-induced mutation
Genomic DNA was isolated as described above. sgRNA target region amplicons were generated by PCR using 100 ng genomic DNA, 100 nM primers and Platinum Taq polymerase (Thermo Fisher Scientific). The amplicons were then sequenced. Oligonucleotide sequences are listed in Supplementary Table S1.
Rescue experiments, piggyBac system
In the rescue experiments, the Bend3 coding sequence was synthesized (GenScript), with silent mutations incorporated within the PAM to prevent cutting by the sgRNAs that may still be expressed in the KO clone cell lines. This CDS sequence was cloned into a piggyBac vector and co-transfected with PB transposase for stable insertion (15). Empty piggyBac vectors served as controls. Transfected cells were selected with 5 μg/ml of blasticidin for 5 days and processed for phenotypic and molecular assays.
RT-qPCR
Total RNA was extracted from cells using the RNeasy Plus Mini kit (Qiagen) and quantified using the Qubit RNA BR Assay kit on a Qubit 2.0 Fluorometer (Thermo Fisher Scientific). One microgram of total RNA was reverse transcribed using SuperScript IV Reverse Transcriptase (Thermo Fisher Scientific) and Oligo dT primers (Promega). Actinb, Ppia and Rplp0 were used for normalization. Oligonucleotide sequences are listed in Supplementary Table S1.
Western blotting
Cells were harvested and lysed in RIPA buffer with protease inhibitors (Thermo Fisher Scientific), sonicated in a Bioruptor (Diagenode) with 30s ON/30s OFF for 5 min, then centrifuged at 16 000g for 5 min at 4°C. The supernatant was collected and quantified by BCA assay (Thermo Fisher Scientific). Thirty micrograms of protein extract per sample were mixed with NuPage 4× LDS Sample Buffer and 10× Sample Reducing Agent (Thermo Fisher Scientific) and denatured at 95°C for 5 min. Samples were separated in precast SDS-PAGE 4–12% gradient gels (Thermo Fisher Scientific) at 120 V electrophoresis for 90 min and blotted onto a nitrocellulose membrane (Millipore). The membranes were blocked with 5% fat-free milk/PBS for 1 h at RT and incubated overnight at 4°C with the appropriate primary antibody; after washing three times with PBS/0.1% Tween20, the membranes were incubated with homologous fluorescent secondary antibodies and visualized with the LI-COR Odyssey-Fc imaging system. The following antibodies were used in this study: α-BEND3 (Abcam #ab220896, 1:1000), α-DAZL (Abcam #ab34139, 1:500), α-DNMT3A (Cell Signaling #3598, 1:500), α-DNMT3B (Novus #56514, 1:500), α-UHRF1 (SCBT #98817; 1.1000): 1000), α-V5 (Abcam #ab206566, 1:1000), α-TUBULIN (Abcam #7291, 1:1000) and α-GAPDH (Abcam #ab9485, 1:1000). The following secondary antibodies were used in this study: IRDye 800CW Donkey α-Rabbit (Licor #926-32213, 1:5000) and IRDye 680RD Donkey α-Mouse (Licor #926-68072, 1:5000).
DNA methylation analysis: luminometry methylation assay (LUMA)
To assess global CpG methylation, 500 ng of genomic DNA was digested in parallel reactions with MspI + EcoRI and HpaII + EcoRI (NEB), with EcoRI as an internal reference. The rate of CpG methylation is defined as the HpaII/MspI ratio. Samples were analyzed using a PyroMark Q24 Advanced Pyrosequencer (Qiagen).
DNA methylation analysis: LC–MS/MS
Genomic DNA was extracted with RNase A as described above, plus an additional digestion step: 1 μg of DNA was treated with 10U DNA Degradase Plus (ZymoResearch) at 37°C for 4 h, followed by inactivation of the enzyme at 70°C for 20 min and then Amicon Ultra-0.5 ml 10 K centrifugal filters (Merck Millipore) were used to filter the solution. The reaction mixture retained on the centrifuge filter was processed for LC–MS/MS analysis; analysis of total 5-mdC and 5-hmdC concentrations was performed using a Q exactive mass spectrometer (Thermo Fisher Scientific). The instrument was equipped with an electrospray ionization source (H-ESI II Probe) coupled to an Ultimate 3000 RS HPLC (Thermo Fisher Scientific). A ThermoFisher Hypersil Gold aQ chromatography column (100 mm × 2.1 mm, 1.9 μm particle size) heated to 30°C was injected with digested DNA. The flow rate was set to 0.3 ml/min and the column was run for 10 min in isocratic eluent consisting of 1% acetonitrile in water containing 0.1% formic acid. Parent ions were fragmented in positive ion mode, parallel reaction monitoring (PRM) mode at 10% normalized collision energy; MS2 resolution was 17 500, AGC target was 2e5, maximum injection time was 50 ms, and separation window was 1.0 m/z. The inclusion list contained the following masses: dC (228.1), 5-mdC (242.1) and 5-hmdC (258.1). Extracted ion chromatograms (±5 ppm) of basic fragments were used for detection and quantification (dC: 112.0506 Da; 5-mdC: 126.0662 Da; 5-hmdC: 142.0609). Calibration curves were previously generated using synthetic standards in the ranges of 0.2–10 pmol injected for dC and 0.02 to 10 pmol for 5mdC and 5hmdC. Results were expressed as % of total dC.
RNA-sequencing: library preparation
A total amount of 1 μg total RNA per sample was used as input material for the RNA sample preparations. RNA samples were spiked with ERCC RNA Spike-In Mix (Thermo Fisher Scientific). Sequencing libraries were generated using NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB) following the manufacturer's recommendations. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5×). First-strand cDNA was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase (RNase H–). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. In the reaction buffer, dNTPs with dTTP were replaced by dUTP. The remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3′ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure was ligated to prepare for hybridization. To select cDNA fragments of preferentially 250–300 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter). Then 3 μl USER Enzyme (NEB) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95°C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. Lastly, products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system.
RNA-sequencing: read alignment
FASTQ reads were trimmed using Trimmomatic (v0.39) and parameters: ILLUMINACLIP:adapters.fa:2:30:10 SLIDINGWINDOW:4:20 MINLEN:36. Read pairs that survived trimming were aligned to the mouse reference genome (build mm10) using STAR (v2.7.5c) and default single-pass parameters. PCR duplicate read alignments were flagged using Picard-tools (2019) MarkDuplicates (v2.23.4). Uniquely aligned, non-PCR-duplicate reads were kept for downstream analysis using Samtools view (v1.10) and parameters: -q 255 -F 1540. Gene expression values were calculated over the mm10 NCBI RefSeq Genes annotation using VisRseq (v0.9.12) and normalized per million aligned reads per transcript length in kilobases (RPKM). Bigwig files were generated using deeptools bamCoverage (v3.3.0) using counts per million (CPM) normalization and visualized in IGV (v2.8.9).
RNA-seq: differential expression and PCA plots
DESeq2 (v1.30.0) was employed using apeglm LFC shrinkage to calculate differential expression (Supplementary Table S3). Genes or transposable elements were categorized as significantly differentially expressed if they showed an expression fold-change ≥2 and associated adjusted P-value <0.01. PCA plots were generated using variance-stabilized RNA-seq read counts using ‘varianceStabilizingTransformation’, and plotted using the DESeq2 function ‘plotPCA’ and visualized using the R package ‘ggplot2’ (v3.3.3) in R (v4.0).
RNA-seq: transposable element quantification
RepeatMasker (last updated 2012-02-06) was downloaded from the UCSC Table Browser. To measure the expression of transposable element families, PCR duplicates were removed and all reads, including uniquely mapped and multi-mapped reads, were enumerated using VisRseq. Multi-mapped reads were counted once, and all individual elements were aggregated to calculate family-wide expression in read count (for differential expression analysis) and RPKM values.
Chromatin immunoprecipitation (ChIP-seq)
Ten million cells were cross-linked (1% formaldehyde, 5 min RT); the cross-linking stopped by adding glycine at 125 mM final). After a PBS wash the cells were resuspended in swelling buffer (0.5% NP-40, 0.85 mM KCl, 1 mM PMSF, 5 mM PIPES pH 8.0; protease inhibitors) and incubated on ice for 20 min. After centrifugation, the cell nuclei were resuspended in IP Buffer (0.1% SDS, 1% Triton X-100, 3 mM EDTA, 1 mM PMSF, 150 mM NaCl, 25 mM Tris–HCl pH 8.0; with protease inhibitors) and sonicated in a Bioruptor Pico (Diagenode) for 5 min (30 s ON/30 s OFF continuously) to generate 200–500 bp fragments. The fragmented chromatin (50 μg) was immunoprecipitated overnight at 4°C in IP buffer with 1 μg of antibody. The following antibodies were used in this study: α-BEND3 (Abcam #ab220896), α-V5 (Abcam #ab206566), IgG (Abcam #ab46540). Subsequent steps were performed using the Pierce Magnetic ChIP Kit (Thermo Fisher Scientific). The resulting DNA was purified using the ChIP DNA Clean&Concentrator Kit (Zymo Research), and libraries were prepared using the KAPA HyperPrep Kit (Roche) according to the manufacturer's instructions.
ChIP-seq analysis
FASTQ reads were trimmed with Trimomatic (v0.39), parameters ILLUMINACLIP:ILLUMINA_adapters.fa:2:30:10 SLIDINGWINDOW:4:20 MINLEN:36 Alignment of trimmed reads was performed in –local mode in Bowtie2 (v2.4.1). After alignment, duplicate-marked bam files were created using Picard (v2.23.4) CleanSam, SamFormatConverter, SortSam and Markduplicates. The resulting bam files were then converted to bigwig with deptools (v3.3.0) Bamcoverage and the option –ignoreDuplicates -normalizeUsing CPM -minMappingQuality 10 -ignoreForNormalization chrX chrY chrM. Peaks were called using MACS2 with default settings and motif enrichment analysis was carried out in HOMER. For datasets that were already publicly available, SRA files were downloaded from NCBI GEO using the Eearch and efetch utilities (v15.9) and converted to FASTQ using SRAtoolkit (v2.8.0). Published datasets were reanalyzed from GSE98149 (16).
Whole-genome-bisulfite sequencing (WGBS)
Genomic DNA was extracted as described for LC–MS/MS; library preparation for WGBS was performed using the previously described tPBAT protocol (17,18). For library preparation, 100 ng of genomic DNA spiked with 1% (w/w) unmethylated lambda DNA (Promega) was used. Sequencing was performed by MacroGen Japan, Inc. using a HiSeq X Ten system; 8 lanes were allocated for analysis of 20 samples. Sequenced reads were mapped in BMap and compiled in an in-house pipeline as previously described, and custom scripts were archived using GitHub (https://github.com/FumihitoMiura/Project-2). Basic metrics for the methylome data are described in Supplementary Table S4. DNA methylation levels on CpGs covered by at least five reads were averaged over the target region of 2 kb bins genome-wide, over CpG island (cpgIslandExt, n = 16 023, last updated 2012-02-09), enhancer elements (Ensembl Regulatory features release 81, n = 73 796), promoters (NCBI Refseq, TSS ± 1 kb, n = 24 371), gene bodies and transposable elements (RepeatMaster, n = 5 147 736).
RESULTS
Principle of the genetic screen
Epigenetic silencing is dynamic in mouse ES cells and responds to culture conditions. In particular, a number of genes are expressed in 2i condition, but are repressed when the cells are switched to serum (Figure 1A). The principle of our screen was to select one such gene and replace it with selectable markers to carry out a CRISPR KO screen in serum, and to positively select the cells in which the KO caused a loss of silencing (Figure 1A). To select the gene in question, we considered three criteria: mechanism(s) of epigenetic repression, biological function(s) of the gene and expression level in repressed and induced condition (to ensure high enough expression of the selection markers). These criteria led us to select Dazl as our target of choice. Dazl encodes an RNA-binding protein which is critical during gametogenesis (19), but it is also expressed during other developmental events, as in the 2-cell-embryo (20). The promoter of Dazl overlaps a large CpG island (1.1 kb); in serum-grown ES cells, the 5′ half of the CpG island is heavily methylated on CpGs, while the 3′ half is mostly void of DNA methylation. In addition, the promoter of Dazl contains both the activating H3K4me3 and the repressive H3K27me3 marks and can therefore be classified as a bivalent promoter (Figures 1B and S1A). Dazl expression is low in serum, however genetic experiments have shown that this repression can be lifted either by impeding DNA methylation (21,22), or by interfering with the ncPRC1.6 complex (23,24). One of our goals was to identify new actors in this process (Figure 1C).
Having chosen Dazl as our reporter gene of interest, we then used CRISPR/Cas9-guided homologous recombination to insert two selectable markers in the gene: mScarlet, encoding a bright red fluorescent protein (25), and HygroR, a bacterial enzyme that renders cells resistant to Hygromycin. The two markers were inserted in-frame in the Dazl exon 6 and were separated by self-cleaving peptides (Supplementary Figure S1B). Recombinant clones were picked and characterized by genomic PCR and sequencing. One heterozygous integrant was further characterized; it contains one wild-type Dazl allele, and the mScarlet-HygroR insertion in the other allele (Figure 1D). We will refer to this line as the DASH (Dazl-Scarlet-Hygro) reporter line.
We next carried out a proof-of-concept experiment to verify that the selectable markers could be induced in serum upon the KO of an epigenetic repressor. For this, we used CRISPR to knock out a key factor in DNA methylation maintenance, UHRF1 (Supplementary Figures S1C and S1E). As expected (26,27), the Uhrf1 KO caused a global loss of DNA methylation as seen by LUMA (Figure 1F). The DASH cells grown in serum became mScarlet-positive and Hygromycin-resistant upon Uhrf1 KO (Figures 1G, H, and S1D), confirming the feasibility of a genetic screen.
Primary screen
We grew the DASH cells in serum and infected them with a lentiviral library of ∼80 000 vectors co-expressing Cas9 and sgRNAs targeting ∼20 000 genes (28). The multiplicity of infection was ∼0.1, the coverage ∼150 infected cells per guide and two independent screens were carried out in parallel. After an initial Puromycin selection to eliminate non-infected cells, Hygromycin selection was applied and the fluorescent cells were purified by FACS (Figure 2A). The sgRNAs present in these cells were amplified by PCR and sequenced by NGS, and the sequencing data statistically analyzed using MAGeCK (29). This procedure yielded an ordered list of candidates, ranked by P-value, and a highly stringent cutoff (P-value lower than 5 × 10−6 and FDR lower than 0.01) yielded a list of 12 hits (Figure 2B). Four of the top 12 hits (ranked #2, 4, 5 and 11 and shown in light blue in Figure 2B) correspond to factors already known to repress Dazl: UHRF1, E2F6, MGA and DNMT1. Their presence on the list was expected and supports the validity of the screen. Five of the hits (in light green in Figure 2B) correspond to genes involved in pluripotency or in Heparan Sulfate synthesis.
Secondary screen on the top hits and focus on BEND3
We next sought to validate the results of the primary screen while also prioritizing the hits for further analysis. For this, we carried out a secondary screen based on gene expression analysis.
We first generated, by lentiviral infection with CRISPR/Cas9 vectors, new KO populations for ten of the 12 top hits (and 3 different negative control populations), as depicted in Supplementary Figure S2A. For technical reasons, Grb2 was excluded from the analysis and, instead of knocking out B4galt7, we targeted B3galt6, a gene in the same complex and with the same function. Rather than carrying out a whole-genome RNA-seq on each hit, we examined the expression of 21 genes of interest (along with 3 normalizers) in the aforementioned knockouts using the Fluidigm technology.
The genes are color coded in Supplementary Figure S2B and belong to different functional classes: besides the mScarlet reporter and Dazl itself, there were other germline genes (Asz1, Hormad1, Fkbp6 and 1700013h16Rik); markers of naive ES cells (Nanog, Prdm14, Zfp42/Rex1); repeated elements (MERVL, MMERVK10C, IAPEz, IAP, MLV, MusD, LINE1 and GLN); and finally epigenetic regulators (Uhrf1, Dnmt1, Dnmt3a, Dnmt3b, Dnmt3l, Tet1 and Tet2).
Unsupervised analysis of the Fluidigm data revealed that the KOs were segregated into 4 clusters. Cluster number 1 contained both the Dnmt1 and Uhrf1 KOs, which expectedly had very similar expression profiles. In accordance with published data (30), these two KOs reactivated the expression of all germline genes, as well as certain repeats, in particular GLN. The Ptpn11 KO was the third and final member of this cluster.
Cluster number 2 included Slc35b2 and B3galt6 KOs, along with Mga and E2f6, which behaved similarly, also as predicted. Unsurprisingly (23,24), all germline genes were reactivated in this cluster.
The third cluster contained the Bend3, Spop and Tcf7l1 KOs. Contrary to cluster #2, the germline genes were either unaffected or only slightly upregulated in this cluster. The fourth cluster contained all three negative controls, along with the Emc1 KO. Dazl was not reactivated in this cluster, and the other genes were similarly unaffected.
This analysis leads to several useful conclusions. First, 10 out of 11 new KO populations (91%) displayed the induction of Dazl that was also observed in the primary screen, therefore validating the hits. The only false positive was Emc1. Secondly, hits affecting DNA methylation or ncPRC1.6 form separate clusters. Third, three hits generate a gene expression pattern that resembles neither DNA methylation mutants nor the ncPRC1.6 mutants. These potentially novel hits are Bend3, Spop and Tcf7l1. Of those three hits, Spop encodes an adaptor protein for an E3 ubiquitin ligase complex (31,32) and as it seemed least likely to act directly at the Dazl promoter we did not pursue it further in this study. Tcf7l1 is a transcription factor that inhibits the naive state in mouse ES cells (33); as a consequence, cells lacking TCF7L1 tend to remain in a 2i-like state, even when grown in serum (34). Because cells in 2i express Dazl, it seemed possible that the upregulated expression of Dazl in the Tcf7l1 KO cells was merely an indirect reflection of a change of cellular identity towards a more naive state. The third hit within the cluster was Bend3, which encodes a chromatin protein (11) known to repress transcription at the rDNA (9). As BEND3 had suitable characteristics for being a new direct repressor of Dazl, we decided to investigate it further.
Isolation of bend3 mutant clones and rescue experiments
To validate the results of the CRISPR screen, we generated new Bend3 mutant populations (Figure 3A). For this, we infected DASH cells with a mix of two distinct lentiviruses, each expressing Cas9 and a specific sgRNA targeting Bend3 (Supplementary Figure S3A). As Bend3 produces two protein isoforms, we targeted the last exon of Bend3, which is included in both isoforms.
Two independent Bend3-mutant populations showed a roughly 10-fold increase of mScarlet-expressing cells relative to control, similar to the magnitude observed for Uhrf1 mutant cells (Supplementary Figure S3B). To carry out more precise experiments, we isolated individual Bend3 mutant clones from these populations (see scheme in Figure 3A). Three clones were isolated, each of them containing a high proportion of mScarlet-expressing cells (25–35%, relative to ∼3% in control cells, Figures 3B and S3C). These clones also expressed more Dazl mRNA (Figure 3C) and DAZL protein (Figure 3D) from the non-targeted allele. Upon sequencing their genomic DNA, we found that two clones had non-sense or frameshift mutations at AA 308, removing about two thirds of the protein. The last clone had a small homozygous in-frame mutation removing four amino acids in the first BEN domain (Supplementary Figure S3D). The BEND3 protein was not detectable in the KO clones (Supplementary Figure S3E), suggesting that the mutant forms are unstable.
We next verified that the phenotypes could be rescued by reintroducing a WT cDNA (see scheme in Figure 3A); for this, we used a PiggyBac transposon-based system (30) driving constitutive expression of V5-tagged BEND3. Expression of the functional Bend3 cDNA in the Bend3 KO cells decreased the number of mScarlet-positive cells to near-background levels (Figures 3E and S3F), rendered the cells sensitive to Hygromycin (Figure 3F), and silenced the expression of the DAZL protein (Figure 3G).
Together these experiments support the conclusion that BEND3 is a bona fide repressor of Dazl expression.
Transcriptome analysis in bend3 mutant cells reveals up- and down-regulated genes
To examine more globally the effects of Bend3 deletion, we performed RNA-seq on the 3 independent mutant clones, grown in serum, which we compared to WT cells (grown in serum or 2i). As shown in Figure 4A, the mutant cells clustered away from the WT cells, and they were highly distinct from 2i grown cells. This supports the idea that the Bend3 KOs do not reactivate Dazl because they enter a 2i-like state.
More genes were repressed than induced in the Bend3 KO cells (635 versus 325 respectively, Figure 4B and S4A). As expected, Dazl was among the induced genes (Figure 4B). A recent publication independently examined the transcriptome of Bend3-/- ES cells (31); we compared our respective RNA-seq results and found that the overlap was significant yet small, with only ∼10–20% upregulated or downregulated genes in common (Supplementary Figure S4B). We speculate this dissimilarity results from different procedures for the mutant cell derivation, combined with different cell culture conditions.
We also examined the behavior of transposable elements in the RNA-seq data (Supplementary Figure S4C). The mutation of Bend3 (left panel) had very little effect on the expression of mouse endogenous retroviruses (ERVs). As a control (right panel), we compared WT cells grown in 2i and serum and predictably found a large induction of ERVs such as IAPs, which validates our analysis pipelines.
Returning to protein-coding genes, we manually inspected some of the differentially expressed (DE) genes and noted that many of them had bivalent marks on their promoter, in line with the recent results of Zhang et al. (13). To determine if this observation could be generalized, we used a list of very stringently curated bivalent genes (35). We cross-referenced this list with our RNA-seq data and found that ∼35% of differentially genes in the Bend3 KO relative to WT are bivalent; this is a highly significant enrichment (P = 8 × 10−54, Figure 4C). The enrichment is stronger for genes that are upregulated in the Bend3 mutants (50% of those are bivalent) than for genes that are downregulated in the Bend3 mutant (27% being bivalent, Figure 4C). Correspondingly, genes induced in the Bend3 mutant were more likely to have a bivalent than a repressed status in the WT cells (Figure 4D). Bivalent chromatin marks often mark developmental genes in mouse ES cells, and in line with this we found that the genes differentially expressed in Bend3 KO cells were markedly enriched for developmental regulators (Supplementary Figure S4D).
Lastly, we performed RT-qPCR in Bend3 KO cells that had been rescued with the V5-BEND3 construct (Figure 4E). Genes that had been downregulated in the KO cells (Cpt1a, Tead4, Klf15, Col4a1, Tbx3) were upregulated in the rescue cells; conversely, genes that had been induced in the KO cells (Dazl, Bend3 (measured with primers that do not amplify the rescue construct), Pax6, Ccdc106) were repressed in the rescue cells (Figure 4E). All of the genes examined in this experiment are directly bound by BEND3, as determined by the ChIP-seq that will be presented later (Figure 5).
Altogether, these data show that removing BEND3 from serum-grown ES cells alters the expression of hundreds of genes, with a high proportion of bivalent genes. Nevertheless, of the ∼3600 bivalent genes, only about 10% are regulated by BEND3 (Figure 4C). To understand this situation, we next sought to identify the direct targets of BEND3.
ChIP-sequencing identifies the direct targets of BEND3
We then used ChIP-seq to identify genes directly regulated by BEND3. Using an antibody directed against endogenous BEND3, we identified ∼7400 binding peaks in the ES cell genome (Figure 5A). This dataset is highly similar to the dataset of Zhang et al. (13), who used the same antibody as we did (Supplementary Figure S5A). We also performed ChIP-seq with an antibody against the V5 tag, in Bend3 KO cells expressing V5-BEND3. This data obtained overlapped extremely well with the ChIP-seq data obtained on the endogenous protein (Supplementary Figure S5B-C). These two lines of evidence strongly support the quality of our ChIP-seq dataset.
Of the ∼7400 peaks detected in the BEND3 ChIP, 47% occur in promoters, 6% in enhancers, 19% in gene bodies and 28% in other intergenic regions (Figure 5A). A de novo search for motifs enriched in the peaks revealed a highly enriched motif with a strong selection for the core sequence CCCACG (Figure 5B). This motif occurred in 70% of all BEND3 peaks, often in multiple instances (Supplementary Figure S5D), hinting that BEND3 binds it directly, which is supported by recent in vitro work (13,36). The overwhelming majority (94%) of promoters bound by BEND3 contain a CpG island, which is a significant enrichment (Figure 5C). BEND3 binding is prevalent at bivalent promoters; in contrast it is extremely rare at promoters bearing H3K9me3 (Supplementary Figure S5E).
We next intersected our RNA-seq data with our BEND3 binding data (Figure 5D). Only a minority of promoters (5.5%) bound by BEND3 become significantly activated or repressed upon KO. Conversely, only about 23% of genes deregulated upon Bend3 deletion have BEND3 binding at their promoter. The genes bound by BEND3 and differentially expressed in the KO, which are likely to be direct targets, are split equally between upregulated genes and downregulated genes (Figure 5E). Such genes include Dazl, which is repressed by BEND3 (Figure 5F), but also Col4a1 and Col4a2, which are activated by BEND3 (Figure 5G). Again, most direct targets of BEND3 (65%) have a bivalent promoter in WT cells (Figure 5E, Supplementary Figure S5E-F); they have a variable number of BEND3 consensus motifs, with 10% having no motif, to 20% having four motifs or more (Figure S5G).
Global DNA hypermethylation and reduced 5-hmC in bend3 mutant cells
The results above establish that many of the genes that are up- or down-regulated when BEND3 is absent are not direct targets of the protein. To investigate what might underpin these indirect effects of BEND3, we examined DNA methylation.
We first measured global DNA methylation with a method based on restriction enzymes, LUMA (37). The WT cells grown in serum and 2i media served as controls and, as expected, global DNA methylation was low in 2i and high in serum (Figure 6A). Additionally, the level of DNA methylation was higher in Bend3 mutant cells than in WT cells (both grown in serum). Importantly, rescuing the loss-of-function mutations with a Bend3-expressing construct brought the level of DNA methylation back to WT levels (Figure 6A).
To validate these observations with an orthogonal approach, we turned to LC–MS/MS: this experiment confirmed that 5mC was more abundant in the Bend3 KO cells (Figure 6B).
We then performed Whole-Genome Bisulfite Sequencing (WGBS) on WT cells, and on each individual Bend3 mutant clone (Supplementary Figure S6A). As anticipated, the 2i-grown WT cells were vastly undermethylated relative to serum-grown WT cells (Figure 6C), and the methylation of genomic tiles in serum-grown WT cells displayed a typical bimodal profile: a minority of tiles unmethylated and most tiles methylated, with a peak around 80% methylation (Figure 6C), lending confidence to our WGBS experiment and analysis. Upon comparison to WT cells, we saw that the Bend3 mutants had a methylation distribution that was still bimodal but in which the methylated compartment was shifted towards higher methylation levels (Figures 6C and S6B). As the three clones behaved similarly, we aggregated their WGBS reads in subsequent analyses.
To carry out a more precise statistical analysis, we narrowed down our calculations to 250-bp tiles containing 5 or more CpGs, for a total of 2.1 million tiles. Of these, 91.3% do not significantly gain or lose methylation upon Bend3 mutation, 8.2% gain DNA methylation and 0.5% lose DNA methylation (Figures 6D, S6C, D).
We then investigated the genomic distribution of the DNA methylation gains. The increase was global and affected promoters, enhancers, gene bodies and transposable elements (Figure 6E). Over the whole genome, 8.2% of the tiles gain DNA methylation, as stated above, however the level of increase differed between genomic elements: enhancers had a marked tendency to gain more DNA methylation, while promoters were slightly under average (Figure 6F). Looking at CpG island, we saw that the CpG islands themselves were rarely hypermethylated, the CpG shores gained slightly less methylation, whereas the CpG shelves gained more methylation than average (Figure 6G). In other words, the gains of DNA methylation occur at a distance from the CpG islands, which is exemplified by the gene Tead4 (Supplementary Figure S6E).
We then crossed our WGBS data with ChIP-seq data reflecting the chromatin state in WT cells (Supplementary Figure S6F). We find that regions that were bivalent or marked by H3K4me3 gain more methylation than average. Regions marked by H3K27Ac, which is typical of enhancers, are also much more likely to gain DNA methylation than the genome average, again suggesting that enhancers are targets for DNA hypermethylation in Bend3 KO cells.
Next, we intersected our WGBS and RNA-seq data. We find that half of transcriptionally downregulated genes in Bend3 KO cells gain DNA methylation in their promoter (Figure 6H). The BEND3-bound peaks themselves, or the genes upregulated in Bend3 KO, gain methylation less than the rest of the genome (Supplementary Figure S6G).
We investigated mechanisms that could possibly underpin the elevated DNA methylation seen in Bend3 mutant cells. Neither our RNA-seq data (Supplementary Figure S6H), nor our western blotting experiments (Supplementary Figure S6I) suggested that BEND3 removal upregulated Uhrf1, Dnmt1, Dnmt3a, or Dnmt3b, which are key actors of DNA methylation. Therefore we speculated that DNA hydroxymethylation might be modified in the Bend3 KOs. To test this idea, we again used LC–MS/MS, and we indeed saw a significant decrease of 5-hmC levels in Bend3 mutant cells (Figure 6I). This effect is not due to transcriptional downregulation of the TET enzymes in the absence of BEND3 (Supplementary Figure S6H), so we considered additional mechanisms that could potentially lead to lower TET activity in Bend3 mutant cells. The TET enzymes are critically dependent on alpha-ketoglutarate, so we tested whether the transcriptional defects occurring in Bend3 KO cells might impinge on the pathways producing this metabolite. A GSEA analysis of our RNA-seq data showed that there was indeed a significant downregulation of the group of genes involved in alpha-ketoglutarate production (Figure 6J). To investigate this finding in more depth, we examined the expression of individual genes belonging to this GO term group, and intersected the data with our BEND3 ChIP-seq results (Figure 6K). We find that several genes of the group are bound by BEND3 and downregulated in its absence. In particular, Idh1 and Idh2, encoding Isocitrate Dehydrogenase, which produces alpha-ketoglutarate, are both directly activated by BEND3. We hypothesize that the lower Idh1/Idh2 expression causes lower alpha-ketoglutarate availability, reducing TET activity, and increasing DNA methylation in cells lacking BEND3.
Integrating the local and global effects of BEND3 on chromatin and transcription
Lastly, we set out to integrate the changes occurring at the chromatin (DNA methylation and histone marks) and transcriptional levels in the Bend3 mutant cells.
For this, we mapped the key histone modification H3K27me3 by CUT&RUN in ES cells with or without BEND3 (Figure 7A). Bivalent promoters were extremely likely to experience changes in H3K27me3 abundance upon BEND3 removal (Figure 7A, left panel), reinforcing the idea that BEND3 has a key role in their homeostasis. Of the promoters directly bound by BEND3, 46% experienced changes in H3K27me3 occupancy, again a significant enrichment (Figure 7A, middle panel). Finally, of the 960 genes differentially expressed upon BEND3 loss, more than half experienced changes in H3K27me3 presence on their promoter (Figure 7A, right panel). The correlations between BEND3 binding, H3K27me3 presence, changes in gene expression and changes in DNA methylation, are explored statistically in Supplementary Figure S7A, and the examples of Dazl and Col4a1/a2 are shown in Supplementary Figure S7B-C.
We also examined some of the interactions more precisely. For instance, we listed the promoters that lose H3K27me3 when BEND3 is absent (Figure 7B). Fifty-four percent of the upregulated genes were in that category, versus only 25% of the downregulated genes. Therefore, transcriptional induction upon Bend3 mutation often involves H3K27me3 loss, in line with many of the target genes having bivalent promoters. The opposite was true for promoters that gain H3K27me3 upon BEND3 loss: they contain more downregulated than upregulated genes (Figure 7C).
Lastly, we crossed our CUT&RUN and WGBS data (Figure 7D). The regions that were occupied by H3K27me3 in WT cells were significantly more likely to gain DNA methylation than the genome average. Similarly, the regions that have lost H327me3 in the Bend3 KO cells were significantly more hypermethylated (Figure 7D). Further examples are given in Supplementary Figure S7D-E.
Altogether these data lead us to propose an integrated model for the local and global effects on BEND3 on chromatin and gene expression (Figure 7E).
DISCUSSION
In this work, we set out to identify new regulators of Dazl using a CRISPR screen in mouse ES cells. Our primary screen identified a shortlist of 12 candidates, of which 4 were known regulators of Dazl. A secondary screen then led us to focus on one particular hit, BEND3. We combined ES cell genetics with various genomic approaches to show that BEND3 regulates not only Dazl, but hundreds of bivalent promoters.
A recent article by Rui-Ming Xu, Bing Zhu and their coworkers, also identified BEND3 as a key regulator of bivalent genes (13). Zhang et al. reached these conclusions by a wholly different approach, which was the generation of Bend3-/- mice and the derivation of ES cells from these mice. We note that, although our conclusions are parallel, some of the specifics are different: for instance, only about 10% of the differentially expressed genes in Bend3 mutant cells relative to WT are the same between our study and Zhang et al. Technical reasons (including genetic background and culture conditions) likely play a role, but a more fundamental difference also exists: Zhang et al. derived their ES cells from Bend3-/- embryos, therefore BEND3 was never present in the ES cells. In contrast, we deleted Bend3 in ES cells that previously expressed the gene. Therefore, the results of Zhang et al. could not discriminate whether BEND3 is involved in bivalent domain establishment, maintenance, or both. Our results establish unequivocally that BEND3 is required for bivalent genes to maintain their poised state, yet they do not rule out an additional role in establishing these domains.
Our ChIP-seq data show that a major determinant of BEND3 distribution on genomic DNA is the presence of consensus sites, with the sequence CCCACGCG. About 70% of all BEND3 peaks contain one or more such sites, and 20% of all BEND3 peaks contain three or more. Our findings are fully consistent with two recent papers showing that the BEN domain 4 of BEND3 recognizes the sequence CCACGC (13), but they differ from another study in which BEND3 was proposed to bind G-quadruplexes (10). The fact that only 10% of all bivalent promoters are regulated by BEND3 is likely directly linked to the presence or absence of the BEND3 consensus sites, but other contributing factors cannot be ruled out. The consensus BEND3 binding motif contains two CpG sites that can potentially be subjected to cytosine DNA methylation. Our data show that bound motifs are unmethylated, while methylated motifs are not bound, similarly to the binding behavior of BANP, a protein related to BEND3 (38). Structurally, CpG methylation has been shown to impair BEND3 binding by blocking an arginine-to-cytosine contact (13,36), which is consistent with our ChIP-seq data.
DNA methylation, therefore, influences BEND3 binding and is conversely regulated by BEND3. Indeed, a striking finding we report is that mouse ES cells lacking BEND3 have decreased DNA hydroxymethylation, combined with increased DNA methylation, suggesting that TET activity is lower in the absence of BEND3. In agreement with this possibility, while the DNA hypermethylation of Bend3 KO cells is global across the genome, it is especially marked at certain regions such as enhancers and CpG shelves, which are normally protected against de novo methylation by the TET enzymes (39–44). We have ruled out a direct transcriptional downregulation of TET1, TET2, or TET3 as the cause for decreased TET activity. However, the TET enzymes are critically dependent on alpha-ketoglutarate (alpha-KG) for catalysis, and the promoters of the two key genes for alpha-KG production, Idh1 and Idh2, are bound and activated by BEND3, therefore we speculate that lower levels of alpha-KG contribute to decreased TET activity in the mutant cells.
Many genes are differentially expressed upon BEND3 removal, even though they are not directly bound. It is likely that some of these indirect targets become repressed because their promoter or enhancer gains abnormal levels of DNA methylation upon BEND3 removal. We and others (13) find that loss of BEND3 leads to a decrease of H3K27me3 at BEND3 sites. This could be explained by a direct effect of BEND3 on H3K27 methylation, such as the recruitment of NuRD, which deacetylates H3K27 as a necessary step towards methylation (10). However, CpG methylation is antagonistic with H3K27 methylation, and our data suggest that increased cytosine methylation in the Bend3 mutant cells could also contribute to the decreased H3K27me3 deposition. Future experiments will help address this possibility.
Supplementary Material
ACKNOWLEDGEMENTS
We are very grateful to the following colleagues for useful advice: Allison Bardin, Sebastian Bultmann, Claire Francastel, Claire Rougeulle, Sophie Polo, Pablo Navarro, Raphael Margueron. We thank the Vectorology platform, Epigenetics platform, Microscopy platform and Bioinformatics/Biostatistics Core Facility (BIBS) at the CNRS Epigenetics and Cell Fate Unit (Université Paris Cité), for providing access and technical advice. We are grateful to Brianna Rodgers for editing the manuscript. We acknowledge the ImagoSeine core facility of the Institut Jacques Monod, member of the France BioImaging (ANR-10-INBS-04) and the support of La ligue contre le Cancer (R03/75–79). Microfluidic RT-qPCR (Fluidigm) analysis was carried out on the qPCR-HD-Genomic Paris Centre Core Facility and was supported by grants from Région Ile-de-France, DIMBIO-RVT-INSERM-ADR-P11 n° 21016711. PAD is supported by Agence Nationale de la Recherche (PRCI INTEGER ANR-19-CE12-0030-01), LabEx ‘Who Am I?’ (ANR-11-LABX-0071), Université de Paris IdEx (ANR-18-IDEX-0001) funded by the French Government through its ‘Investments for the Future’ program, Fondation pour la Recherche Médicale, Fondation ARC (Programme Labellisé PGA1/RF20180206807). PAD, AS and GC were supported by grant RETROMET, ANR-16-CE12-0020, from Agence Nationale de la Recherche. JRA and MVCG were supported by Laboratoire d’excellence Who Am I? (Labex 11-LABX-0071) Emerging Teams Grant and by the European Research Council (ERC-StG-2019 DyNAmecs). This research was supported by Platform Project for Supporting Drug Discovery and Life Science Research (Basis for Supporting Innovative Drug Discovery and Life Science Research (BINDS)) from AMED under Grant Number JP20am0101103 (support number 2652).
Author contributions: N.G. and P.A.D. conceived the project. L.Y., N.G., P.T. and P.A.D. planned the experiments. L.Y., A.A., N.G. and S.B. performed experiments and analyzed the data. F.M. performed WGBS. L.F. performed LUMA, ChIP-seq and CUT&RUN. F.B. performed mass spectrometry. J.R.A. performed RNA-seq and ChIP-seq analysis. E.B. and O.K. performed bioinformatic analyses on WGBS and other data. A.S., K.Y. and G.C. performed mutation analysis. L.Y. and P.A.D. wrote the manuscript. P.A.D., T.I., L.Y. and N.G. supervised the project. P.B.A., G.C., M.V.C.G., T.I. and P.A.D. acquired funding. All authors reviewed the manuscript.
Notes
Present address: Nikhil Gupta, Joint AstraZeneca-Cancer Research Horizons Functional Genomics Centre, Jeffrey Cheah Biomedical Centre, Cambridge, UK.
Contributor Information
Lounis Yakhou, Université Paris Cité, CNRS, Epigenetics and Cell Fate, F-75013 Paris, France.
Anaelle Azogui, Université Paris Cité, CNRS, Epigenetics and Cell Fate, F-75013 Paris, France.
Nikhil Gupta, Université Paris Cité, CNRS, Epigenetics and Cell Fate, F-75013 Paris, France.
Julien Richard Albert, Université Paris Cité, CNRS, Institut Jacques Monod, F-75013 Paris, France.
Fumihito Miura, Department of Biochemistry, Kyushu University Graduate School of Medical Sciences, Fukuoka, Fukuoka 812-8582, Japan.
Laure Ferry, Université Paris Cité, CNRS, Epigenetics and Cell Fate, F-75013 Paris, France.
Kosuke Yamaguchi, Université Paris Cité, CNRS, Epigenetics and Cell Fate, F-75013 Paris, France.
Sarah Battault, Université Paris Cité, CNRS, Epigenetics and Cell Fate, F-75013 Paris, France.
Pierre Therizols, Université Paris Cité, CNRS, Epigenetics and Cell Fate, F-75013 Paris, France.
Frédéric Bonhomme, Institut Pasteur, Université Paris Cité, CNRS, Epigenetic Chemical Biology, UMR 3523, F-75724 Paris, France.
Elouan Bethuel, Université Paris Cité, CNRS, Epigenetics and Cell Fate, F-75013 Paris, France.
Arpita Sarkar, Université Côte d’Azur, Inserm, CNRS, IRCAN, Nice, France.
Maxim V C Greenberg, Université Paris Cité, CNRS, Institut Jacques Monod, F-75013 Paris, France.
Paola B Arimondo, Institut Pasteur, Université Paris Cité, CNRS, Epigenetic Chemical Biology, UMR 3523, F-75724 Paris, France.
Gael Cristofari, Université Côte d’Azur, Inserm, CNRS, IRCAN, Nice, France.
Olivier Kirsh, Université Paris Cité, CNRS, Epigenetics and Cell Fate, F-75013 Paris, France.
Takashi Ito, Department of Biochemistry, Kyushu University Graduate School of Medical Sciences, Fukuoka, Fukuoka 812-8582, Japan.
Pierre-Antoine Defossez, Université Paris Cité, CNRS, Epigenetics and Cell Fate, F-75013 Paris, France.
DATA AVAILABILITY
The accession number for the deep sequencing data reported in this paper are: GSE221710 (CRISPR screen, Mutational analysis) and GSE224077 (RNA-seq, WGBS, ChIP-seq, CUT&RUN-seq). The scripts used are available on Github: https://github.com/FumihitoMiura/Project-2 (permanent DOI: 10.6084/m9.figshare.23899266) and https://github.com/kirsho/DASH (permanent DOI: 10.6084/m9.figshare.23899254).
SUPPLEMENTARY DATA
Supplementary Data are available at NAR Online.
FUNDING
LabEx “Who Am I?” [ANR-11-LABX-0071]; Agence Nationale de la Recherche [INTEGER ANR-19-CE12-0030-01, REMEDY ANR-21-CE12-0015-03, RETROMET ANR-16-CE12-0020]; Université de Paris [ANR-18-IDEX-0001]; Platform Project for Supporting Drug Discovery and Life Science Research; European Research Council [ERC-StG-2019]; Fondation ARC [PGA1/RF20180206807]. The open access publication charge for this paper has been waived by Oxford University Press – NAR Editorial Board members are entitled to one free paper per year in recognition of their work on behalf of the journal.
Conflict of interest statement. None declared.
REFERENCES
- 1. Peter I.S., Davidson E.H.. Implications of developmental gene regulatory networks inside and outside developmental biology. Curr. Top. Dev. Biol. 2016; 117:237–251. [DOI] [PubMed] [Google Scholar]
- 2. Janssen S.M., Lorincz M.C.. Interplay between chromatin marks in development and disease. Nat. Rev. Genet. 2022; 23:137–153. [DOI] [PubMed] [Google Scholar]
- 3. Cavalli G., Heard E.. Advances in epigenetics link genetics to the environment and disease. Nature. 2019; 571:489–499. [DOI] [PubMed] [Google Scholar]
- 4. Blanco E., González-Ramírez M., Alcaine-Colet A., Aranda S., Di Croce L.. The bivalent genome: characterization, structure, and regulation. Trends Genet. 2020; 36:118–131. [DOI] [PubMed] [Google Scholar]
- 5. Bernstein B.E., Mikkelsen T.S., Xie X., Kamal M., Huebert D.J., Cuff J., Fry B., Meissner A., Wernig M., Plath K.et al.. A bivalent chromatin structure marks key developmental genes in embryonic stem cells. Cell. 2006; 125:315–326. [DOI] [PubMed] [Google Scholar]
- 6. Azuara V., Perry P., Sauer S., Spivakov M., Jørgensen H.F., John R.M., Gouti M., Casanova M., Warnes G., Merkenschlager M.et al.. Chromatin signatures of pluripotent cell lines. Nat. Cell Biol. 2006; 8:532–538. [DOI] [PubMed] [Google Scholar]
- 7. Mikkelsen T.S., Ku M., Jaffe D.B., Issac B., Lieberman E., Giannoukos G., Alvarez P., Brockman W., Kim T.-K., Koche R.P.et al.. Genome-wide maps of chromatin state in pluripotent and lineage-committed cells. Nature. 2007; 448:553–560. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 8. Hanna R.E., Doench J.G.. Design and analysis of CRISPR-Cas experiments. Nat. Biotechnol. 2020; 38:813–823. [DOI] [PubMed] [Google Scholar]
- 9. Khan A., Giri S., Wang Y., Chakraborty A., Ghosh A.K., Anantharaman A., Aggarwal V., Sathyan K.M., Ha T., Prasanth K.V.et al.. BEND3 represses rDNA transcription by stabilizing a NoRC component via USP21 deubiquitinase. Proc. Natl. Acad. Sci. U.S.A. 2015; 112:8338–8343. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 10. Kurniawan F., Chetlangia N., Kamran M., Redon C.E., Pongor L., Sun Q., Lin Y.-C., Mohan V., Shaqildi O., Asoudegi D.et al.. BEND3 safeguards pluripotency by repressing differentiation-associated genes. Proc. Natl. Acad. Sci. U.S.A. 2022; 119:e2107406119. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 11. Saksouk N., Barth T.K., Ziegler-Birling C., Olova N., Nowak A., Rey E., Mateos-Langerak J., Urbach S., Reik W., Torres-Padilla M.-E.et al.. Redundant mechanisms to form silent chromatin at pericentromeric regions rely on BEND3 and DNA methylation. Mol. Cell. 2014; 56:580–594. [DOI] [PubMed] [Google Scholar]
- 12. Pailles M., Hirlemann M., Brochard V., Chebrout M., Oudin J.-F., Marks H., Jouneau A., Bonnet-Garnier A.. H3K27me3 at pericentromeric heterochromatin is a defining feature of the early mouse blastocyst. Sci. Rep. 2022; 12:13908. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 13. Zhang J., Zhang Y., You Q., Huang C., Zhang T., Wang M., Zhang T., Yang X., Xiong J., Li Y.et al.. Highly enriched BEND3 prevents the premature activation of bivalent genes during differentiation. Science. 2022; 375:1053–1058. [DOI] [PubMed] [Google Scholar]
- 14. Gupta N., Yakhou L., Richard Albert J., Azogui A., Ferry L., Kirsh O., Miura F., Battault S., Yamaguchi K., Laisné M.et al.. A genome-wide screen reveals new regulators of the 2-cell-like cell state. Nat. Struct. Mol. Biol. 2023; 30:1105–1118. [DOI] [PubMed] [Google Scholar]
- 15. Fukuda K., Okuda A., Yusa K., Shinkai Y.. A CRISPR knockout screen identifies SETDB1-target retroelement silencing factors in embryonic stem cells. Genome Res. 2018; 28:846–858. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 16. Wang C., Liu X., Gao Y., Yang L., Li C., Liu W., Chen C., Kou X., Zhao Y., Chen J.et al.. Reprogramming of H3K9me3-dependent heterochromatin during mammalian embryo development. Nat. Cell Biol. 2018; 20:620–631. [DOI] [PubMed] [Google Scholar]
- 17. Miura F., Enomoto Y., Dairiki R., Ito T.. Amplification-free whole-genome bisulfite sequencing by post-bisulfite adaptor tagging. Nucleic Acids Res. 2012; 40:e136. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 18. Miura F., Shibata Y., Miura M., Sangatsuda Y., Hisano O., Araki H., Ito T.. Highly efficient single-stranded DNA ligation technique improves low-input whole-genome bisulfite sequencing by post-bisulfite adaptor tagging. Nucleic Acids Res. 2019; 47:e85. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 19. Zagore L.L., Sweet T.J., Hannigan M.M., Weyn-Vanhentenryck S.M., Jobava R., Hatzoglou M., Zhang C., Licatalosi D.D.. DAZL regulates germ cell survival through a network of PolyA-proximal mRNA interactions. Cell Rep. 2018; 25:1225–1240. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 20. Rodriguez-Terrones D., Gaume X., Ishiuchi T., Weiss A., Kopp A., Kruse K., Penning A., Vaquerizas J.M., Brino L., Torres-Padilla M.-E.. A molecular roadmap for the emergence of early-embryonic-like cells in culture. Nat. Genet. 2018; 50:106–119. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 21. Borgel J., Guibert S., Li Y., Chiba H., Schübeler D., Sasaki H., Forné T., Weber M.. Targets and dynamics of promoter DNA methylation during early mouse development. Nat. Genet. 2010; 42:1093–1100. [DOI] [PubMed] [Google Scholar]
- 22. Hackett J.A., Reddington J.P., Nestor C.E., Dunican D.S., Branco M.R., Reichmann J., Reik W., Surani M.A., Adams I.R., Meehan R.R.. Promoter DNA methylation couples genome-defence mechanisms to epigenetic reprogramming in the mouse germline. Development. 2012; 139:3623–3632. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 23. Endoh M., Endo T.A., Shinga J., Hayashi K., Farcas A., Ma K.-W., Ito S., Sharif J., Endoh T., Onaga N.et al.. PCGF6-PRC1 suppresses premature differentiation of mouse embryonic stem cells by regulating germ cell-related genes. eLife. 2017; 6:e21064. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 24. Stielow B., Finkernagel F., Stiewe T., Nist A., Suske G.. MGA, L3MBTL2 and E2F6 determine genomic binding of the non-canonical polycomb repressive complex PRC1.6. PLoS Genet. 2018; 14:e1007193. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 25. Bindels D.S., Haarbosch L., van Weeren L., Postma M., Wiese K.E., Mastop M., Aumonier S., Gotthard G., Royant A., Hink M.A.et al.. mScarlet: a bright monomeric red fluorescent protein for cellular imaging. Nat. Methods. 2017; 14:53–56. [DOI] [PubMed] [Google Scholar]
- 26. Ferry L., Fournier A., Tsusaka T., Adelmant G., Shimazu T., Matano S., Kirsh O., Amouroux R., Dohmae N., Suzuki T.et al.. Methylation of DNA ligase 1 by G9a/GLP recruits UHRF1 to replicating DNA and regulates DNA methylation. Mol. Cell. 2017; 67:550–565. [DOI] [PubMed] [Google Scholar]
- 27. Petryk N., Bultmann S., Bartke T., Defossez P.-A.. Staying true to yourself: mechanisms of DNA methylation maintenance in mammals. Nucleic Acids Res. 2020; 49:3032–3020. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 28. Doench J.G., Fusi N., Sullender M., Hegde M., Vaimberg E.W., Donovan K.F., Smith I., Tothova Z., Wilen C., Orchard R.et al.. Optimized sgRNA design to maximize activity and minimize off-target effects of CRISPR-Cas9. Nat. Biotechnol. 2016; 34:184–191. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 29. Li W., Xu H., Xiao T., Cong L., Love M.I., Zhang F., Irizarry R.A., Liu J.S., Brown M., Liu X.S.. MAGeCK enables robust identification of essential genes from genome-scale CRISPR/Cas9 knockout screens. Genome Biol. 2014; 15:554. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 30. Karimi M.M., Goyal P., Maksakova I.A., Bilenky M., Leung D., Tang J.X., Shinkai Y., Mager D.L., Jones S., Hirst M.et al.. DNA methylation and SETDB1/H3K9me3 regulate predominantly distinct sets of genes, retroelements, and chimeric transcripts in mESCs. Cell Stem Cell. 2011; 8:676–687. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 31. Perez-Torrado R., Yamada D., Defossez P.-A.. Born to bind: the BTB protein-protein interaction domain. Bioessays. 2006; 28:1194–1202. [DOI] [PubMed] [Google Scholar]
- 32. Cuneo M.J., Mittag T.. The ubiquitin ligase adaptor SPOP in cancer. FEBS J. 2019; 286:3946–3958. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 33. Hoffman J.A., Wu C.-I., Merrill B.J.. Tcf7l1 prepares epiblast cells in the gastrulating mouse embryo for lineage specification. Development. 2013; 140:1665–1675. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 34. Li M., Yu J.S.L., Tilgner K., Ong S.H., Koike-Yusa H., Yusa K.. Genome-wide CRISPR-KO screen uncovers mTORC1-mediated Gsk3 regulation in naive pluripotency maintenance and dissolution. Cell Rep. 2018; 24:489–502. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 35. Mas G., Blanco E., Ballaré C., Sansó M., Spill Y.G., Hu D., Aoi Y., Le Dily F., Shilatifard A., Marti-Renom M.A.et al.. Promoter bivalency favors an open chromatin architecture in embryonic stem cells. Nat. Genet. 2018; 50:1452–1462. [DOI] [PubMed] [Google Scholar]
- 36. Zheng L., Liu J., Niu L., Kamran M., Yang A.W.H., Jolma A., Dai Q., Hughes T.R., Patel D.J., Zhang L.et al.. Distinct structural bases for sequence-specific DNA binding by mammalian BEN domain proteins. Genes Dev. 2022; 36:225–240. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 37. Karimi M., Johansson S., Ekström T.J.. Using LUMA: a Luminometric-based assay for global DNA-methylation. Epigenetics. 2006; 1:45–48. [DOI] [PubMed] [Google Scholar]
- 38. Grand R.S., Burger L., Gräwe C., Michael A.K., Isbel L., Hess D., Hoerner L., Iesmantavicius V., Durdu S., Pregnolato M.et al.. BANP opens chromatin and activates CpG-island-regulated genes. Nature. 2021; 596:133–137. [DOI] [PubMed] [Google Scholar]
- 39. Lu F., Liu Y., Jiang L., Yamaguchi S., Zhang Y.. Role of tet proteins in enhancer activity and telomere elongation. Genes Dev. 2014; 28:2103–2119. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 40. Manzo M., Wirz J., Ambrosi C., Villaseñor R., Roschitzki B., Baubec T.. Isoform-specific localization of DNMT3A regulates DNA methylation fidelity at bivalent CpG islands. EMBO J. 2017; 36:3421–3434. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 41. Gu T., Lin X., Cullen S.M., Luo M., Jeong M., Estecio M., Shen J., Hardikar S., Sun D., Su J.et al.. DNMT3A and TET1 cooperate to regulate promoter epigenetic landscapes in mouse embryonic stem cells. Genome Biol. 2018; 19:88. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 42. Li Y., Zheng H., Wang Q., Zhou C., Wei L., Liu X., Zhang W., Zhang Y., Du Z., Wang X.et al.. Genome-wide analyses reveal a role of Polycomb in promoting hypomethylation of DNA methylation valleys. Genome Biol. 2018; 19:18. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 43. Charlton J., Jung E.J., Mattei A.L., Bailly N., Liao J., Martin E.J., Giesselmann P., Brändl B., Stamenova E.K., Müller F.-J.et al.. TETs compete with DNMT3 activity in pluripotent cells at thousands of methylated somatic enhancers. Nat. Genet. 2020; 52:819–827. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 44. Dixon G., Pan H., Yang D., Rosen B.P., Jashari T., Verma N., Pulecio J., Caspi I., Lee K., Stransky S.et al.. QSER1 protects DNA methylation valleys from de novo methylation. Science. 2021; 372:eabd0875. [DOI] [PMC free article] [PubMed] [Google Scholar]
Associated Data
This section collects any data citations, data availability statements, or supplementary materials included in this article.
Supplementary Materials
Data Availability Statement
The accession number for the deep sequencing data reported in this paper are: GSE221710 (CRISPR screen, Mutational analysis) and GSE224077 (RNA-seq, WGBS, ChIP-seq, CUT&RUN-seq). The scripts used are available on Github: https://github.com/FumihitoMiura/Project-2 (permanent DOI: 10.6084/m9.figshare.23899266) and https://github.com/kirsho/DASH (permanent DOI: 10.6084/m9.figshare.23899254).