1. Stem Cells and Regenerative Medicine
  2. Chromosomes and Gene Expression
Download icon

DNA methylation directs genomic localization of Mbd2 and Mbd3 in embryonic stem cells

  1. Sarah J Hainer
  2. Kurtis N McCannell
  3. Jun Yu
  4. Ly-Sha Ee
  5. Lihua J Zhu
  6. Oliver J Rando
  7. Thomas G Fazzio  Is a corresponding author
  1. University of Massachusetts Medical School, United States
Research Article
  • Cited 5
  • Views 1,650
  • Annotations
Cite this article as: eLife 2016;5:e21964 doi: 10.7554/eLife.21964

Abstract

Cytosine methylation is an epigenetic and regulatory mark that functions in part through recruitment of chromatin remodeling complexes containing methyl-CpG binding domain (MBD) proteins. Two MBD proteins, Mbd2 and Mbd3, were previously shown to bind methylated or hydroxymethylated DNA, respectively; however, both of these findings have been disputed. Here, we investigated this controversy using experimental approaches and re-analysis of published data and find no evidence for methylation-independent functions of Mbd2 or Mbd3. We show that chromatin localization of Mbd2 and Mbd3 is highly overlapping and, unexpectedly, we find Mbd2 and Mbd3 are interdependent for chromatin association. Further investigation reveals that both proteins are required for normal levels of cytosine methylation and hydroxymethylation in murine embryonic stem cells. Furthermore, Mbd2 and Mbd3 regulate overlapping sets of genes that are also regulated by DNA methylation/hydroxymethylation factors. These findings reveal an interdependent regulatory mechanism mediated by the DNA methylation machinery and its readers.

https://doi.org/10.7554/eLife.21964.001

Introduction

In mammalian cells, 5-methylcytosine (5mC), a heritable epigenetic modification on DNA, occurs mainly at CpG dinucleotides through the actions of DNA methyltransferase (DNMT) enzymes (Bester, 1988). The members of the DNMT family that modify cytosine on DNA include Dnmt1, Dnmt3a, and Dnmt3b (Okano et al., 1998). Dnmt3a and Dnmt3b are de novo DNA methyltransferases, which establish methylation during development. On the other hand, the more abundant Dnmt1 predominantly methylates hemimethylated CpG dinucleotides to restore methylation patterns following genomic replication (although it also has some de novo methylation activity) and therefore, is considered to be the 'maintenance methyltransferase' (Egger et al., 2006; Kho et al., 1998; Pradhan et al., 1999; Probst et al., 2009). Dnmt1 is responsible for the majority of methylation in murine embryonic stem (ES) cells, as Dnmt1 knockout (KO) ES cells carry only about 20% of normal methylation levels (Lei et al., 1996).

