Coordinate H3K9 and DNA methylation silencing of ZNFs in toxicant-induced malignant transformation.

Genome-wide disruption of the epigenetic code is a hallmark of malignancy that encompasses many distinct, highly interactive modifications. Delineating the aberrant epigenome produced during toxicant-mediated malignant transformation will help identify the underlying epigenetic drivers of environmental toxicant-induced carcinogenesis. Gene promoter DNA methylation and gene expression profiling of arsenite-transformed prostate epithelial cells showed a negative correlation between gene expression changes and DNA methylation changes; however, less than 10% of the genes with increased promoter methylation were downregulated. Studies described herein confirm that a majority of the DNA hypermethylation events occur at H3K27me3 marked genes that were already transcriptionally repressed. In contrast to aberrant DNA methylation targeting H3K27me3 pre-marked silent genes, we found that actively expressed C2H2 zinc finger genes (ZNFs) marked with H3K9me3 on their 3' ends, were the favored targets of DNA methylation linked gene silencing. DNA methylation coupled, H3K9me3 mediated gene silencing of ZNF genes was widespread, occurring at individual ZNF genes on multiple chromosomes and across ZNF gene family clusters. At ZNF gene promoters, H3K9me3 and DNA hypermethylation replaced H3K4me3, resulting in a widespread downregulation of ZNF gene expression, which accounted for 8% of all the downregulated genes in the arsenical-transformed cells. In summary, these studies associate toxicant exposure with widespread silencing of ZNF genes by DNA hypermethylation-linked H3K9me3 spreading, further implicating epigenetic dysfunction as a driver of toxicant associated carcinogenesis.

Introduction

One important function of epigenetics is to provide a mechanism of genomic plasticity so that cells may respond to a changing and potentially harmful environment. Global disruption of the epigenetic code is a hallmark of malignancy that encompasses many distinct modifications that are highly interactive. For instance, during normal differentiation, histone H3 lysine 9 trimethylation (H3K9me3) acts as a substrate for heterochromatin binding proteins that recruit DNA methyltransferases thereby directing DNA methylation to specific genomic elements., Another repressive histone modification, H3 lysine 27 trimethylation (H3K27me3), predetermines many of the aberrant de novo DNA methylation sites in cancer most likely through protein-protein interactions which recruit DNA methyltransferases., On the other hand, trimethylation of histone H3 lysine 4 (H3K4me3) is found at the promoters of actively transcribed genes where it aids in the formation of an active chromatin conformation and prevents DNA methylation. Although some may function as individual enzymes, most epigenetic modifying proteins work as members of higher order protein complexes which often contain several epigenetic modifying enzymes, collectively controlling multiple distinct epigenetic marks. Generally, when one epigenetic mark is pathologically disrupted, several others are likely to be altered as a result. By characterizing the interplay between epigenetic modifications in toxicant-associated malignant transformation we can begin to uncover the underlying pathways frequently disrupted in cancer that drive widespread epigenetic malfunction.

The C2H2 zinc finger genes (ZNFs) make up one of the largest families of sequence specific DNA binding transcriptional repressors. At least a third of ZNF genes contain the KRAB domain (KRAB-ZNFs), which is known to interact with TRIM28 (also known as KAP-1) and indirectly recruit histone and DNA methyltransferases to specific genomic sequences that become marked with H3K9me3 and DNA methylation as a result.- In non-malignant cells, this targeted repressive activity of ZNFs participates in important anti-cancer defenses including regulation of cell proliferation and silencing of endogenous retroviruses.,, Because ZNFs are indirectly responsible for recruitment of distinct epigenetic modifying enzymes they represent a potential common pathway that when disrupted would alter multiple epigenetic marks at specific genomic loci.

