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
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.
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.
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).
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.
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).
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.
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.
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).
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.
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.
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.
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.
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.
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).
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-puroR (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.
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.
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).
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.
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.
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.
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).
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.
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.
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 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.
Transcriptional network analysis in muscle reveals AP-1 as a partner of PGC-1α in the regulation of the hypoxic gene programMolecular and Cellular Biology 34:2996–3012.https://doi.org/10.1128/MCB.01710-13
Probing the dynamic distribution of bound states for methylcytosine-binding domains on DNAJournal of Biological Chemistry 289:1294–1302.https://doi.org/10.1074/jbc.M113.512236
A family of human zinc finger proteins that bind methylated DNA and repress transcriptionMolecular and Cellular Biology 26:169–181.https://doi.org/10.1128/MCB.26.1.169-181.2006
Suppression of pervasive noncoding transcription in embryonic stem cells by esBAFGenes & Development 29:362–378.https://doi.org/10.1101/gad.253534.114
Identification and characterization of a family of mammalian methyl-CpG binding proteinsMolecular and Cellular Biology 18:6538–6547.https://doi.org/10.1128/MCB.18.11.6538
TET1 is a maintenance DNA demethylase that prevents methylation spreading in differentiated cellsNucleic Acids Research 42:6956–6971.https://doi.org/10.1093/nar/gku372
The NuRD component Mbd3 is required for pluripotency of embryonic stem cellsNature Cell Biology 8:285–292.https://doi.org/10.1038/ncb1372
Stalling of human DNA (cytosine-5) methyltransferase at single-strand conformers from a site of dynamic mutationJournal of Molecular Biology 275:67–79.https://doi.org/10.1006/jmbi.1997.1430
Genomic DNA methylation: the mark and its mediatorsTrends in Biochemical Sciences 31:89–97.https://doi.org/10.1016/j.tibs.2005.12.008
MBD2/NuRD and MBD3/NuRD, two distinct complexes with different biochemical and functional propertiesMolecular and Cellular Biology 26:843–851.https://doi.org/10.1128/MCB.26.3.843-851.2006
MBD3, a component of the NuRD complex, facilitates chromatin alteration and deposition of epigenetic marksMolecular and Cellular Biology 28:5912–5923.https://doi.org/10.1128/MCB.00467-08
Dissection of the methyl-CpG binding domain from the chromosomal protein MeCP2Nucleic Acids Research 21:4886–4892.https://doi.org/10.1093/nar/21.21.4886
Recombinant human DNA (cytosine-5) methyltransferase. I. Expression, purification, and comparison of de novo and maintenance methylationJournal of Biological Chemistry 274:33002–33010.https://doi.org/10.1074/jbc.274.46.33002
The p120 catenin partner Kaiso is a DNA methylation-dependent transcriptional repressorGenes & Development 15:1613–1618.https://doi.org/10.1101/gad.198501
Java Treeview--extensible visualization of microarray dataBioinformatics 20:3246–3248.https://doi.org/10.1093/bioinformatics/bth349
Transcriptional enhancers: from properties to genome-wide predictionsNature Reviews Genetics 15:272–286.https://doi.org/10.1038/nrg3682
Statistical issues in cDNA microarray data analysisMethods in Molecular Biology 224:111–136.https://doi.org/10.1385/1-59259-364-X:111
Integrative analysis of ChIP-chip and ChIP-seq datasetMethods in Molecular Biology 1067:105–124.https://doi.org/10.1007/978-1-62703-607-8_8
Jeannie T LeeReviewing 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.).
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?)
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.
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
- Thomas G Fazzio
- Sarah J Hainer
- Sarah J Hainer
- Oliver J Rando
- Thomas G Fazzio
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
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.
- Jeannie T Lee, Reviewing Editor, Massachusetts General Hospital, United States
- Received: September 29, 2016
- Accepted: November 1, 2016
- Version of Record published: November 16, 2016 (version 1)
© 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.