Active demethylation of 5mC involves a relatively complex series of reactions that starts with oxidation by the ten-eleven translocation (TET) family of dioxygenases (including Tet1, Tet2, and Tet3; (Lu et al., 2015), which actively demethylate DNA by oxidizing the 5-methyl group of 5mC to form 5-hydroxymethylcytosine (5hmC) (Tahiliani et al., 2009). Further oxidation can occur through conversion of 5hmC into 5-formylcytosine (5fC) and 5-carboxylcytosine (5caC) (He et al., 2011; Ito et al., 2011). The level of 5hmC is approximately 10% the level of 5mC in ES cells (Tahiliani et al., 2009), whereas 5fC and 5caC are much less abundant (Ito et al., 2011). Together, Tet1 and Tet2 are responsible for essentially all the 5hmC present in ES cells (Koh et al., 2011). However, Tet1 is responsible for 5hmC production at promoter-proximal regions, whereas Tet2 is poorly chromatin-associated and primarily acts within gene bodies (Huang et al., 2014; Vella et al., 2013). In addition, knockdown (KD) or KO of Tet1 or Tet2 skews the profile of ES cell differentiation (Dawlaty et al., 2011; Ficz et al., 2011; Koh et al., 2011), and some reports suggest Tet1 KD also leads to a defect in self-renewal (Freudenberg et al., 2012; Ito et al., 2010). However, the functions of Tet proteins during development remain incompletely resolved as Tet1-/- Tet2-/-double KO mice have been shown to be viable with only modest phenotypes (Dawlaty et al., 2013).

The methylation status of DNA influences many biological processes during mammalian development, and cytosine methylation of promoter-proximal regulatory sequences in mammals is generally considered repressive for gene expression (Klose and Bird, 2006). One of the best-studied families of cytosine methylation effectors consists of the methyl-CpG-binding domain (MBD) family of proteins, most of which specifically bind methylated CpG dinucleotides and play important roles in determining the transcriptional state of the mammalian epigenome (Filion et al., 2006; Hendrich and Bird, 1998; Prokhortchouk et al., 2001; Unoki et al., 2004). MBD proteins are transcriptional regulators, primarily thought to act as repressors, which play a major role in coordinating crosstalk between cytosine methylation and chromatin structure to regulate transcription (Denslow and Wade, 2007). The MBD has the ability to bind single symmetrically methylated CpG dinucleotides (Nan et al., 1993; Ohki et al., 2001), and the functions of MeCP2 and the Mbd1-4 family members have been well-studied. With the exception of Mbd3, these proteins preferentially bind 5mC over unmethylated cytosine (Hendrich and Bird, 1998; Zhang et al., 1999). Consistent with this binding activity, Mbd2 occupies methylated CpG island–containing promoters of inactive genes (Chatagnon et al., 2011; Menafra et al., 2014). Although Mbd3 does not bind 5mC, MeCP2 and Mbd3 can bind 5hmC in vitro and Mbd3 is enriched in regions of elevated 5hmC in ES cells (Mellén et al., 2012; Yildirim et al., 2011). Consistent with this reported localization, knockdown (KD) of Mbd3 leads to misregulation of 5hmC-marked genes (Yildirim et al., 2011). Furthermore, Mbd3 binding in ESCs was strongly reduced upon RNAi-mediated KD of Tet1 (Yildirim et al., 2011) However, arguing against the above studies, other in vitro studies using short DNA probes containing a single symmetric 5hmCpG find poor binding of 5hmC by MBD family members (Cramer et al., 2014; Spruijt et al., 2013).

Mbd2 and Mbd3 are highly similar in amino acid sequence and are components of mutually exclusive versions of the nucleosome remodeling and deacetylase (NuRD) complex (Hendrich and Bird, 1998; Le Guezennec et al., 2006; Wade et al., 1999; Zhang et al., 1999). Furthermore, Mbd2 and Mbd3 exhibit partially overlapping localization profiles at some methylated regions in vivo (Günther et al., 2013), consistent with the possibility that these factors bind to DNA enriched for 5mC or its derivative, 5hmC. However, these highly similar complexes play distinct biological roles in vivo. Mbd3/NuRD is necessary for ES cell pluripotency and differentiation, as well as embryonic development, whereas Mbd2 KO mice are viable and fertile (Hendrich et al., 2001; Kaji et al., 2006; Reynolds et al., 2012). In addition, Mbd3/NuRD coordinates cytosine methylation by recruiting DNA methyltransferases to the promoters of tumor suppressor genes in colon cancer cell and leukemia cell lines (Cai et al., 2014; Choi et al., 2013; Morey et al., 2008), and depletion of Mbd3 results in reduced DNA methylation levels at some locations in ES cells (Latos et al., 2012).

Recently, evidence has arisen questioning the dependence of Mbd2 and Mbd3 on cytosine methylation for genomic localization (Baubec et al., 2013). The authors of this study reported that the enrichments of Mbd2 and Mbd3 at LMRs (low-methylated regions that are enriched for transcription factor-binding sites and exhibit approximately 30% methylation on average) (Stadler et al., 2011) were minimally altered in Dnmt1-/- Dnmt3a-/- Dnmt3b-/- triple knockout (TKO) ES cells. Here, we sought to resolve the conflicting data addressing the dependence of Mbd2 and Mbd3 localization on 5mC and 5hmC. Analyses of ChIP-seq data from Baubec et al. as well as multiple new ChIP-seq datasets reported here demonstrate methylation-dependence of Mbd2 and Mbd3 binding throughout the genome. Interestingly, we show that Mbd2 and Mbd3 exhibit significantly overlapping localization in vivo and find that Tet1 activity is required for normal chromatin association by both Mbd3 and Mbd2. Furthermore, we show that Mbd3 and Mbd2 are each required for the binding of the other, as well as for normal levels of 5mC and 5hmC. Finally, we find that individual KD of Mbd2, Mbd3, Dnmt1, or Tet1 results in highly concordant changes in gene expression. Together these data reveal interdependence among Mbd2/NuRD, Mbd3/NuRD, 5mC, and 5hmC to maintain normal chromatin structure and gene regulation in ESCs.

Results

Evidence of methylation-dependent localization of Mbd3 and Mbd2 in bio-MBD datasets

Multiple studies have established that Mbd2, but not Mbd3, binds methylated DNA with higher affinity than unmethylated DNA (Hendrich and Bird, 1998; Zhang et al., 1999). We previously found that Mbd3 occupies genomic regions containing 5hmC in ES cells (Yildirim et al., 2011). Furthermore, we found that KD of Tet1 results in a reduction of Mbd3 occupancy, and that KD of Mbd3 results in a reduction in bulk 5hmC levels (Yildirim et al., 2011), demonstrating interdependence between Mbd3 and 5hmC. These findings also lead to a straightforward prediction that cytosine methylation would be required for Mbd3 binding, since 5hmC is generated by oxidation of 5mC. However, Baubec et al. reported that both Mbd3 and Mbd2 occupancies at LMRs are minimally altered in TKO cells (Baubec et al., 2013), arguing against our previous findings regarding Mbd3 and long-standing models of Mbd2 function.

To assess the discrepancies between these studies, we first re-analyzed the datasets from Baubec et al. (Figure 1, Figure 1—figure supplement 1, and Figure 1—figure supplement 2). There were two ChIP-seq replicates for the biotin-tagged Mbd3a isoform (bio-Mbd3a) in wild-type (WT) cells and a single replicate reported for bio-Mbd3a ChIP-seq in TKO cells (see Table 1). We initially focused on genomic binding at the subset of LMRs that are larger than 150 bp and at least 3 kb away from another LMR or an unmethylated region (hereafter denoted as 'LMR subset') that was utilized in the prior study. We observed only a modest reduction (~28%) in bio-Mbd3a occupancy at the LMR subset upon loss of DNA methylation in TKO cells (Figure 1A), similar to the original observation (see Figure 7A of Baubec et al. (2013)). However, when we assessed the change in overall bio-Mbd3a enrichment by examining all bio-Mbd3a peaks (called from pooled WT bio-Mbd3a ChIP-seq reads), we observed a ~75% reduction in bio-Mbd3a binding in TKO cells relative to WT cells (Figure 1B–C and Figure 1—figure supplement 1A–B). These data demonstrate that Mbd3a relies heavily on DNA methylation in order to bind chromatin in ES cells and that confining analysis to a subset of LMRs obscures the global impact of methylation on Mbd3 localization.

Figure 1 with 2 supplements see all
Analysis of bio-Mbd3 and bio-Mbd2 ChIP-seq datasets reveals methylation-dependent chromatin localization in ES cells.

Mbd3 and Mbd2 ChIP-seq datasets from (Baubec et al., 2013) were obtained from GEO (GSE39610; see Table 1). (A–C) Analysis of ChIP-seq data for biotin-tagged Mbd3 (bio-Mbd3) in WT or TKO cells. The two available WT replicates and one available TKO replicate were mapped over the LMR subset (adapted from Stadler et al. [2011]) ± 2 kb (A) or bio-Mbd3 peaks from WT cells ± 2 kb (B). (C) Heatmaps for bio-Mbd3 ChIP-seq in WT or TKO cells over bio-Mbd3 peaks ± 2 kb, sorted by bio-Mbd3 occupancy (high to low) in replicate 1 WT cells. (D–H) Analysis of ChIP-seq data for biotin-tagged Mbd2 (bio-Mbd2) in WT or TKO cells. (D–E) Data were analyzed as described in (Baubec et al., 2013) where replicates were pooled and mapped over LMRs that are > 150 bp and at least 3 kb away from other LMRs or unmethylated regions (LMR subset) (D) or all LMRs (E). (F–G) Biological replicates were analyzed separately and mapped over the LMR subset ± 2 kb (F) or bio-Mbd2 peaks from WT cells ± 2 kb (G). (H) Heatmaps for bio-Mbd2 ChIP-seq in WT or TKO cells over bio-Mbd2 peaks ± 2 kb, sorted by bio-Mbd2 occupancy in replicate 1 WT cells. Peaks of bio-Mbd3 and bio-Mbd2 were called as described in Supplemental Experimental Procedures.

https://doi.org/10.7554/eLife.21964.002

While Baubec et al. found that bio-Mbd2 requires DNA methylation for proper binding to chromatin at highly methylated CpG islands (Figure 4G–H of Baubec et al. (2013)), they observed minimal effect of loss of methylation on bio-Mbd2 binding at the LMR subset (Figure 7A from Baubec et al. (2013)). To replicate this prior analysis, we pooled the mapped reads for each of three bio-Mbd2 ChIP-seq replicates from WT or two bio-Mbd2 ChIP-seq replicates from TKO cells (see Table 1) and measured their enrichment over the LMR subset. We observed a modest reduction (~15%) in bio-Mbd2 occupancy relative to WT (Figure 1D), consistent with the prior report. In contrast, when we analyzed bio-Mbd2 binding over all LMRs (referred to throughout as 'total LMRs'), we observed a stronger reduction (~40%) in bio-Mbd2 occupancy in the pooled TKO libraries relative to WT (Figure 1E).

Table 1

Related to Figure 1. SRA file numbers for re-analyzed Baubec et al. (2013) datasets.

https://doi.org/10.7554/eLife.21964.005

Name

SRA file number

Number of mapped reads (Bowtie, up to three mismatches)

WT Mbd3 ChIP replicate 1

SRR696667

11,601,021

WT Mbd3 ChIP replicate 2

SRR696673

14,166,151

 TKO Mbd3 ChIP

SRR769560

34,666,789

WT Mbd2 ChIP replicate 1

SRR527128

9,441,721

WT Mbd2 ChIP replicate 2

SRR527129

15,116,380

WT Mbd2 ChIP replicate 3

SRR696658

12,545,659

TKO Mbd2 ChIP replicate 1

SRR527161

24,282,816

TKO Mbd2 ChIP replicate 2

SRR696682

7,605,171

WT Mbd1a ChIP replicate 1

SRR527126

11,566,092

WT Mbd1a ChIP replicate 2

SRR527127

19,528,778

 TKO Mbd1a ChIP

SRR527159

22,070,693

WT Mbd1b ChIP replicate 1

SRR527147

14,320,842

WT Mbd1b ChIP replicate 2

SRR527148

22,371,797

 TKO Mbd1b ChIP

SRR527160

7,503,068

WT Mbd4 ChIP replicate 1

SRR527131

13,864,108

WT Mbd4 ChIP replicate 2

SRR527132

23,655,753

TKO Mbd4 ChIP replicate 1

SRR669321

13,670,155

TKO Mbd4 ChIP replicate 2

SRR669322

11,251,853

WT MeCP2 ChIP replicate 1

SRR527133

19,229,199

WT MeCP2 ChIP replicate 2

SRR527134

11,771,499

WT MeCP2 ChIP replicate 3

SRR696681

12,118,572

TKO MeCP2 ChIP replicate 1

SRR527162

17,003,233

TKO MeCP2 ChIP replicate 2

SRR696683

28,586,006

To better evaluate the effects of cytosine methylation on bio-Mbd2 binding, we analyzed each bio-Mbd2 ChIP-seq replicate separately. We observed high variation among the three replicates for bio-Mbd2 in WT cells (Figure 1F–H and Figure 1—figure supplement 1C–I). For example, while two WT replicates are modestly enriched for bio-Mbd2 binding relative to the two TKO replicates (~10% reduction) at the LMR subset, the third bio-Mbd2 ChIP-seq experiment in WT cells has a much higher peak, implying a much larger (~50%) reduction in occupancy in TKO cells (Figure 1F). Similar disparities were observed at other regions of the genome (Figure 1—figure supplement 1C–F). Importantly, when we examined the effect of 5mC loss over all peaks of bio-Mbd2 binding, bio-Mbd2 enrichment was reduced to near background levels in TKO cells (Figure 1G). Furthermore, we observed considerable gene-by-gene differences in the binding profiles among WT bio-Mbd2 replicates, where individual locations that are lowly bound in one replicate are highly bound in another replicate (Figure 1H, Figure 1—figure supplement 1G–I). Taken together, these analyses demonstrate bio-Mbd3 and bio-Mbd2 strongly depend on cytosine methylation for normal chromatin association. Furthermore, poor concordance between biological replicates was also found for several other MBD proteins profiled in Baubec et al. (Figure 1—figure supplement 2), preventing straightforward interpretation of the role of cytosine methylation in localization of these factors.

Dnmt1 and Tet1 are required for Mbd3 and Mbd2 occupancy in ES cells

To further address the discrepancies discussed above, and to better elucidate the roles of DNA methylation and hydroxymethylation in chromatin binding by Mbd3 and Mbd2, we performed biological duplicate ChIP-seq experiments for endogenous Mbd3 and Mbd2 in ES cells acutely depleted of Dnmt1 or Tet1 (Figure 2). Upon efficient KD of these factors (Figure 2—figure supplement 1A–B), we observed a significant reduction in promoter-proximal Mbd3 and Mbd2 binding relative to the binding observed in control (EGFP KD) cells (examples shown in Figure 2A, Figure 2—figure supplement 1C–D). The reductions in Mbd3 and Mbd2 binding observed upon Dnmt1 KD or Tet1 KD were not due to altered expression of these proteins (Figure 2—figure supplement 1E). When we compared the locations bound by Mbd3 and Mbd2 throughout the genome, we found that Mbd3 and Mbd2 co-occupy numerous regions throughout the genome of ES cells, and this overlap is particularly enriched at promoter-proximal regions (Figure 2B).

Figure 2 with 4 supplements see all
Dnmt1 and Tet1 are required for Mbd3 and Mbd2 binding in ES cells.

(A) Genome browser tracks of replicate ChIP-seq experiments examining endogenous Mbd3 or Mbd2 occupancy in control (EGFP KD), Dnmt1 KD, and Tet1 KD ES cells over indicated loci (Hoxd cluster). (B) Overlap of Mbd3 and Mbd2 binding. Shown are Venn diagrams delineating the overlap between all genomic locations (top panel) or TSS-proximal locations (−1 kb to +100 bp; bottom panel) bound by Mbd3 and Mbd2. (C–D) Aggregation plots of Mbd3 (C) or Mbd2 (D) ChIP-seq data showing occupancy over ICPs (top left panel, from [Weber et al., 2007]), annotated TSSs (bottom left panel), the LMR subset (top middle panel, from [Stadler et al., 2011]), gene-distal DHSs (bottom middle panel, from GSM1014154 with TSSs removed) Mbd2 peaks (top right panel, called from EGFP KD ChIP-seq experiments), and Mbd3 peaks (bottom right panel, called from EGFP KD ChIP-seq experiments) ± 2 kb in control (EGFP KD), Dnmt1 KD, or Tet1 KD ES cells. (E–F) Heatmaps of Mbd3 enrichment over Mbd3 binding sites sorted by Mbd2 occupancy (E) and Mbd2 enrichment over Mbd2 binding sites sorted by Mbd3 occupancy (F) in control (EGFP KD), Dnmt1 KD, or Tet1 KD ES cells. The profiles shown in aggregation plots and heatmaps represent the average of two biological replicates.

https://doi.org/10.7554/eLife.21964.006

Next, we quantified average binding profiles of Mbd3 and Mbd2 over six genomic features – intermediate CpG density promoters (ICPs; from (Weber et al., 2007), transcription start sites (TSSs), the LMR subset described above, gene-distal DNaseI hypersensitive sites (DHSs; from ENCODE; GSM1014154), Mbd2 peaks (from EGFP KD cells), and Mbd3 peaks (from EGFP KD cells). We observed a similar pattern over each genomic feature: Dnmt1 KD and Tet1 KD each resulted in a dramatic reduction in Mbd3 and Mbd2 occupancies (Figure 2C–F). Importantly, our duplicate ChIP-seq datasets were highly reproducible (Figure 2A; Figure 2—figure supplement 1C–D and Figure 2—figure supplement 2A–D). We validated our observations at eight genomic locations by performing three biological replicate ChIP-qPCRs for Mbd3 and Mbd2 in Dnmt1 KD and Tet1 KD ES cells (Figure 2—figure supplement 2E) and in Dnmt1 KO cells (Figure 2—figure supplement 2F–G).

Four key controls demonstrate the specificity of these ChIP data. First, we performed ChIP-qPCR for Mbd3 in Mbd3 KD cells and for Mbd2 in Mbd2 KD cells, confirming the reduced ChIP signals expected in both cases (Figure 2—figure supplement 2E). To validate our findings in a manner independent of the antibodies targeting the endogenous proteins, we genetically FLAG-tagged endogenous Mbd2 (hereafter Mbd2-3XFLAG) or Mbd3 ([Yildirim et al., 2011]; hereafter Mbd3abc-3XFLAG) and performed ChIP-qPCR or ChIP-seq, respectively, in EGFP KD, Dnmt1 KD, and Tet1 KD cells (Figure 2—figure supplement 3A–D). We again observed a reduction in Mbd2-3XFLAG and Mbd3abc-3XFLAG enrichment following knockdown of Dnmt1 or Tet1 (Figure 2—figure supplement 3B–D), despite lower overall enrichment than we observed with antibodies against the endogenous proteins. FLAG ChIP-seq and ChIP-qPCR for Mbd3abc-3XFLAG in Mbd3 KD cells and for Mbd2-3XFLAG in Mbd2 KD cells confirmed the reduced ChIP signals expected in both cases (Figure 2—figure supplement 3C–D; Figure 2—figure supplement 4A–B). Finally, to test whether overexpression of Mbd3a (analogous to the bio-Mbd3a ChIP-seq studies in [Baubec et al., 2013]) resulted in altered genomic binding relative to endogenously expressed Mbd3, and whether chromatin binding by overexpressed Mbd3a was altered upon Dnmt1 KD or Tet1 KD, we overexpressed Mbd3a-3XFLAG in ES cells and performed ChIP-qPCR experiments upon KD of Dnmt1, Tet1, or EGFP (Figure 2—figure supplement 3C). Although Mbd3a-3XFLAG enrichment was reduced relative to endogenously expressed Mbd3abc-3XFLAG, there was a decrease in Mbd3a-3XFLAG binding upon Dnmt1 or Tet1 depletion at all eight promoter-proximal regions examined (Figure 2—figure supplement 3C). Together, these data demonstrate a requirement for Dnmt1 and Tet1 for Mbd3 and Mbd2 binding throughout the genome of ES cells. Although the dependence of Mbd3 and Mbd2 on DNA methylation (or Dnmt1) and the requirement of Mbd3 for hydroxymethylation (or Tet1) were consistent with our previous findings and those of others (Chatagnon et al., 2011; Menafra et al., 2014; Yildirim et al., 2011), the requirement of Mbd2 for Tet1 has not been previously observed. We will address this finding below.

One argument suggesting the data in Yildirim et al. may not be valid was that Mbd3 binding was not observed at enhancers (Menafra and Stunnenberg, 2014). In this study, we used H3K4me1 as a marker for enhancer regions (Yildirim et al., 2011), which is an imperfect marker for putative enhancer elements (Shlyueva et al., 2014). We therefore re-analyzed the Yildirim data at gene-distal DNaseI hypersensitive sites (DHSs), which represent active enhancers and some additional regulatory features, and observed a peak of Mbd3 binding that was strongly reduced upon Tet1 depletion (Figure 2—figure supplement 4C). Similarly, our new Mbd3 and Mbd2 ChIP-seq datasets exhibit robust enrichment at gene-distal DHSs that is reduced upon Dnmt1 KD or Tet1 KD (Figure 2C–D), confirming Mbd3 binds to these regions in a manner dependent on the cytosine methylation and hydroxymethylation machinery.

The catalytic activity of Tet1 is required for proper chromatin localization by Mbd3 and Mbd2 in ES cells

Although our KD studies demonstrated Tet1 is necessary for chromatin binding by Mbd3 and Mbd2, the Tet1 protein has been shown, in some cases, to regulate gene expression independently of its catalytic activity (Jin et al., 2014; Kaas et al., 2013; Williams et al., 2011). These findings leave open the possibility that Tet1 functions in Mbd3 and Mbd2 binding independently of its roles in DNA hydroxymethylation and demethylation. To test whether the catalytic activity of Tet1 is necessary for Mbd3 and Mbd2 binding in ES cells, we generated homozygous catalytically inactive Tet1 mutant (Tet1c.i.) ES cells. We knocked-in two amino acid substitutions, H1652Y and D1654A (Figure 3—figure supplement 1), which were previously demonstrated to eliminate DNA demethylation activity by Tet1 (Ito et al., 2010; Tahiliani et al., 2009). We performed ChIP-seq for Mbd3 and Mbd2 in two independent Tet1c.i. lines (Figure 3). We found that relative to WT ES cells, ChIP-seq of Mbd3 and Mbd2 in Tet1c.i. cells shows decreased occupancies of both Mbd3 and Mbd2 at ICPs, TSSs, the LMR subset, gene-distal DHSs, Mbd2 peaks, and Mbd3 peaks (Figure 3). These data demonstrate that the catalytic activity of Tet1 is required to promote binding of Mbd3 and Mbd2 to multiple genomic locations in ES cells.

Figure 3 with 1 supplement see all
The catalytic activity of Tet1 is required for occupancy of Mbd3 and Mbd2 in ES cells.

Aggregation plots of Mbd3 or Mbd2 ChIP-seq showing association over ICPs (top left panel), annotated TSSs (bottom left panel), the LMR subset (top middle panel), gene-distal DHSs (bottom middle panel), Mbd2 peaks (top right panel), and Mbd3 peaks (bottom right panel) ± 2 kb in wild-type (WT) and Tet1 catalytically inactive (Tet1c.i.) ES cells. Aggregation plots represent the average of two biological replicates.

https://doi.org/10.7554/eLife.21964.011

Loss of Mbd3 or Mbd2 results in reduced 5mC and 5hmC at regulatory regions

In order to better understand the dependence of Mbd3 and Mbd2 on cytosine methylation and hydroxymethylation, we performed duplicate MeDIP-seq and hMeDIP-seq experiments (which enrich for DNA harboring 5mC or 5hmC, respectively) on EGFP KD, Dnmt1 KD, Mbd3 KD, and Mbd2 KD ES cells (Figure 4A–B). We examined the effects of depletion of these factors on methylation (Figure 4A) and hydroxymethylation (Figure 4B) at ICPs, TSSs, the LMR subset, gene-distal DHSs, Mbd2 peaks, and Mbd3 peaks. As expected, Dnmt1 KD resulted in decreased 5mC and 5hmC. Furthermore, both Mbd3 KD and Mbd2 KD also resulted in decreased 5mC and 5hmC levels. These data are consistent with previous findings showing that Mbd3 KO cells have decreased 5mC levels at some locations (Latos et al., 2012), and reduced 5hmC levels overall (Yildirim et al., 2011) in ES cells. Interestingly, when we examined the methylation and hydroxymethylation profiles over total LMRs (Figure 4—figure supplement 1A–B), we observed a much stronger peak of enrichment in control (EGFP KD) ES cells than found at the LMR subset. Further investigation revealed that the LMRs excluded from the LMR subset were overrepresented at promoter proximal regions of genes (Figure 4—figure supplement 1E), and exhibited higher average methylation and hydroxymethylation (Figure 4—figure supplement 1C–D).

Figure 4 with 1 supplement see all
Dnmt1, Mbd3, and Mbd2 are required for 5mC and 5hmC in ES cells.

(A–B) Aggregation plots of 5mC (A) or 5hmC (B) enrichment over ICPs (top left panel), annotated TSSs (bottom left panel), the LMR subset (top middle panel), gene-distal DHSs (bottom middle panel) Mbd2 peaks (top right panel), and Mbd3 peaks (bottom right panel) ± 2 kb in control (EGFP KD), Dnmt1 KD, Mbd3 KD, or Mbd2 KD ES cells. Aggregation plots represent the average of two biological replicates. (C–D) Dot blot of bulk 5mC (C) and 5hmC (D) levels in EGFP KD, Dnmt1 KD, Mbd3 KD, or Mbd2 KD ES cells. Left panel shows 5mC or 5hmC levels and right panel shows methylene blue staining, which serves as a loading control. (E) Thin layer chromatography separation of radioactively end-labeled bases from MspI digested genomic DNA. Quantification of 5mC or 5hmC is expressed relative to T in each KD. Levels in EGFP KD are set to 100%. (F) Restriction enzyme digest of genomic DNA from EGFP KD, Dnmt1 KD, Mbd3 KD, or Mbd2 KD ES cells using an enzyme blocked by CpG methylation/hydroxymethylation (HpaII) or a CpG methylation insensitive enzyme (MspI) that cuts the same site (CCGG).

https://doi.org/10.7554/eLife.21964.013

To validate these data, we performed dot blotting, thin layer chromatography (TLC), and digests with methylation sensitive (HpaII) or insensitive (MspI) restriction enzymes on DNA isolated from Mbd3 KD, Mbd2 KD, or control cells (EGFP KD and Dnmt1 KD). Consistent with the MeDIP- and hMeDIP-seq data, bulk 5mC and 5hmC levels were reduced upon Mbd3 KD or Mbd2 KD, as measured by dot blotting and TLC, while digestion with HpaII was enhanced (Figure 4C–F). Together these data demonstrate that levels of cytosine methylation and hydroxymethylation are regulated not only by the DNA methyltransferase enzymes, but also by factors that depend on these modifications for proper chromatin localization.

Interdependence of Mbd3 and Mbd2 occupancy

The reduction in 5mC and 5hmC levels observed upon depletion of Mbd3 and Mbd2 raises the possibility that Mbd3 and Mbd2 may each be necessary (either directly or indirectly) for normal chromatin binding by the other. To address this possibility, we investigated the effect of Mbd3 KD or Mbd2 KD on chromatin association by Mbd2 or Mbd3, respectively. Relative to control (EGFP KD) cells, we found that when Mbd3 was depleted there was a decrease in Mbd2 occupancy, and, conversely, when Mbd2 was depleted there was a decrease in Mbd3 occupancy (example loci shown in Figure 5A; Figure 5—figure supplement 1A–B). We extended this analysis to the whole genome, confirming a significant reduction in binding of both Mbd3 and Mbd2 upon KD of the other MBD protein at ICPs, TSSs, the LMR subset, gene-distal DHSs, Mbd2 peaks, and Mbd3 peaks (Figure 5B–C). Depletion of Mbd3 does not alter Mbd2 levels, nor does depletion of Mbd2 alter levels of Mbd3 (Figure 2—figure supplement 1E). We validated our observations by performing three biological replicate ChIP-qPCRs for Mbd3 and Mbd2 in Mbd3 KD and Mbd2 KD ES cells (Figure 2—figure supplement 2E) and by ChIP-seq for Mbd3abc-3XFLAG and Mbd2-3XFLAG in Mbd3 KD and Mbd2 KD ES cells (Figure 2—figure supplement 4A–B). Together these data demonstrate interdependence between Mbd3 and Mbd2, whereby these mutually exclusive constituents of the NuRD complex regulate the genomic localization of one another.

Figure 5 with 2 supplements see all
Mbd3 and Mbd2 are required for each other’s binding in ES cells.

(A) Genome browser tracks of ChIP-seq experiments examining Mbd3 or Mbd2 occupancy in control (EGFP KD), Mbd2 KD, or Mbd3 KD ES cells over one example locus (Pitx1). (B) Aggregation plots of Mbd3 or Mbd2 ChIP-seq data showing occupancy over ICPs (top left panel), annotated TSSs (bottom left panel), LMR subset (top middle panel), gene-distal DHSs (bottom middle panel), Mbd2 peaks (top right panel), and Mbd3 peaks (bottom right panel) ± 2 kb in control (EGFP KD), Mbd3 KD, or Mbd2 KD ES cells. (C) Heatmaps of Mbd3 enrichment over Mbd3 binding sites ± 2 kb sorted by Mbd2 occupancy (top panel) and Mbd2 enrichment over Mbd2 binding sites ± 2 kb sorted by Mbd3 occupancy (bottom panel) in control (EGFP KD), Mbd2 KD, or Mbd3 KD ES cells. The profiles shown in the browser tracks, aggregation plots, and heatmaps represent the average of two biological replicates.

https://doi.org/10.7554/eLife.21964.015

Mbd3/NuRD was previously shown to play a role in methylation of promoter-proximal regions in cancer cell lines by helping recruit DNA methyltransferase enzymes to specific loci (Cai et al., 2014; Choi et al., 2013; Morey et al., 2008). To investigate whether this potential mechanism of regulation was also operational in ES cells, we examined the binding profile of Dnmt1 in Mbd3 KO cells (Figure 5—figure supplement 2). We found that Mbd3 KO cells (Figure 5—figure supplement 2A) showed a general loss of Dnmt1 binding, with Dnmt1 occupancy decreased relative to WT cells at all genomic locations examined (Figure 5—figure supplement 2B). These data show Dnmt1 relies on Mbd3 for proper binding. Therefore, the decreased methylation, hydroxymethylation, and Mbd2 binding observed upon depletion of Mbd3 may be attributed (at least in part) to loss of chromatin binding by Dnmt1.

Overlapping effects of Dnmt1, Tet1, Mbd3, or Mbd2 loss on gene expression

The various interdependencies detailed above make an additional functional prediction – that Dnmt1, Tet1, Mbd3, and Mbd2 should have overlapping roles in gene regulation in ES cells. To test this, we used DNA microarrays to examine the changes in mRNA levels upon Dnmt1 KD, Tet1 KD, Mbd3 KD, Mbd2 KD, and (since it is a key constituent of both NuRD complexes in ES cells) Chd4 KD. We observed well-correlated gene expression profiles in ES cells knocked down individually for these factors (Figure 6—figure supplement 1, Table 3), suggesting a significant overlap in their sets of target genes, although Tet1 KD and Mbd2 KD generally had weaker effects on gene expression (Figure 6A), consistent with the less severe phenotypes observed upon Mbd2 loss (Günther et al., 2013; Hendrich et al., 2001). Together, these data are consistent with a model in which these chromatin remodeling factors and DNA modifications coordinate to control gene activity in ES cells.

Figure 6 with 1 supplement see all
KD of DNA methylation/hydroxymethylation machinery or its readers results in similar effects on gene expression.

(A) K-means clustering (k = 3) of genes misregulated (adjusted p<0.05) upon Dnmt1 KD, Tet1 KD, Mbd3 KD, Mbd2 KD, or Chd4 KD. Upregulated genes are indicated in yellow and downregulated genes are indicated in blue. (B) Model for interdependent regulatory mechanism mediated by DNA methylation/hydroxymethylation and readers of 5mC/5hmC. The DNA methylation/hydroxymethylation machinery (Dnmt1 and Tet1) and its readers (Mbd3 and Mbd2) form a regulatory loop where disruption of one factor results in altered methylation patterns and MBD binding. Binding of a gene activator to unmethylated site, loss of methylation by another mechanism, or loss of a reader of 5mC/5hmC can break the regulatory loop. YFG (your favorite gene) indicates a generic gene body.

https://doi.org/10.7554/eLife.21964.018

Discussion

The role of cytosine methylation in the regulation of Mbd3 and Mbd2 binding to chromatin has been controversial. Here, we sought to address these conflicting data. Several differences in study design may help explain some of the discrepancies between our previous studies and those reported by Baubec et al. First, whereas endogenous Mbd3 was examined in ChIP studies performed in Yildirim et al., the bio-tagged Mbd3 in Baubec et al. was overexpressed [see Figure S1F from Baubec et al. (2013)]. Since overexpression of proteins can promote promiscuous binding as assayed by ChIP (Baresic et al., 2014; Fernandez et al., 2003), it is possible some of the peaks of Mbd3 binding identified by this method are non-physiological. Second, all Mbd3 isoforms were immunoprecipitated in the Yildirim study, rather than the single largest isoform examined in Baubec et al. Third, whereas Yildirim et al. examined Mbd3 binding to promoter-proximal regions of genes, Baubec et al. focused mainly on a subset of LMRs and ICPs. Regardless of these differences, comprehensive analysis of the data available from Baubec et al. revealed a strong reduction in Mbd3 binding in TKO cells at the majority of bio-Mbd3 binding sites (Figure 1B–C). Similarly, our re-analysis of the bio-Mbd2 ChIP-seq data also revealed a strong reliance on DNA methylation (Figure 1D–E).

To more thoroughly address the role of DNA methylation in the regulation of Mbd3 and Mbd2 binding in ES cells, we performed a series of ChIP-seq experiments for endogenous Mbd3 and Mbd2 in Dnmt1 or Tet1 KD cells. Duplicate ChIP-seq experiments performed with polyclonal antibodies targeting the endogenous MBD proteins (Figure 2), independent triplicate ChIP-qPCR experiments performed with the same polyclonal antibodies in Dnmt1 KD, Tet1 KD, or Dnmt1 KO cells (Figure 2—figure supplement 2), and ChIP experiments performed on endogenously C-terminal FLAG-tagged versions of Mbd2 or Mbd3 (Figure 2—figure supplement 3) all yielded the same result: depletion of Dnmt1 or Tet1 results in reduced binding of Mbd3 and Mbd2 throughout the genome. Furthermore, experiments utilizing point mutations that eliminate the catalytic activity of Tet1 demonstrated the DNA demethylase activity of Tet1 is necessary for Mbd3 and Mbd2 binding (Figure 3) and (along with data from Dnmt1 KO cells) show that our findings are not an artifact of potential off-target effects from RNAi. These data, together with a re-evaluation of published data, conclusively demonstrates the requirement for DNA methylation/hydroxymethylation for normal chromatin binding by Mbd3 and Mbd2 in ES cells.

Having confirmed a role for cytosine methylation in Mbd3 and Mbd2 binding to DNA in ES cells, we sought to further investigate the mechanism through which DNA methylation regulates binding of these factors. Previous studies revealed that Mbd3 is required for 5mC at some locations in ES cells (Latos et al., 2012), but may prevent methylation at other sites (Cui and Irudayaraj, 2015). Our data examining 5mC and 5hmC enrichment genome-wide reveals a decrease in global 5mC and 5hmC levels in ES cells depleted of Mbd3 or Mbd2 (Figure 4). This decrease is likely due to the reduction in Dnmt1 binding throughout the genome we observed in Mbd3 KO cells (Figure 5—figure supplement 2), consistent with findings in cancer cell lines (Cai et al., 2014; Choi et al., 2013; Morey et al., 2008), although other mechanisms may also contribute to this phenomenon.

Consistent with previous studies (Günther et al., 2013), we found that Mbd3 and Mbd2 binding overlap at many locations. This observation, along with the fact that Mbd3 is necessary for normal DNA methylation, suggested Mbd3 might be indirectly necessary for chromatin binding by Mbd2. We therefore investigated the potential interplay between these factors and found that Mbd2 requires Mbd3 for binding, and vice versa, throughout the genome (Figure 5). Together these data support a model in which Mbd3 and Mbd2 rely on DNA methylation and hydroxymethylation for binding, while these factors are also required for DNA methylation and hydroxymethylation. Such a model must be consistent with the following observations: (1) DNA methylation is necessary for proper Mbd2 and Mbd3 localization (as shown by Dnmt1 KD and KO experiments described in this study). (2) DNA demethylation is required for proper Mbd2 and Mbd3 localization (Tet1 KD and Tet1c.i. studies here, and data from Yildirim et al. (2011)). (3) Mbd2 binds DNA harboring 5mC (Hendrich and Bird, 1998; Zhang et al., 1999) and Mbd3 binds DNA binding 5hmC much more strongly than DNA containing 5mC (Yildirim et al., 2011). (4) Mbd2 and Mbd3 are required for normal levels of both 5mC and 5hmC (MeDIP/hMeDIP, TLC, and other experiments described in this study). (5) Dnmt1 localization depends on Mbd3 (ChIP studies described in this study). (6) RNAi-mediated depletion of Mbd2, Mbd3, Dnmt1, Tet1, or Chd4 results in misregulation of many of the same genes (expression profiling described in this study). Therefore, we propose that these two MBD proteins and DNA modifications all contribute to a regulatory loop that modulates gene expression (Figure 6B). Consequently, upon disruption of any component of this regulatory loop (through depletion, KO, or mutation of the catalytic domain), the gene regulatory outcome is similar in ES cells: binding of both NuRD complexes is reduced, altering gene expression. While the molecular mechanisms underlying these interdependencies are unclear, physical interactions between Mbd3 and DNA methylation/hydroxymethylation proteins have previously been observed (Cai et al., 2014; Yildirim et al., 2011). Furthermore, predictive studies utilizing published ChIP-seq datasets suggest interdependency between the Mbd2 and Tet1 proteins in ES cells (Juan et al., 2016).

How widespread are the interdependencies among Mbd2, Mbd3, and the DNA modification machinery? Mbd2 KO mice are viable and fertile, whereas Mbd3 KO mice are lethal (Hendrich et al., 2001; Kaji et al., 2006; Reynolds et al., 2012). Mbd3 is required for ES cell differentiation, suggesting Mbd3 may be sufficiently functional in the absence of Mbd2 for its essential roles in vivo. Consistent with this possibility, Mbd3 KD had a stronger effect on some shared Mbd2/Mbd3 target genes than Mbd2 KD in ES cells. Alternatively, Mbd2/Mbd3 interdependence may be lost during differentiation, or within specific lineages, which may partially account for the different phenotypes of Mbd2 and Mbd3 KO mice. Further studies will be required to distinguish between these possibilities.

Materials and methods

Cell culture

Mouse ES cells were derived from E14 (Hooper et al., 1987 RRID:CVCL_C320) and were cultured as previously described (Chen et al., 2013). Cells have been verified that they are of male mouse origin through sequencing performed in this and previous studies, and were previously tested to ensure they were free of mycoplasma. RNAi-mediated KD was performed for 48 hr with endoribonuclease III digested siRNAs (esiRNAs) as previously described (Fazzio et al., 2008). The Mbd3abc-3XFLAG cell line was described previously (Yildirim et al., 2011).

Generation of cell lines

To generate the Mbd3a-3XFLAG lines, the mouse Mbd3a coding sequence was amplified from cDNA, a 6XHis + 3 XFLAG tag was added to the 3’ end, and it was cloned into pLJM1-EGFP (courtesy of David Sabatini’s lab; Addgene plasmid # 19319) at the AfeI and BstBI sites, replacing the EGFP gene. 293T cells were transfected with pLJM1-Mbd3a-3XFLAG and Lenti-X plasmids (Clontech, Mountain View, CA, USA) via X-tremeGENE (Roche, Branford, CT, USA) to produce lentivirus. Virus was collected three times between 40 and 64 hr and concentrated as described previously (Chen et al., 2013). Virus was introduced to E14 ES cells with 8 µg/mL hexadimethrine bromide (Polybrene; Sigma, St. Louis, MO, USA). Cells were selected with puromycin, diluted to single clones, and screened for Mbd3a-3XFLAG expression by western blotting of whole cell lysates.

To generate the endogenously tagged Mbd2-3XFLAG lines, sequence corresponding to Mbd2 exon six was synthesized (IDT, Coralville, IA, USA) with a 6XHis + 3 XFLAG tag before the stop codon and cloned into pBluescript before an internal ribosome entry site followed by a NeoR gene and polyadenylation site, with ~2 kb of homology to either side of the Mbd2 target site, which was amplified from genomic DNA. Oligonucleotides (F: CACCGTGAGGCGTAAGAATATGATC, R: AAACGATCATATTCTTACGCCTCAC) were cloned into pX330-puro(Cong et al., 2013) for expression of a guide RNA that targets Cas9 to the C-terminus of Mbd2, in order to increase efficiency of targeted homologous recombination with the homology template. The guide RNA was designed over the stop codon of the gene so that it could not target the homology construct or a recombined Mbd2-3XFLAG. The two plasmids were transfected together into E14 ES cells by electroporation, and the cells were subsequently selected with G418. Clones were screened by PCR on genomic DNA and by Western blot for FLAG-tagged Mbd2 on whole cell lysates.

Dnmt1 and Mbd3abc knockout (KO) lines were created by CRISPR/Cas9 targeted double-strand breaks and error-prone DNA repair to generate frameshift mutations (Cong et al., 2013; Wang et al., 2013). Oligonucleotides for guide RNAs specific to Dnmt1 exons 15 (F: CACCGGGTTAGGGTCGTCTAGGTGC, R: AAACGCACCTAGACGACCCTAACCC) or 20 (F: CACCGAAACTGGGCGTGGCGTAAGA, R: AAACTCTTACGCCACGCCCAGTTTC) or to Mbd3 exon 5 (F: CACCGGTGTGTAGAGCACTCGCAA, R: AAACTTGCGAGTGCTCTACACACC) were cloned into pX330-puro. The plasmids were separately transfected into E14 ES cells by electroporation (for Dnmt1 KO) or using FuGENE HD (Promega, Madison, WI, USA, for Mbd3abc KO) and clones were selected with puromycin. Candidates were screened by Western blotting of Dnmt1 or Mbd3 protein from whole cell lysates and by sequencing of PCRs on genomic DNA at the target sites.

Tet1 catalytically inactive (Tet1c.i.) lines were generated by CRISPR/Cas9 targeted homologous recombination. Oligonucleotides (F: CACCGATTTTTGTGCCCATTCTCACA R: AAACTGTGAGAATGGGCACAAAAATC) were cloned into pX330-puro for a guide RNA that targets the Tet1 catalytic site for Cas9 cleavage. A sequence harboring the Tet1 H1652Y and D1654A mutations (previously demonstrated to eliminate Tet1 catalytic activity (Ito et al., 2010; Tahiliani et al., 2009)) was synthesized (IDT) with ~500 bp of homology to the genomic DNA on each side and cloned into pCR2.1 (Invitrogen, Grand Island, NY, USA). The two plasmids were transiently transfected together into E14 ES cells with FuGENE HD (Promega), and the cells were subsequently puromycin selected for ~40 hr to enrich for transfected cells. Clones were screened for incorporation of the mutations at both alleles by specific digestion by PsiI (NEB, Ipswich, MA, USA) of PCRs from genomic DNA. PCR products were sequenced. Western blots of whole cell lysates from the lines were performed to ensure that mutant Tet1 protein levels were unchanged relative to wild-type levels.

RT-qPCR

RNA was isolated from cells using TRIzol (Invitrogen) and used to synthesize cDNA with random hexamers (Promega) as described (Hainer et al., 2015). cDNA was used in quantitative PCR reactions with Dnmt1, Tet1, Mbd3, Mbd2, Chd4, or GAPDH specific primers (see Table 2) and a FAST SYBR mix (KAPA Biosystems, Woburn, MA, USA) on an Eppendorf Realplex.

Table 2

Related to Figure 2—figure supplement 13. Primer list for qPCR.

https://doi.org/10.7554/eLife.21964.020

Gene

Forward

Reverse

Dnmt1

TGTTCTGTCGTCTGCAACCT

GCCATCTCTTTCCAAGTCTTT

Tet1

TCACAGGCACAGGTTACAAAAG

TCCTTACATTTTCAAGGGGATG

Mbd3

GGCCACAGGGATGTCTTTTACT

CTTGACCTGGTTGGAAGAATCA

Mbd2

AACCAAATTCACGAACCACC

CCTTGTAGCCTCTTCTCCCA

Chd4

AAGTTTGCAGAGATGGAAGAGC

GGTCGTAGTCCTGAATCTCCAC

Farp1

AACTGCAAGTCATTCTAAATCTCG

GGTATTCAATGCCAGAGACACA

Ppp2r2c

CGAATTATCCAGCTCTGCCTTA

TGGAGGAGAGACTTAGGGGTGT

Gpr83

GAGCCACCTTACTGTAGGGAATG

CACGCTCACCAGCTTTCTGTA

Cldn7

ACCTTTGGAAGAGCAGTCAGTG

CCTTCTCCATCCACACACTTTC

Tgfb1

AAGTCAGAGACGTGGGGACTTCTTG

AGTCTTCGCGGGAGGCGGGGT

Trh

TAATGCCTCTGACCTGGGATC

CCCACATCCTAATTCCAAAGTG

Table 3

Related to Figure 6−figure supplement 1. Correlation coefficients for pairwise combinations of factor KD. Tip60 KD (a histone acetyltransferase with repressive functions in ES cells) and Brg1 KD (an ATP-dependent chromatin remodeling enzyme) are included for comparison.

https://doi.org/10.7554/eLife.21964.021

Dnmt1

Tet1

Mbd3

Mbd2

Chd4

Tip60

Brg1

KD:

1.0

0.508

0.693

0.612

0.645

-0.033

0.051

Dnmt1

1.0

0.441

0.551

0.475

-0.003

0.104

Tet1

1.0

0.712

0.831

-0.101

0.308

Mbd3

1.0

0.656

-0.027

0.033

Mbd2

1.0

-0.094

0.269

Chd4

1.0

0.199

Tip60

1.0

Brg1

Western blotting

Western blotting was performed as previously described (Hainer et al., 2015). From whole cell lysates, 30 μg of protein were separated by SDS-PAGE, transferred to nitrocellulose (Life Sciences, St Petersburg, FL, USA), and assayed by immunoblotting. The antibodies used to detect proteins were: anti-Dnmt1 (1:1000, Sigma D4692, RRID:AB_262096), anti-Tet1 (1:1000, Millipore, Billerica, MA, USA, 09–872, RRID:AB_10806199), anti-Mbd3 (1:1000, Bethyl Labs, Montgomery, TX, USA, A302-529A; RRID:AB_1998976), anti-Mbd2 (1:1000, Bethyl Labs A301-632A, RRID:AB_1211478), anti-Chd4 (1:1000, Bethyl Labs A301-082A, RRID:AB_873002), anti-FLAG M2 (1:10,000, Sigma F1804, RRID:AB_262044) and anti-β-actin (1:50,000, Sigma A1978, RRID:AB_476692).

Chromatin immunoprecipitation (ChIP)-qPCR and ChIP-seq

ChIP experiments were performed as previously described (Hainer et al., 2015). Ten million cells were fixed, washed with ice-cold PBS, pelleted, lysed through sonication on high in a Bioruptor UCD-200 (Diagenode, Delville, NJ, USA) in SDS Lysis Buffer (1% SDS, 10 mM EDTA, 50 mM Tris-HCl pH 8.0), and supernatants were saved. Input samples were stored (30 μL) at 4°C while the remainder of the sheared chromatin was combined with antibody-coupled protein A magnetic beads (NEB) and incubated at 4°C overnight with constant rotation. Mbd3 antibody (Bethyl Labs A302-529A, RRID:AB_1998976), Mbd2 antibody (Bethyl Labs A301-632A, RRID:AB_1211478), or Dnmt1 antibody (Sigma D4692, RRID:AB_262096) coupled protein A magnetic beads (NEB) were blocked with 5 mg/mL BSA overnight at 4°C, prior to incubation with sheared chromatin. Magnetic beads were washed, material was eluted at 65°C on a thermomixer, and the eluent was incubated at 65°C overnight to reverse crosslinking. Input DNA was diluted with 170 μL elution buffer (20 mM Tris-HCl pH 8.0, 100 mM NaCl, 20 mM EDTA, 1% SDS) and treated similarly. Samples were treated with RNaseA/T1 (Ambion, Carlsbad, CA, USA) followed by proteinase K (Ambion) and then PCI extracted. Ethanol precipitated ChIP-enriched DNA was then used for library construction or as a template for qPCR reactions.

FLAG-ChIP experiments were performed as previously described (Chen et al., 2013). Cells were fixed, washed with PBS, and pelleted. After resuspension and incubation with Lysis buffer 1 (50 mM HEPES KOH pH 7.6, 140 mM NaC, 1 mM EDTA, 10% (w/v) glycerol, 0.5% NP-40, 0.25% Triton X100), pellets were sheared in Lysis buffer 2 (10 mM Tris-HCl pH 8.0, 200 mM NaCl, 1 mM EDTA, 0.5 mM EGTA) for sonication in a Bioruptor. Sheared chromatin was incubated at 4°C overnight with protein G magnetic beads (NEB) pre-blocked and coupled with anti-Flag M2 (Sigma F1804, RRID:AB_262044). Magnetic beads were washed, material was eluted at 65°C on a thermomixer, and the eluent was incubated at 65°C overnight to reverse crosslinking. Input DNA was diluted with 170 μL elution buffer and treated similarly. Samples were treated with RNaseA/T1 (Ambion) followed by proteinase K (Ambion) and then PCI extracted. Ethanol precipitated, ChIP-enriched DNA was then used for library construction or as a template for qPCR reactions.

Library construction

Libraries of ChIP-enriched DNA were prepared as described previously (Chen et al., 2013). Samples were end-repaired, A-tailed, and adaptor-ligated with DNA purification over a column between each step. DNA was PCR amplified with KAPA HiFi polymerase using 16 cycles of PCR. Each library was size-selected on a 1% agarose gel, its concentration determined, and the integrity was confirmed by sequencing ~10 fragments from each library. Libraries were sequenced on an Illumina HiSeq2000 using single-end sequencing at the UMass Medical School deep sequencing core facility.

Data analysis

Single-end fastq reads were split by barcode adapter sequences, adapter sequences were removed, and reads were mapped to the mm9 genome using bowtie (RRID:SCR_005476), allowing up to three mismatches. Aligned reads were processed in HOMER (Heinz et al., 2010, RRID:SCR_010881). Genome browser tracks were generated from mapped reads using the ‘makeUCSCfile’ command. Peaks were called using the ‘findPeaks’ command after reads were mapped. Peaks were called individually from replicate datasets and those peaks enriched in both datasets were retained. Mapped reads were aligned over specific regions using the ‘annotatePeaks’ command to make 20 bp bins over regions of interest and sum the reads within each window. Experiments were aligned over the following datasets: TSS reference sites (from HOMER software), intermediate CpG promoters (ICPs, defined in (Weber et al., 2007), low methylated regions (LMRs from (Stadler et al., 2011) with or without additional selection criteria used by Baubec et al. – exclusion of LMRs under 150 bp and within 3 kb of other LMRs, referred to as 'LMR subset'), DNaseI hypersensitive sites (DHSs) from mouse ENCODE data (GSM1014154) with TSS locations removed, Mbd3 peaks (peaks called from both EGFP KD Mbd3 ChIP-seq libraries using the ‘findPeaks’ command in HOMER with input as a control) and Mbd2 peaks (peaks called from both EGFP KD Mbd2 ChIP-seq libraries using the ‘findPeaks’ command in HOMER with input as a control). After anchoring mapped reads over reference sites, aggregation plots were generated by averaging data obtained from biological replicates. Overlapping peaks for Figure 2B were identified using the ‘mergePeaks’ command.

Analysis of data from (Baubec et al., 2013) was performed similarly. Data were downloaded from GSE39610 (see Table 1 for SRA numbers) and converted to fastq files using SRA dump. Barcodes were trimmed where necessary, and reads were mapped to the mm9 genome using bowtie (RRID:SCR_005476), allowing up to three mismatches. Aligned reads were processed in HOMER (Heinz et al., 2010, RRID:SCR_010881) and analysis was confirmed independently using R and Macs2. Peaks for bio-Mbd2, bio-Mbd3, bio-Mbd1a, bio-Mbd1b, bio-Mbd4, or bio-MeCP2 were called using the ‘findPeaks’ command after reads were mapped and separately called using Macs2 to confirm validity of called peaks. Mapped reads were aligned over specific regions using the ‘annotatePeaks’ command to make 20 bp bins over regions of interest and sum the reads within each window. Experiments were aligned over the same genomic features described above and bio-Mbd3 peaks (called from bio-Mbd3 ChIP-seq in WT cells using ‘findPeaks’ command in HOMER with pooled input as a control), bio-Mbd2 peaks (called from bio-Mbd2 ChIP-seq in WT cells using ‘findPeaks’ command in HOMER with pooled input as a control), bio-Mbd1a peaks (called from bio-Mbd1a ChIP-seq in WT cells using ‘findPeaks’ command in HOMER with pooled input as a control), bio-Mbd1b peaks (called from bio-Mbd1b ChIP-seq in WT cells using ‘findPeaks’ command in HOMER with pooled input as a control), bio-Mbd4 peaks (called from bio-Mbd4 ChIP-seq in WT cells using ‘findPeaks’ command in HOMER with pooled input as a control), or bio-MeCP2 peaks (called from bio-MeCP2 ChIP-seq in WT cells using ‘findPeaks’ command in HOMER with pooled input as a control). To confirm this analysis, bio-Mbd2 and bio-Mbd3 data were downloaded and mapped independently and mapped reads were annotated and summarized using the Bioconductor package ChIPpeakAnno (Zhu et al., 2010; Zhu, 2013, RRID:SCR_012828). Furthermore, bio-Mbd2 and bio-Mbd3 peaks were confirmed independently by calling peaks from WT bio-Mbd2 or WT bio-Mbd3 ChIP-seq experiments, respectively, using Macs2. Finally, we re-analyzed a subset of these data after re-aligning reads (using bowtie, RRID:SCR_005476) while allowing only two mismatches instead of three. The use of two or three mismatches during alignment had no effect on the results.

Analysis of LMRs was performed as follows. LMR locations were downloaded from Stadler et al. (2011) (referred throughout as 'total LMRs'). As per the description in Baubec et al. (2013), LMRs < 150 bp and within 3 kb of other LMRs or unmethylated regions were removed (referred throughout as 'LMR subset'). The LMRs not included in this subset are referred throughout as 'removed LMRs'. To assess the genomic locations of these LMRs, the peaks were annotated in HOMER (Heinz et al., 2010, RRID:SCR_010881) using the ‘annotatePeaks’ command.

Mbd3 ChIP-seq data from (Yildirim et al., 2011) was downloaded and converted to fastq using SRA-dump from GSE31690. Reads were mapped to the mm9 genome using bowtie, allowing up to three mismatches. Aligned reads were processed in HOMER (Heinz et al., 2010). Mapped reads were aligned over gene-distal DHSs. Replicates for Mbd3 ChIP-seq in control cells were examined and resulted in similar profiles. Only replicate 2 (which has higher read coverage) is shown in Figure 2—figure supplement 4E.

Methylated and hydroxymethylated DNA immunoprecipitation sequencing (MeDIP-seq and hMeDIP-seq)

MeDIPs and hMeDIPs were performed as previously described (Mohn et al., 2009) with slight modifications. Genomic DNA was isolated from ES cells by resuspending cells in ES cell lysis buffer (10 mM Tris-HCl pH 7.5, 10 mM EDTA, 10 mM NaCL, 0.5% sarkosyl), treating with RNase A/T1 (Ambion), and incubating with 1 μg/μL proteinase K overnight at 55°C. DNA was cleaned through PCI extraction and isopropanol precipitation. Genomic DNA was then diluted to 20 ug/mL in TE buffer and sonicated in a Bioruptor UCD-200 (Diagenode) on low for 20 min with intervals of 15 s on/15 s off at 4°C. Sheared DNA (ranging in size from 200 to 600 bp) was prepared for deep sequencing and further processed by end-repairing, A-tailing, and adaptor-ligating as previously described for MNase-Seq libraries (Hainer et al., 2015) with DNA purification through PCI extraction and ethanol precipitation between each step. Of the prepared DNA, 1.5 μg was diluted in TE buffer up to 450 μL, denatured for 10 min at 95°C and immediately transferred to ice for 10 min. Of the prepared DNA, 1.5 μg was stored for input samples. Denatured DNA was then incubated with 50 μL 10X DNA IP buffer (100 mM Na-HPO4 pH 7.0, 1.4 M NaCl, 0.5% Triton X-100) and either 1.5 μL anti-5mC (Eurogentec Fremont, CA, USA, BI-MECY, RRID:AB_2616058) or 1.5 μL anti-5hmC (Active Motif, Carlsbad, CA, USA, 39791, RRID:AB_2630381) overnight at 4°C with constant rotation. 40 μL of magnetic beads (anti-mouse IgG for 5mC and anti-rabbit IgG for 5hmC, Life Technologies, Waltham, MA, USA) were pre-blocked with BSA and then resuspended in 1X DNA IP buffer. Beads were then added to the DNA/antibody mixture and incubated for 2 hr at 4°C with constant rotation. Beads were then washed three times with 700 μL 1X DNA IP buffer, and DNA was eluted from the beads by incubating at 50°C on a thermomixer with 250 μL 1X DNA IP buffer and 4 μL proteinase K (20 μg/μL) for 3 hr. Input samples were diluted up to 250 μL with DNA IP buffer and treated similarly. DNA from both input and IP samples were cleaned through PCI extraction and ethanol precipitation. DNA was PCR amplified with KAPA HiFi polymerase using 16 cycles of PCR. Each library was size-selected on a 1% agarose gel, its concentration was determined, and the integrity was confirmed by TOPO-cloning (Invitrogen) a portion and sequencing ~10 fragments from each library. Libraries were sequenced on an Illumina HiSeq2000 using single-end sequencing at the UMass Medical School deep sequencing core facility. Resulting fastq sequences were processed and analyzed as described for ChIP-seq libraries using the HOMER software (Heinz et al., 2010, RRID:SCR_010881).

Dot blotting

Twofold serial dilutions of sheared genomic DNA (200–400 bp) starting at 300 ng were denatured at 95°C for 10 min, then put on ice immediately for 10 min. Denatured DNA samples were spotted onto Amersham Hybond N+ nylon membrane (GE Healthcare, Uppsala, Sweden) and the membranes were UV crosslinked. For 5mC, the membrane was incubated in 0.1% SDS overnight, washed five times with PBS-T, blocked with 5% nonfat milk and 3% BSA for 4 hr, washed three times with PBS-T, incubated with anti-5mC (1:1000, Eurogentec BI-MECY, RRID:AB_2616058) for 1 hr, washed three times with PBS-T, incubated with HRP conjugated anti-mouse secondary (1:10,000, Bio-Rad, Hercules, CA, USA, 170–6516, RRID:AB_11125547) for 1 hr, washed three times with PBS-T, and detected with enhanced chemiluminescence. For 5hmC, the membrane was blocked with 5% nonfat milk for 4 hr, washed three times with PBS-T, incubated overnight with anti-5hmC (1:1000, Active Motif 39791, RRID:AB_2630381), washed three times with PBS-T, incubated with HRP conjugated anti-rabbit secondary (1:10,000, Bio-Rad 170–6515, RRID:AB_11125142) for 1 hr, washed three times with PBS-T, and detected with enhanced chemiluminescece. For loading, sheared gDNA samples were diluted simultaneously, spotted directly onto Amersham Hybond N+ nylon membrane (GE Healthcare) and the membranes were UV crosslinked. Membranes were incubated with 0.2% methylene blue for 5 min and washed five times with water.

Thin layer chromatography (TLC)

TLC was performed as previously described (Ficz et al., 2011; Tahiliani et al., 2009). 2 μg of genomic DNA isolated from control, Dnmt1 KD, Mbd3 KD, or Mbd2 KD ES cells was digested overnight with 50 units of MspI (NEB) and 10 μg of RNaseA/T1 (Ambion) at 37°C. The reaction was then heat inactivated at 65°C, DNA was dephosphorylated with 10 units of rSAP (NEB) for 1 hr at 37°C and cleaned over a DNA clean and concentrator column (Zymo Research, Irvine, CA, USA). 400 ng of DNA was incubated with 10 μCi 32P-ATP and 10 U of T4 PNK for 1 hr at 37°C. Labeled DNA was precipitated, resuspended in 18 μL 30 mM Tris pH 8.9, 15 mM MgCl2, 2 mM CaCl, and fragmented to single nucleotides with 10 units of DNaseI (NEB) and 10 μg SVPD (Worthington, Lakewood, NJ, USA) for 3 hr at 37°C. 3 μL of samples were spotted onto PEI cellulose F TLC plates (Millipore) and developed in isobutyric acid:H2O:NH3 (66:20:1 v/v/v) for 15 hr. Plates were dried and exposed to film for 30 hr. Plate was then exposed to a phosphorimager screen for 110 hr and scanned on a Typhoon FLA 700 (GE Healthcare). Pixels were quantitated using Fiji ImageJ software.

Methylation-sensitive restriction digests

Five micrograms of genomic DNA isolated from control, Dnmt1 KD, Mbd3 KD, or Mbd2 KD ES cells was digested at 37°C for 12 hr with 50 units of MspI or HpaII (both from NEB) followed by 20 min of heat inactivation at 65°C. MspI and HpaII recognize the 4 bp sequence CCGG, but HpaII is blocked by CpG methylation whereas MspI is insensitive to CpG methylation. 15 μL of the 50 μL reaction was run on a 1% agarose gel to visualize the cutting efficiency of each enzyme.

Microarray analysis

Microarray analysis was performed as previously described (Chen et al., 2013). Total RNA from control (EGFP KD), Dnmt1 KD, Tet1 KD, Mbd3 KD, Mbd2 KD, or Chd4 KD ES cells was subjected to RNA amplification and labeling using the Low Input Quick Amp Labeling Kit (Agilent, Santa Clara, CA, USA). Following this procedure, 50 picomoles of cRNA was used for fragmentation and hybridization on Agilent 4 × 44K mouse whole-genome microarrays. Slides were scanned on a DNA microarray scanner (Agilent G2565CA), and fluorescence data were obtained using Agilent Feature Extraction software. The expression data from biological replicates were fit to a linear model using the Bioconductor package limma (Smyth et al., 2005; 2003; Wettenhall and Smyth, 2004, RRID:SCR_010943), and analyzed as previously described (Yildirim et al., 2011). Datasets were visualized through Java TreeView (Saldanha, 2004, RRID:SCR_013503). Tip60 KD and Brg1 KD microarray data was downloaded from GSE31008 (Yildirim et al., 2011) and compared to averaged microarray data for fold change of Dnmt1 KD, Tet1 KD, Mbd3 KD, Mbd2 KD, and Chd4 KD. Correlation coefficients for all genes in pairwise factor KD were calculated in R.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
    De novo DNA cytosine methyltransferase activities in mouse embryonic stem cells
    1. H Lei
    2. SP Oh
    3. M Okano
    4. R Jüttermann
    5. KA Goss
    6. R Jaenisch
    7. E Li
    (1996)
    Development 122:3195–3205.
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70

Decision letter

  1. Jeannie T Lee
    Reviewing Editor; Massachusetts General Hospital, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

[Editors’ note: a previous version of this study was rejected after peer review, but the authors submitted for reconsideration and the paper was subsequently accepted for publication. The first decision letter after peer review is shown below.]

Thank you for submitting your work entitled "DNA methylation directs genomic localization of Mbd2 and Mbd3 in ES cells" for consideration by eLife. Your article has been favorably evaluated by Kevin Struhl (Senior Editor) and three reviewers, one of whom is a member of our Board of Reviewing Editors. Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews below, we regret to inform you that we cannot publish your manuscript in its present form.

All reviewers agree that your manuscript addresses an important area of research, that it adds important information to a high-profile controversy, and that the new experiments and re-analysis of older data are a valuable contribution to the field. As such, your study remains of high interest to eLife such that a suitably revised manuscript could be reconsidered. There are a number of important technical concerns that we think might take some additional time to address. One of the most important pertains to Figure 6A, as discussed by the first two reviewers. A more objective and quantitative analysis of the microarray expression profile is necessary to make the strong statements about high correlation. Additionally, we agreed that antibody specificity could be an issue and that this study would add positively to the debate if the authors could confirm their IP data in two ways. First, with FLAG ChIP of tagged MBDs to validate their MBD ChIPseq data. Second, with a parallel approach such as TLC to validate their MeDIP and hMeDIP analysis. Finally, all reviewers agreed that the prose could be toned down when comparing the Baubec et al. paper to their previous paper (Yildirim et al.).

Reviewer #1:

DNA methylation 'readers' include MBD2, which binds 5-methyl-cytosine, and MBD3, which binds its oxidized form, 5-hydroxymethyl-cytosine. Hainer et al. take a new look at previously published data by Baubec et al. (Cell 153:480, 2013) and perform new experiments, and they cast doubt on the likelihood of DNA methylation-independent binding of MBD2 and MBD3. Further, they find that, surprisingly, MBD2 binding depends on Tet1 gene function, and that MBD2 and MBD3 binding to DNA is mutually dependent.

Overall, this is an important area of research and an area of active debate. The new data provided by the authors, along with their reexamination of the existing data add significantly to the debate. On the other hand, the manuscript does not significantly advance understanding of mechanism. For publication in eLife, one might ask the authors to provide a better understanding of MBD function. These points may be open for discussion.

Specific comments in roughly descending order of importance:

Figure 6A. The strong assertion is made that gene expression profiles upon knockdown of Dnmt1, Tet1, Mbd3, Mbd2 and Chd4 are all "well-correlated." This seems to be based on visual inspection of the heat maps of microarray results, i.e., the top rows across the different treatments are mostly unregulated genes, and the bottom rows are mostly down regulated. But it seems there should be a more quantitative way to assess patterns of gene regulation. I'd like to see how the expression profiles of Mbd/Tet knockdowns would compare to perturbing a different repressive mechanism, such as knockdown of PRC1 or PRC2 components. Cluster analysis could be done of the columns (different treatments) of the heat map rather than just the rows (genes). You'd expect that, say, profiles of Eed1 and Ezh1 knockdowns would sort together as an "outgroup" compared to the Mbd-Tet-etc treatments.

Mutual dependence of Mbd2, Mbd3, Tet1, 5mC/5hmC: Figure 6B is a nicely done visualization of a possible regulatory cycle that explains the data well. I'd like to see more discussion however of why the different phenotypes vary so much in intact mice for Tet1, Dnmt1, Mbd2 and Mbd3 (e.g. the molecular phenotypes and regulatory cycle here may apply only to ES cells). Does the model predict that the MBDs recruit TET enzymes? Can they be Co-IP'd?

LMR terminology should be very explicit at its first mention, subsection “Unbiased analyses of published data reveal methylation-dependent localization of Mbd3 and Mbd2”, second paragraph. It is confusing later on: in the third paragraph of the aforementioned subsection, and subsection “Loss of Mbd3 or Mbd2 results in reduced 5mC and 5hmC at regulatory regions”. I believe that plain "LMR" is meant strictly in the sense of Baubec et al., 2013, who in turn (I think) use it the sense of Stadler et al., Nature 480:490, 2011, meaning regions >150bp with ~30% methylation. (Does "enrichment of TF binding sites" also apply here?) The authors contrast "LMR" with "total LMRs", which should be defined explicitly. (What is the minimal length? Does the methylations status differ from ~30% etc.)?

Figure 5B: light yellow trace for Mbd2 KD Mbd3 ChIP is not visible at all in first two panels ("ICPs" and "LMRs"), and it's hard to see in the other panels. (Presumably it's obscured by brown Mbd3 KD Mbd2 ChIP trace?)

Reviewer #2:

Proteins with methyl-binding domains (MBDs) constitute a major class of methylated DNA "readers". The manuscript by Hainer et al. addresses conflicting reports concerning function and binding specificity of two members of that protein family, Mbd2 and Mbd3. There is a disagreement between the findings published in 2011 by the same group and subsequent reports from the Schuebeler lab and Spruijt et al. on the role of CpG methylation and hydroxymethylation in the binding of these two proteins.

This is an intriguing area of research and a high-profile controversy, and the new experiments and re-analysis of the existing data presented by Hainer et al. are a valuable contribution to the field. The authors argue (convincingly, in my view) that many apparent discrepancies are due to differences in study design and analysis. These include particular genomic features analyzed in ChIP experiments, to whether binding specificity was measured using synthetic methylated DNA probes. Among the new experiments, the demonstration that the Tet1 catalytic activity was needed for Mbd2/3 occupancy (Figure 3) in particular served to support authors' argument.

Some of these differences seem highly intriguing but are not pursued any further. For instance, Figure 1 shows that Mbd2 binding in "total LMRs" shows much stronger dependence on methylation than in the subset of LMRs analyzed in Baubec et al. The apparent implication is that most of the methylation dependence is concentrated in those "other" regions which constitute the difference between those two sets. What is special about those regions?

In weighing strengths and weaknesses for eLife publication one has to balance the overall quality of the descriptive data presented here with the limited insights into mechanisms and function. I am on the fence as to whether this work presents a sufficiently major advance regarding processes involved in the Mbd2/3 function.

By far the weakest arguments of this paper are presented in Figure 6. Similarity of all expression profiles shown in Figure 6A is not self-evident. It looks like there are some downregulated genes in the top half of the heatmap and upregulated genes in the bottom half, but it's hard to see at the used resolution and color scheme. A more formal analysis of similarity would strengthen the case the authors make. In addition, a change in the heatmap color scheme might help underscore consistency of up- and down-regulated gene sets.

The "regulatory loop" scheme (Figure 6B) as formulated in Discussion is strongly overstated. In what sense "the outcome is the same" for disruption of any component of the proposed loop when the phenotypes of Mbd2, Mbd3, and Tet KO ESCs and mice are different? Accordingly, the statement (in the Abstract) that the authors' findings "describe a regulatory loop" is too strong.

Reviewer #3:

Summary:

In mammalian cells 5-methylcytosine (5mC) is a well-known heritable epigenetic modification, with well-studied effector proteins. 5-hydroxymethylcytosine (5hmC), initially identified as transient by product of active demethylation process, is less understood. Recent studies show high steady-state levels of 5hmC and effector proteins that preferentially bind to this modification and imply that 5hmC is an abundant form of cytosine modification in mammalian cells. The extent to which Mbd2 and Mbd3 interact with these modifications has been controversial.

In this study the authors investigate the interdependency of Mbd2 and Mbd3 for chromatin localization through binding to DNA methylation and hydroxymethylation sites respectively; they also study how the levels of these two Mbd proteins impact the stability of these two cytosine modifications. The requirement of Tet1 catalytic activity for the localization of Mbd3 and Mbd2 and effect of Mbd3 depletion on Dnmt1 localization are important and interesting findings. The paper offers an advance by providing a good argument that there is strong interdependence between these two binding proteins and the two covalent modifications that they recognize.

I found this paper interesting, but have a few points that are mainly related to helping to resolve the controversies in this area.

1) The interdependence of 5mC and 5hmC levels relative to Mbd2 and Mbd3 was demonstrated with MeDIP and hMeDIP experiments (Figure 4). Since these methods depend on antibody specificity, it would be useful to support these findings with another method such as TLC on overall levels of 5mC and 5hmC.

2) The interdependence of Mbd2 and Mbd3 binding was shown by ChIP with the antibodies that recognize endogenous proteins (Figure 5). Since Mbd2 and Mbd3 have 80% sequence similarity, it would be nice to see the genome-wide effect of Mbd2 and Mbd3 depletion on Mbd3-Flag and Mbd2-Flag occupancy with flag-tag ChIP. This would rule out cross-specificity artifacts.

3) I think the paper would read better if the authors toned down the prose when comparing the Baubec et al. paper to their previous paper (Yildirim et al.) The point they are making is clear, but perhaps too clear as it has a bit of a sledgehammer aspect to it.