Our past studies have identified global DNA methylation and histone modification changes that occurred in toxicant-treated, malignant cell lines relative to their untreated, immortal counterparts.- These models suggested that non-mutagenic environmental toxicants, such as arsenic, could induce aberrant agglomerative DNA methylation. We also found that in all the malignant cell lines and tumor biopsies tested, DNA hypermethylation occurred disproportionately within human embryonic stem cell (hESC) H3K27me3 and H3K9me3 domains with a striking correlation between densely packed hESC H3K9me3 domains and aberrantly methylated gene family clusters. Finding that DNA hypermethylation occurred disproportionately within H3K27me3 and H3K9me3 domains was not surprising, since these histone modifications are known to be linked to DNA methylation in certain scenarios., In the context of gene promoters, H3K27me3, H3K9me3 and DNA methylation are each known to participate in transcriptional repression, which led us to investigate the relationships between H3K27me3, H3K9me3, aberrant DNA hypermethylation and gene expression changes in environmental toxicant-associated malignant transformation.

In this pursuit, we have expanded on earlier studies by incorporating gene expression and selected chromatin immunoprecipitation analyses of the immortal prostate epithelial cell line RWPE-1 and its arsenic-exposed malignant variant, CAsE-PE, to better understand the relationships between aberrant DNA methylation, H3K9me3, H3K27me3 and gene expression of individual genes and gene family clusters in toxicant-associated malignant transformation. Our results confirm that H3K27me3 marked genes are targeted for aberrant DNA methylation and remain inactive. Contrarily, in immortalized cells we found that ZNFs in particular are frequently marked with H3K9me3 on their 3′ ends and H3K4me3 at the 5′ ends while actively transcribed. However, after malignant transformation, H3K9me3 had replaced H3K4me3 at the promoters of ZNFs and was accompanied by aberrant DNA hypermethylation, resulting in a widespread downregulation of ZNF gene expression.

Results

To interpret global DNA methylation changes associated with malignant transformation we performed MeDIP-on-Chip on 2 replicates of CAsE-PE and 2 replicates of RWPE-1 using Agilent Human promoter microarrays. This microarray probes the genomic regions between 5.5 kb upstream and 2.5 kb downstream of the transcription start sites (TSSs) of approximately 17,000 genes (Fig. S1A). Therefore, the probes are well positioned to detect DNA methylation changes that occur in close proximity to TSSs of protein coding genes but any genomic elements which are far from gene promoters have minimal coverage. With this probe distribution, the microarray is more likely to detect hypermethylation events associated with malignant transformation, since the DNA of normal cells is generally unmethylated near TSSs. In agreement with the array distribution and the general profile of DNA methylation at gene promoters in cancer, our final MeDIP-on-Chip data set of CAsE-PE relative to RWPE-1 contained predominately DNA hypermethylation events occurring near TSSs, while DNA hypomethylation events were fewer and more evenly distributed through the regions analyzed (Fig. S1A).

Hypermethylated genes share common ontologies

Overall, the MeDIP-on-Chip detected 869 hypermethylated and 166 hypomethylated genes in CAsE-PE cells (Supplemental Data Set 1). Most of the genes hypermethylated in CAsE-PE are H3K27me3 marked in human embryonic stem cells (hESC), consistent with data published previously (Fig. S1B). Several of the aberrant DNA methylation events detected by MeDIP-on-Chip were confirmed by MassARRAY technology (Fig. S2). To determine if particular classes of genes were targeted for hypermethylation in CAsE-PE, we performed a gene set enrichment analysis. Results showed that genes involved in developmental processes, and genes that are DNA or chromatin binding transcription factors were targeted for DNA hypermethylation (Table S1). This term enrichment could mean that these genes were selected for hypermethylation by nature of their function, which when silenced by methylation conferred some growth advantage to the cells. When taking into account that H3K27me3 domains of stem cells are targeted for aberrant DNA methylation in cancer, a more likely explanation is that these functional categories of genes are aberrantly methylated because they are normally regulated by the same epigenetic mechanisms, including H3K27me3 (Fig. S1B).

Widespread downregulation of ZNFs in CAsE-PE