4) Figure 5 shows a strong interdependency of Mbd2 and 3. How do the authors explain the phenotypic differences between Mbd3 and Mbd2 null mice?

https://doi.org/10.7554/eLife.21964.030

Author response

[Editors’ note: the author responses to the first round of peer review follow.]

[…]

Reviewer #1:

DNA methylation 'readers' include MBD2, which binds 5-methyl-cytosine, and MBD3, which binds its oxidized form, 5-hydroxymethyl-cytosine. Hainer et al. take a new look at previously published data by Baubec et al. (Cell 153:480, 2013) and perform new experiments, and they cast doubt on the likelihood of DNA methylation-independent binding of MBD2 and MBD3. Further, they find that, surprisingly, MBD2 binding depends on Tet1 gene function, and that MBD2 and MBD3 binding to DNA is mutually dependent.

Overall, this is an important area of research and an area of active debate. The new data provided by the authors, along with their reexamination of the existing data add significantly to the debate. On the other hand, the manuscript does not significantly advance understanding of mechanism. For publication in eLife, one might ask the authors to provide a better understanding of MBD function. These points may be open for discussion.

Specific comments in roughly descending order of importance:

Figure 6A. The strong assertion is made that gene expression profiles upon knockdown of Dnmt1, Tet1, Mbd3, Mbd2 and Chd4 are all "well-correlated." This seems to be based on visual inspection of the heat maps of microarray results, i.e., the top rows across the different treatments are mostly unregulated genes, and the bottom rows are mostly down regulated. But it seems there should be a more quantitative way to assess patterns of gene regulation. I'd like to see how the expression profiles of Mbd/Tet knockdowns would compare to perturbing a different repressive mechanism, such as knockdown of PRC1 or PRC2 components. Cluster analysis could be done of the columns (different treatments) of the heat map rather than just the rows (genes). You'd expect that, say, profiles of Eed1 and Ezh1 knockdowns would sort together as an "outgroup" compared to the Mbd-Tet-etc treatments.