In addition to DNA methylation profiling, we also performed gene expression profiling to identify changes that were associated with the malignant phenotype of CAsE-PE. We measured global gene expression levels in 3 replicates of CAsE-PE and 3 replicates of RWPE-1 using the Affymetrix Human Gene ST 1.0 microarray platform. Based on the criteria of fold change ≥ 2 and adjusted p-value ≤ 0.05, we found 428 downregulated and 514 upregulated genes in the malignant CAsE-PE cells relative to the non-malignant RWPE-1 (Supplemental Data Set 2). We found several gene expression changes that have been linked to cancer progression of the prostate and other tissues and we used QRT-PCR to confirm a small representative set (Fig. S3). Surprisingly, there was a highly significant downregulation of ZNF and KRAB-ZNF genes, comprising almost 8% of all the downregulated genes. Thirty-four out of the total 519 ZNF genes covered by the expression array were significantly downregulated and 26 of these were KRAB-ZNF genes (p < 1 × 10−12, hypergeometric test of ZNFs). Over-representation of the ZNFs and particularly the KRAB-ZNFs among the downregulated genes suggests that they were silenced in a coordinated fashion which may be driven by a common epigenetic mechanism and that their broad silencing could have a substantial impact on the phenotype of CAsE-PE cells.

A negative correlation between DNA methylation and gene expression is driven by hypermethylation of ZNF gene promoters

Cancer-associated, DNA hypermethylation occurs disproportionately in H3K27me3 and H3K9me3 domains, leading us to question how each of these potential targeting mechanisms would affect cellular phenotype since both these histone modifications and DNA methylation at gene promoters repress transcription. To begin to address this question we determined to what degree the DNA methylation changes correlated with gene expression changes. We merged the MeDIP-on-Chip data with the expression data by calculating the methylation changes within 2 kb windows centered on TSSs. Plotting the DNA methylation changes relative to the gene expression changes showed that they were negatively correlated which indicates that the hypermethylation did not only occur at genes that were already silenced by H3K27me3 or H3K9me3 (Fig. 1). Forty-six genes were significantly downregulated and had a log2 fold DNA methylation increase greater than 0.5 (Fig. 1, blue text). Some of the genes that met these criteria are LAMB2, ANXA6, SLC4A11, UCHL1, FGF5, PTGS2, SPARC, and BDNF. Upregulation of only three genes was linked to decreased methylation of the TSS (S100P, INSL6, and ADRB2, Fig. 1). In all, about 9% of the observed gene expression changes were linked to DNA methylation changes, with almost all the changes being repressive.

Figure 1. Expression by methylation plot. The blue-shaded background density cloud shows the distribution of all the genes covered by both arrays. Differentially expressed genes (553 genes, represented as black dots) were used to calculate a correlation between methylation change and expression change. Spearman rank correlation coefficient = -0.37. Genes labeled with blue text had a significant decrease in expression and a log2 increase in methylation greater than 0.5. Not all the genes that meet these criteria are labeled due to space limitations. Selected genes labeled with green text had minimal expression change and a log2 increase in methylation greater than 0.5. Genes labeled with red text had a significant increase in expression and a log2 decrease in methylation greater than 0.5. Selected genes labeled with black text were also analyzed by QRT-PCR or MassARRAYTM (Figs. S2 and S3). CASP14 is labeled due to its high expression in CAsE-PE relative to RWPE-1.

It was readily apparent that many ZNFs were among the genes that were both hypermethylated and downregulated. Twelve ZNF genes (ZNF542, ZNF300, ZNF560, ZNF454, ZNF471, ZNF667, ZFP28, ZFP42, ZSCAN18, ZNF677, ZIK1, and ZNF582), 9 of which are KRAB-ZNFs, demonstrated significant downregulation and increased DNA methylation (> 0.5 log2 fold change). Overall, silencing of ZNF genes represented a substantial portion (12 of 46, ~26%) of all the significantly downregulated genes in CAsE-PE that were correlated with increased DNA methylation.