As suggested, in the revised manuscript we quantitatively compare the gene expression profiles of the various KDs and make use of KDs of unrelated factors as outgroups. These analyses show the degree to which the MBD, DNMT, and TET KDs have highly overlapping effects on gene expression, and verify that these alterations are much more similar to each other than the outgroups. The finding that NuRD mediated deacetylation of H3K27 is necessary for PRC2- mediated H3K27 methylation at many loci (Reynolds et al., 2012) complicates use of PRC2 KDs as outgroups. As alternatives, we compared our data with Tip60 KD, representing a chromatin remodeling complex known to mainly repress gene expression in ES cells (Fazzio et al., 2008), and Brg1 KD, representing a (mostly) activating chromatin remodeling complex. We downloaded previously published microarray data for both Tip60 KD and Brg1 KD (Yildirim et al., 2011) and included these analyses in our revised manuscript (Figure 6—figure supplement 1A). This new analysis reveals that gene expression changes in Tip60 KD or Brg1 KD ES cells are dissimilar to those of Dnmt1, Tet1, Mbd3, Mbd2 or Chd4 KDs. We further show individual KDs of Dnmt1, Tet1, Mbd3, Mbd2, and Chd4 correlate well with each other, but correlate poorly with Tip60 or Brg1 KDs (Figure 6—figure supplement 1B).