H3K27me3 domains are targets of aberrant DNA methylation

After merging the expression and methylation data sets it was apparent that many of the DNA hypermethylation events were not associated with changes in gene expression, since the target genes themselves were not expressed in parental RWPE-1 cells. In fact, genes associated with some of the most extensive DNA hypermethylation events, such as PHOX2B and DBC1, experienced no expression change (Fig. 1, green text). VSX1, SLITRK3, HOXC12, and PAX3 are other examples of substantial DNA hypermethylation events with no expression change, which is understandable since this group of genes is expressed at low levels in the parental RWPE-1 cells (Supplemental Data Set 2). Chromatin immunoprecipitations determined that VSX1, HOXC12 and PHOX2B were marked with H3K27me3 in RWPE-1 and CAsE-PE (Fig. 2). Overall, 92% of the genes with increased promoter methylation were not downregulated. Of these, 61% correspond to hESC H3K27me3 domains and 11% correspond to hESC H3K9me3 domains. While only a select few genes were verified as H3K27me3 targets in our models, based on the large overlap between hESC H3K27me3 domains and CAsE-PE hypermethylated genes (Fig. S1B) it is likely that a high proportion of the hypermethylated genes are H3K27me3 marked in the parental RWPE-1. These data further support that H3K27me3 is a major factor in directing aberrant DNA methylation during malignant progression, but the phenotypic effects of this phenomenon are unknown since it has not been associated with gene expression changes.

Figure 2. H3K27me3 target genes. H3K27me3 ChIPs were quantified by QRT-PCR. Error bars represent the SEM of three experiments. GAPDH primers were used as a control to show ChIP levels of a H3K27me3 negative region.

Zinc finger gene clusters associate with hESC H3K9me3 domains and are silenced by DNA methylation

In contrast to aberrant DNA methylation targeting H3K27me3 pre-marked silent genes, we found that actively expressed ZNF genes that overlap hESC H3K9me3 domains are preferred targets of DNA methylation linked gene silencing. We found a major overlap of hESC H3K9me3 domains with ZNF genes, many of which were downregulated and DNA hypermethylated. Most of the human ZNF genes reside within one of the 81 ZNF gene clusters across the human genome that are known to associate with H3K9me3 and heterochromatin binding proteins., Chromosome 19 is particularly rich with ZNF genes, accounting for about 44% of all the human ZNF genes including the largest ZNF gene cluster (19.13) with 76 ZNF genes. We found that 9 of the ZNF genes in cluster 19.13 were significantly hypermethylated while 10 were significantly downregulated in CAsE-PE relative to RWPE-1 (Fig. 3). Of the downregulated ZNFs in this cluster, 5 are hypermethylated, while ZFP28 and ZNF667 have increased promoter DNA methylation that did not reach the strict DMR requirement (Fig. 3). ZNF cluster 19.13 and most other ZNF gene clusters are hESC H3K9me3 domains (Fig. 4A). ZNF gene cluster 19.7 was also significantly downregulated and associated with H3K9me3 stem cell domains, but the methylation data for this cluster is limited because it was poorly covered by the MeDIP-on-Chip microarray (Fig. 4A). Four other ZNF clusters on chromosome 19 (19.12, 19.11, 19.6, and 19.5) also had downregulated and hypermethylated ZNFs. Downregulated, hypermethylated ZNFs were also found on other chromosomes, for example ZNF454, which is part of ZNF gene cluster 5.1 on chromosome 5, was hypermethylated, downregulated and associated with H3K9me3 stem cell domains (Fig. 4B). Likewise, ZNF300 on chromosome 5 was hypermethylated and downregulated, but does not belong to a ZNF gene cluster (Fig. 4B). In summary, while the data are limited, they suggest that ZNFs with unique connections to hESC H3K9me3 domains, are widely targeted for aberrant DNA hypermethylation and downregulation during toxicant-induced malignant transformation.

Figure 3. Venn diagram of ZNF genes in ZNF gene cluster 19.13 that are hypermethylated (i.e., contain a hyper DMR) or are downregulated.

Figure 4. Expression and Methylation Changes in CAsE-PE. Significantly upregulated genes or downregulated genes are plotted in red or blue respectively while other genes are shown in gray. Purple stars show the methylation change at differentially expressed genes where available, measured within a 2kb window centered on the TSS. Grey boxes mark the ZNF gene family clusters and green boxes mark H3K9me3 stem cell domains. (A) Chromosome 19. (B) Portion of Chromosome 5q.

Aberrant spreading of H3K9me3 is linked to DNA methylation and gene silencing

One would expect little to no transcription from genes associated with H3K9me3 since this histone modification induces heterochromatin formation and gene silencing. Therefore, it would seem contradictory for ZNF genes to be expressed in RWPE-1 cells and associated with domains commonly marked with H3K9me3 in other cell types. As a resolution to this contradiction, work by Blahnik et al. showed that actively expressed ZNF genes are marked by H3K9me3 primarily at their 3′ exons rather than through the entire gene. Based on this information we analyzed the histone methylation patterns at the 5′ and 3′ ends of several ZNF genes that were silenced by DNA methylation in CAsE-PE. Histone methylation levels were measured by ChIP at ZNF542, ZNF471 and ZFP28, all of which pertain to ZNF gene cluster 19.13 (Fig. 5). In parental RWPE-1 cells we found that the 3′ end of each of these genes was marked with H3K9me3 whereas at the 5′ ends of the genes H3K4me3 was present and H3K9me3 was absent. In the malignantly transformed CAsE-PE cells, H3K4me3 had been replaced by H3K9me3 at the 5′ ends of all three genes (Fig. 5). This same pattern was found at the ends of ZNF560 in ZNF cluster 19.5 (Fig. 6). This pattern of histone methylation suggests that during arsenic induced malignant transformation, H3K9me3 spreads from the 3′ exons of expressed ZNF genes to the 5′ end near the promoter leading to gene silencing and DNA methylation. Similar changes to histone methylation patterns were also observed at ZNF genes on chromosome 5 which experienced a replacement of H3K4me3 with H3K9me3 at their 5′ ends, DNA hypermethylation and downregulation in CAsE-PE (Fig. 6). Together, these data show that in malignant cells, individual ZNF genes and entire ZNF gene clusters across different chromosomes are silenced by spreading of H3K9me3 and DNA hypermethylation.

Figure 5. ZNF cluster 19.13 ChIPs. ChIPs of 5′ and 3′ ends of three ZNF genes in cluster 19.13 were quantified by QRT-PCR. (Top) H3K4me3. (Middle) H3K9me3. Error bars represent SEM of three experiments (*p = 0.05, Wilcoxon rank sum test). (Bottom) Locations of the ChIP amplicons.

Figure 6. ChIPs of 5′ and 3′ ends of ZNF560 and 5′ ends of ZNF300 and ZNF454 were quantified by QRT-PCR. (Top) H3K4me3. (Bottom) H3K9me3. Error bars represent the SEM of three experiments (*p = 0.05, Wilcoxon rank sum test)

Discussion

In this study we measured 1035 differentially methylated and 942 differentially expressed genes in a prostate epithelial model of malignant transformation by an environmental toxicant. Most of the hypermethylated genes had no gene expression change and this is most likely a result of H3K27me3 targeted aberrant DNA methylation. In the malignant prostate epithelial cells there was widespread downregulation of ZNF genes that was mediated by H3K9me3 spreading from the 3′ ends of ZNF genes to the ZNF gene promoters. Loss of expression of ZNF genes was associated with a replacement of H3K4me3 with H3K9me3 and DNA methylation at the ZNF gene promoters. These findings suggest a mechanism in which the boundaries of H3K9me3 domains are disrupted, allowing H3K9me3 to spread and direct de novo DNA methylation to the promoters of ZNF genes in particular.