Mutual dependence of Mbd2, Mbd3, Tet1, 5mC/5hmC: Figure 6B is a nicely done visualization of a possible regulatory cycle that explains the data well. I'd like to see more discussion however of why the different phenotypes vary so much in intact mice for Tet1, Dnmt1, Mbd2 and Mbd3 (e.g. the molecular phenotypes and regulatory cycle here may apply only to ES cells). Does the model predict that the MBDs recruit TET enzymes? Can they be Co-IP'd?

Previously, we identified a weak interaction between Mbd3/NuRD and Tet1 (Yildirim et al. 2011), which was later confirmed by another group (Shi, Kim et al., 2013). We have attempted co-IP experiments between Mbd2 and Tet1 but have not detected an interaction. We have added discussion of these findings in reference to our model in the fourth paragraph of the Discussion and modified our model to remove the arrows between Mbd2/3 and Tet1 to avoid conflating hypothetical functional connections (that are unnecessary for the model) with connections that are established from data in this manuscript and the published literature. Furthermore, we added more discussion regarding the phenotypes observed in mice and ES cells and their implications with regards to our findings (Discussion, last paragraph).

LMR terminology should be very explicit at its first mention, subsection “Unbiased analyses of published data reveal methylation-dependent localization of Mbd3 and Mbd2”, second paragraph. It is confusing later on: in the third paragraph of the aforementioned subsection, and subsection “Loss of Mbd3 or Mbd2 results in reduced 5mC and 5hmC at regulatory regions”. I believe that plain "LMR" is meant strictly in the sense of Baubec at et., 2013, who in turn (I think) use it the sense of Stadler et al., Nature 480:490, 2011, meaning regions >150bp with ~30% methylation. (Does "enrichment of TF binding sites" also apply here?) The authors contrast "LMR" with "total LMRs", which should be defined explicitly. (What is the minimal length? Does the methylations status differ from ~30% etc.)?

In the revised manuscript, we have attempted to clarify the difference by referring to the subset of LMRs examined by Baubec et al. as the “LMR subset” throughout the text. We utilized the same criteria defined by Baubec et al. to generate the LMR subset (remove LMRs smaller than 150bp and those within 3kb of another LMR or a UMR) from the total set of LMRs defined in Stadler et al., 2011 (which we now refer to as “total LMRs”). By the definition of LMRs in Stadler et al., total LMRs have on average 30% methylation, although they range in methylation from 10- 50%. We found that 72% of the LMRs excluded by Baubec et al. (i.e., present within total LMRs but not the LMR subset) were excluded because they are less than 150bp in size, and are therefore smaller than those included in the LMR subset. To address whether these “removed LMRs” are enriched for TF binding sites, we examined the frequencies of 319 known TF motifs in each of the LMR subset, removed LMRs, or total LMR locations. We found on average that the LMR subset had 12.5 motifs per LMR, removed LMRs had 13, and total LMRs had 12.7, suggesting that all three classes of LMRs have similar representation of TF binding sites. Closer examination of the differences between total LMRs and the LMR subset showed that removed LMRs were more highly methylated and hydroxymethylated (compare levels of 5mC and 5hmC in EGFP KD for removed LMRs in Figure 4—figure supplement 1C-D to total LMRs shown in Figure 4—figure supplement 1A-B and the LMR subset shown in Figure 4A-B; we observe a 5mC peak height of 1.48, 0.96, and 2.64 for total, subset, and removed LMRs, respectively and a 5hmC peak height of 1.51, 1.16, and 2.25 for total, subset, and removed LMRs, respectively). In addition, removed LMRs were more enriched for promoter proximal locations (Figure 4—figure supplement 1E in revised manuscript).