While H3K27me3 may be linked to more DNA hypermethylation events than H3K9me3 (Fig. S1B), the effect of aberrant DNA hypermethylation at H3K27me3 marked genes is unclear. One reasonable possibility is that aberrant DNA hypermethylation adds a more permanent level of repression to already silent genes, causing them to become non-inducible in cancer. In contrast to H3K27me3, H3K9me3 linked aberrant DNA methylation is clearly connected with gene expression changes and therefore could be an important contributor to the malignant phenotype by reducing global expression of ZNF genes.

Individual ZNFs and even some clusters of ZNFs have often been found hypermethylated and silenced in several tumor types, suggesting that our findings represent a commonly disrupted epigenetic pathway in cancer progression that is not a rare event.- Perhaps the most compelling evidence that silencing of ZNFs in the arsenite-transformed prostate epithelial model is a relevant event in human carcinomas are the publicly available TCGA gene expression and DNA methylation data sets from normal prostate and prostate adenocarcinoma tissues (http://cancergenome.nih.gov/). The TCGA data show that DNA methylation-linked silencing of ZNFs is common (Fig. S4). In addition, large deletions along chromosome 19 are common in several malignancies including cervical cancer, esophageal squamous cell carcinoma and nasopharyngeal carcinoma.- Together, these data suggest that silencing of ZNFs is an important event in the progression of various cancers and is not an event unique to arsenic-associated malignancies.

In the present study we identified 34 downregulated ZNFs and therefore we can only speculate on what the consequences of such a broad silencing event would be. Several studies have identified individual ZNFs as potential tumor suppressors (including ZNFs 382, 569, 383, 649, 540, 446, 411, 418, and 322) that regulate cellular proliferation by inhibiting MAPK signaling and repressing various oncogenes such as MYC, MITF, HMGA2, and CDK6.,,- Together, these studies suggest that widespread silencing of ZNFs could lead to increases in MAPK signaling and proliferation. Another interesting possibility arises from recent work which has shown that ZNF genes co-evolved with human endogenous retroviruses (HERVs) and that one of their functions is to recruit epigenetic machinery to HERV sequences in the genome, which leads to HERV silencing by histone and DNA methylation.,- This suggests that one possible effect of widespread ZNF gene silencing could be hypomethylation and de-repression of HERVs, events that are commonly found in tumors (reviewed in ref. 41).

Finally, results from our ChIP assays indicated that in malignant cells, H3K9me3 had spread into previously unoccupied genomic regions. One potential mechanism underlying aberrant H3K9me3 spreading in malignant transformation involves the CTCF/PARP-1 insulator elements that are found at the boundaries between inactive and active genomic regions. Previous studies have shown that CTCF and PARP-1 together, protect certain genomic regions from acquiring repressive epigenetic marks., Oxidative stress, a common effect of arsenic exposure, disrupts the CTCF/PARP-1 interaction and may lead to increases in DNA methylation and other repressive epigenetic marks at the corresponding genomic loci., Therefore, disruption of CTCF/PARP-1 insulator elements by arsenical induced ROS is one possible explanation of how H3K9me3 could spread from the 3′UTR of ZNF genes to the promoter regions, but future studies will have to show whether this hypothesis holds any merit.

Materials and Methods

Cell culture and sample acquisition

The samples used in this study along with their culture conditions were described previously. The non-tumorigenic RWPE-1 cell line was derived from prostate epithelial cells that were infected with human papillomavirus 18 to induce immortalization. CAsE-PE (arsenite transformed RWPE-1) cells were derived and cultured as described previously., STR profiling confirmed that RWPE-1 and CAsE-PE are from the same individual.

Nucleic acid isolation

Genomic DNA was isolated using the DNeasy Blood and Tissue Kit according to the manufacturer’s protocol (Qiagen). Total RNA was harvested using TRIzol and purified using the miRNeasy Kit (Qiagen). The quantity of each sample was assessed using absorbance at 260 nm on the NanoDrop 1000 Spectrophotometer.

MeDIP coupled to microarray (MeDIP-on-Chip)

Methyl-cytosine DNA Immunoprecipitations (MeDIP) were performed as previously described. MeDIP and input samples were amplified, labeled, co-hybridized to Human Promoter ChIP-on-Chip Microarrays Set 244K (Agilent Technologies, G4489A) and washed as previously described.

MeDIP-on-Chip data analysis

Microarray data were imported to R and normalized as previously described. M values (log2 ratios of immunoprecipitated and input channel) were used for further analysis as a measure of enrichment of a region centered on individual probes. For calculation of differentially methylated regions (DMRs), M values were analyzed in a sliding window of 1200 bp, with a step of one probe. DMRs were defined as regions of at least three consecutive probes less than 600 bp apart where the mean difference of ratios was at least 1.5 fold. A p-value cut was selected to maintain the false discovery rate (FDR) ≤ 5%. The FDR was determined by analysis of permutated data. DMRs were calculated between RWPE-1 and CAsE-PE samples. The coverage of the Human Promoter ChIP-on-Chip Microarray Set 244K was annotated to individual genes and a gene was considered differentially methylated when its annotated region contained a DMR. Gene set enrichment analyses of hypermethylated genes were performed using the GOstats package in R.

Analysis of gene expression change

Three biological replicates of RWPE-1 and 3 replicates of CAsE-PE were grown until they reached approximately 80% confluence in 25 cm2 flasks. Total RNA was isolated from each flask and concentrations and quality were determined using absorbance at A260 and A280 followed by analysis on an Agilent Nanochip (Agilent Technologies). RNA samples were used to produce labeled target, hybridized to Affymetrix GeneChip® Human Gene 1.0 ST Arrays, and read using an Affymetrix scanner according to the manufacturer’s protocols. Raw data (CEL files) were normalized and summarized according to Irizarry et al., as implemented in the package aroma.affymetrics. Differential expression was tested using the package LIMMA in the R programming environment. All p-values were adjusted according to Benjamini and Hochberg’s method to control the false discovery rate using a global method. Differential gene expression was defined as a fold change ≥ 2 and adjusted p-value ≤ 0.05. Several gene expression changes were confirmed by quantitative real-time (QRT)-PCR as described previously.

Merging DNA methylation data with gene expression data

MeDIP-on-Chip data were merged with gene expression microarray data by calculating the average methylation differences within 2 kb windows centered on TSSs of genes that had 5 or more MeDIP-on-Chip probes within the 2 kb window.

Chromatin immunoprecipitation

Chromatin immunoprecipitations (ChIPs) were performed using the MAGnify Chromatin Immunoprecipitation system (Life Technologies, 49-2024) according to the manufacturer’s protocol with the following adjustments. Antibodies against H3K4me3 (Millipore, 05-745R), H3K27me3 (Cell Signaling, C36B11) and H3K9me3 (Diagenode, mAb-146-050) were used for chromatin immunoprecipitation. Cells were grown in 10-cm dishes, washed with Hank’s Balanced Salt Solution (HBSS) and then crosslinked by adding 5 mL of HBSS with 1% formaldehyde for 8 min at room temperature. The crosslinking reaction was stopped by adding 500 μL of 1.25 M glycine and incubating 5 min. Dishes were washed 2× with 5 mL HBSS and then the cells were scraped from the dishes in two steps using 2.5 mL of HBSS with 2% FBS, 0.1% EDTA and protease inhibitors for each scraping. Chromatin was sheared by sonication in a chilled Bioruptor (Diagenode) set to high for 42 cycles of 30 sec on/off. Approximately 200,000 cells were used for each ChIP. Input DNA was purified from a sample aliquot equal to 10% of the total cells used for each ChIP. Equal volumes of ChIP and input DNA were quantified by QRT-PCR using the ABI 7500 Real-Time Detection System. Primer sequences are provided (Table S2).