Figure 5B: light yellow trace for Mbd2 KD Mbd3 ChIP is not visible at all in first two panels ("ICPs" and "LMRs"), and it's hard to see in the other panels. (Presumably it's obscured by brown Mbd3 KD Mbd2 ChIP trace?)

We thank the reviewer for pointing this out. We changed the thickness of the yellow trace to increase its visibility.

Reviewer #2:

Proteins with methyl-binding domains (MBDs) constitute a major class of methylated DNA "readers". The manuscript by Hainer et al. addresses conflicting reports concerning function and binding specificity of two members of that protein family, Mbd2 and Mbd3. There is a disagreement between the findings published in 2011 by the same group and subsequent reports from the Schuebeler lab and Spruijt et al. on the role of CpG methylation and hydroxymethylation in the binding of these two proteins.

This is an intriguing area of research and a high-profile controversy, and the new experiments and re-analysis of the existing data presented by Hainer et al. are a valuable contribution to the field. The authors argue (convincingly, in my view) that many apparent discrepancies are due to differences in study design and analysis. These include particular genomic features analyzed in ChIP experiments, to whether binding specificity was measured using synthetic methylated DNA probes. Among the new experiments, the demonstration that the Tet1 catalytic activity was needed for Mbd2/3 occupancy (Figure 3) in particular served to support authors' argument.

Some of these differences seem highly intriguing but are not pursued any further. For instance, Figure 1 shows that Mbd2 binding in "total LMRs" shows much stronger dependence on methylation than in the subset of LMRs analyzed in Baubec et al. The apparent implication is that most of the methylation dependence is concentrated in those "other" regions which constitute the difference between those two sets. What is special about those regions?

As noted in the response to reviewer #1, most of the LMRs excluded by Baubec et al. were excluded due to their size (<150bp), while a smaller portion were excluded due to their location within 3kb of another LMR or a UMR. After examining the differences between total LMRs and the LMR subset taken by Baubec et al., we found that the LMRs that were not included in the LMR subset (“removed LMRs”) were enriched for promoter proximal locations, and had higher overall methylation and hydroxymethylation than the LMR subset analyzed in Baubec et al. (Figure 4—figure supplement 1C-E). The higher 5mC and 5hmC levels in the removed LMRs can be observed by comparing the peak heights of MeDIP and hMeDIP in EGFP KD cells over total LMRs and the LMR subset (Figure 4—figure Supplement 1A-B and Figure 4A-B) to EGFP KD cells at removed LMRs (Figure 4—figure supplement 1C-D).

In weighing strengths and weaknesses for eLife publication one has to balance the overall quality of the descriptive data presented here with the limited insights into mechanisms and function. I am on the fence as to whether this work presents a sufficiently major advance regarding processes involved in the Mbd2/3 function.

We feel that our findings represent a significant advance in this field. Not only do these data and analyses clarify a major error in the published literature, they also reveal novel findings regarding the interdependence between Mbd2, Mbd3, and the DNA modification machinery. Although previous studies in human cells have shown some overlap in Mbd2 and Mbd3 binding at specific genes, the interdependence between these factors has remained unexplored until now. Along with finding that Mbd2 and Mbd3 are necessary for normal chromatin localization by the other, we also found that both factors are necessary for normal levels of DNA methylation and hydroxymethylation – a very surprising finding in our opinion, and a finding that will have to be considered when mutants/KDs of MBD proteins are used in future studies. Therefore, although the first half of our manuscript is devoted to correcting an error in the literature, we feel that the novel insights into the interdependence between these marks and their “readers”, presented in the latter half of the manuscript, represent important new findings.

By far the weakest arguments of this paper are presented in Figure 6. Similarity of all expression profiles shown in Figure 6A is not self-evident. It looks like there are some downregulated genes in the top half of the heatmap and upregulated genes in the bottom half, but it's hard to see at the used resolution and color scheme. A more formal analysis of similarity would strengthen the case the authors make. In addition, a change in the heatmap color scheme might help underscore consistency of up- and down-regulated gene sets.

As suggested, we quantitatively assessed the alterations in gene expression in these KDs and brightened the colors in the heatmap to aid visibility. Our analyses of these data now include quantification of the degree to which the gene expression changes in each KD correlate, and we now include data from two unrelated KDs as outgroups for comparison (see response to the first point of reviewer #1 for rationale). These analyses are shown in Figure 6—figure supplement 1, and reveal significantly correlated changes among Dnmt1, Tet1, Mbd3, Mbd2, and Chd4 KDs.

Regarding the colors of the heatmap, we wanted to maintain a color scheme that is visible to individuals with red-green colorblindness, and have found the yellow-blue scheme is among the best in this regard, but we brightened the blue and yellow colors to increase visibility.

The "regulatory loop" scheme (Figure 6B) as formulated in Discussion is strongly overstated. In what sense "the outcome is the same" for disruption of any component of the proposed loop when the phenotypes of Mbd2, Mbd3, and Tet KO ESCs and mice are different? Accordingly, the statement (in the Abstract) that the authors' findings "describe a regulatory loop" is too strong.

In the revised manuscript, we toned down this statement, altering the last sentence in the Abstract along with the description of the model in the Discussion. In particular, we now discuss the functional overlap of these factors specifically in ES cell gene regulation, and acknowledge that this functional overlap is likely to change during differentiation and development (when 5hmC levels are reduced), which we do not address in this study. In addition, we have expanded our discussion regarding the phenotypes of Mbd2 and Mbd3 KO mice (please also see our response to the last comment from reviewer #3 and the revised text in the last two paragraphs of the Discussion). Unfortunately, there is some controversy regarding Tet phenotypes in the field, and since our data do not directly address this controversy, we discuss these issues only briefly in the Introduction. To clarify, phenotypic alterations upon KD or KO of Tet genes have been described in several studies. In two studies, KD of Tet1 impaired ES cell self-renewal, in part because of Nanog down-regulation (Ito et al., 2010, Freudenberg et al., 2012). However, separate studies showed no effect of Tet1 KO on self-renewal (Dawlaty et al., 2011; Koh et al., 2011; Ficz et al., 2011). On the other hand, Tet1 single KO and Tet1 Tet2 double KO ESCs are reported to exhibit a skewed differentiation profile (Dawlaty et al., 2011, Dawlaty et al., 2013), although Tet1 Tet2 double KO mice develop to term and have only modest phenotypes (Dawlaty et al., 2013). Because of these somewhat confusing results and the fact that our ES cell data do not address any of these differentiation or developmental phenotypes, we have confined our conclusions in the revised manuscript to the roles of these factors in gene regulation in ES cells.

Reviewer #3: I found this paper interesting, but have a few points that are mainly related to helping to resolve the controversies in this area.

1) The interdependence of 5mC and 5hmC levels relative to Mbd2 and Mbd3 was demonstrated with MeDIP and hMeDIP experiments (Figure 4). Since these methods depend on antibody specificity, it would be useful to support these findings with another method such as TLC on overall levels of 5mC and 5hmC.

As suggested, we performed TLC (Figure 4E) to analyze total levels and 5mC and 5hmC without the use of an antibody, as well as digests of genomic DNA with methylation sensitive or insensitive restriction enzymes (Figure 4F) to look at global cytosine methylation/hydroxymethylation levels. These new data, along with dot blot assays examining the total levels of 5mC and 5hmC (Figure 4C-D), are included in Figure 4 of the revised manuscript. The data from the antibody independent assays are consistent with the MeDIP-seq, hMeDIP-seq, and dot blots, all of which show a decrease in bulk 5mC and 5hmC levels upon depletion of Mbd2 or Mbd3.

2) The interdependence of Mbd2 and Mbd3 binding was shown by ChIP with the antibodies that recognize endogenous proteins (Figure 5). Since Mbd2 and Mbd3 have 80% sequence similarity, it would be nice to see the genome-wide effect of Mbd2 and Mbd3 depletion on Mbd3-Flag and Mbd2-Flag occupancy with flag-tag ChIP. This would rule out cross-specificity artifacts.

In the revised manuscript, we now include ChIP-seq data for Mbd3abc-3XFLAG and Mbd2- 3XFLAG in cells depleted of either Mbd3 or Mbd2, as requested (Figure 2—figure supplement 4 in revised manuscript). These data are consistent with our ChIP-qPCR data for Mbd3/Mbd3abc- 3XFLAG or Mbd2/Mbd2-3XFLAG in cells depleted of Mbd3 or Mbd2, respectively (Figure 2—figure supplement 2E; Figure 2—figure supplement 3C-D). Together, these data confirm the interdependence of Mbd2 and Mbd3 as well as the specificity of the antibodies against the endogenous proteins.

3) I think the paper would read better if the authors toned down the prose when comparing the Baubec et al. paper to their previous paper (Yildirim et al.) The point they are making is clear, but perhaps too clear as it has a bit of a sledgehammer aspect to it.

As suggested, we have toned down the manuscript to more delicately explain the differences between our data and those of Baubec et al., as well as the possible reasons underlying these differences.

4) Figure 5 shows a strong interdependency of Mbd2 and 3. How do the authors explain the phenotypic differences between Mbd3 and Mbd2 null mice?

While chromatin binding and gene regulation by Mbd2 and Mbd3 is largely interdependent in ES cells, this interdependence is not absolute – some Mbd2 and Mbd3 binding remains upon depletion of the other (Figure 5), and their gene expression profiles do not overlap 100% (Figure 6). In addition, ES cells are a transient cell type and we do not know whether the interdependency observed in ES cells is found in various differentiated cell types. Because the functional overlap is not absolute and may not extend to all cell types, this may explain their different phenotypes during development. We have added more discussion to the manuscript regarding the differences between these phenotypes, and potential explanations for these differences (Discussion, last two paragraphs).

https://doi.org/10.7554/eLife.21964.031

Article and author information

Author details

  1. Sarah J Hainer

    Department of Molecular, Cell, and Cancer Biology, University of Massachusetts Medical School, Worcester, United States
    Contribution
    SJH, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  2. Kurtis N McCannell

    Department of Molecular, Cell, and Cancer Biology, University of Massachusetts Medical School, Worcester, United States
    Contribution
    KNM, Acquisition of data, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.
  3. Jun Yu

    Department of Molecular, Cell, and Cancer Biology, University of Massachusetts Medical School, Worcester, United States
    Contribution
    JY, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  4. Ly-Sha Ee

    Department of Molecular, Cell, and Cancer Biology, University of Massachusetts Medical School, Worcester, United States
    Contribution
    L-SE, Drafting or revising the article, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.
  5. Lihua J Zhu

    1. Department of Molecular, Cell, and Cancer Biology, University of Massachusetts Medical School, Worcester, United States
    2. Program in Bioinformatics and Integrative Biology, University of Massachusetts Medical School, Worcester, United States
    3. Program in Molecular Medicine, University of Massachusetts Medical School, Worcester, United States
    Contribution
    LJZ, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  6. Oliver J Rando

    Department of Biochemistry and Molecular Pharmacology, University of Massachusetts Medical School, Worcester, United States
    Contribution
    OJR, Conception and design, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  7. Thomas G Fazzio

    1. Department of Molecular, Cell, and Cancer Biology, University of Massachusetts Medical School, Worcester, United States
    2. Program in Molecular Medicine, University of Massachusetts Medical School, Worcester, United States
    Contribution
    TGF, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article
    For correspondence
    thomas.fazzio@umassmed.edu
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0353-7466

Funding

Eunice Kennedy Shriver National Institute of Child Health and Human Development (HD072122)

  • Thomas G Fazzio

National Cancer Institute (T32CA130807)

  • Sarah J Hainer

Leukemia and Lymphoma Society (Special Fellows Award)

  • Sarah J Hainer

Eunice Kennedy Shriver National Institute of Child Health and Human Development (HD080224)

  • Oliver J Rando

Leukemia and Lymphoma Society (Scholars Award)

  • Thomas G Fazzio

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

Acknowledgements

We thank I Bach, E Torres, J Benanti, and members of the Fazzio lab for insightful discussions regarding this manuscript. This work was supported by the T32CA130807 training grant and a Leukemia and Lymphoma postdoctoral fellowship to SJH, NIH HD072122 to TGF, and NIH HD080224 to OJR TGF is a Leukemia and Lymphoma Society Scholar. SJH is currently a Special Fellow of the Leukemia and Lymphoma Society. The funders had no role in study design, data collection, analysis, decision to publish, or preparation of the manuscript. All deep sequencing was performed at the UMass Medical School Core facility on a HiSeq2000 supported by 1S10RR027052-01. Genomic datasets have been deposited within GEO (accession number GSE79771). The authors declare no conflict of interest.

Reviewing Editor

  1. Jeannie T Lee, Massachusetts General Hospital, United States

Publication history

  1. Received: September 29, 2016
  2. Accepted: November 1, 2016
  3. Version of Record published: November 16, 2016 (version 1)

Copyright

© 2016, Hainer et al.

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

Metrics

  • 1,650
    Page views
  • 344
    Downloads
  • 5
    Citations

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

Download links

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

Downloads (link to download the article as PDF)

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

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

Further reading

    1. Developmental Biology
    2. Stem Cells and Regenerative Medicine
    Anna Yoney et al.
    Research Article
    1. Developmental Biology
    2. Stem Cells and Regenerative Medicine
    Ashley RG Libby et al.
    Research Article