Global Analysis of DNA Methylation by Methyl-Capture Sequencing Reveals Epigenetic Control of Cisplatin Resistance in Ovarian Cancer Cell

Cisplatin resistance is one of the major reasons leading to the high death rate of ovarian cancer. Methyl-Capture sequencing (MethylCap-seq), which combines precipitation of methylated DNA by recombinant methyl-CpG binding domain of MBD2 protein with NGS, global and unbiased analysis of global DNA methylation patterns. We applied MethylCap-seq to analyze genome-wide DNA methylation profile of cisplatin sensitive ovarian cancer cell line A2780 and its isogenic derivative resistant line A2780CP. We obtained 21,763,035 raw reads for the drug resistant cell line A2780CP and 18,821,061reads for the sensitive cell line A2780. We identified 1224 hyper-methylated and 1216 hypomethylated DMRs (differentially methylated region) in A2780CP compared to A2780. Our MethylCap-seq data on this ovarian cancer cisplatin resistant model provided a good resource for the research community. We also found that A2780CP, compared to A2780, has lower observed to expected methylated CpG ratios, suggesting a lower global CpG methylation in A2780CP cells. Methylation specific PCR and bisulfite sequencing confirmed hypermethylation of PTK6, PRKCE and BCL2L1 in A2780 compared with A2780CP. Furthermore, treatment with the demethylation reagent 5-aza-dC in A2780 cells demethylated the promoters and restored the expression of PTK6, PRKCE and BCL2L1.


Introduction
Drug resistance is the major reason leading to the high death rate of ovarian cancer. The chemotherapeutic agent cisplatin (cisdiamminedi-chloroplatinum(II)) is particularly effective against ovarian carcinoma with an initial response rate of up to 70% [1]. However, ovarian cancer eventually develops resistant to cisplatin, and the 5-year survival rate for patients is only 15-20% [2]. Cisplatin mediates its actions by forming DNA adducts-primarily intra-strand crosslink adducts [3] and activates several signal transduction pathways include the ATR, p53, p73, and MAPK pathways, resulting in apoptosis [3,4].
DNA methylation is often associated with transcriptional repression of gene expression [5] and with responses to chemotherapy [6,7]. One classical example is that methylation of the MGMT (O6-methylguanine-DNA methyltransferase) promoter in gliomas is a useful predictor of the responsiveness of the tumors to alkylating agent carmustine [1,3-bis(2-chloroethyl)-1nitrosourea], as well as of overall and disease-free survival in gliomas [6]. In ovarian cancers, Taniguchi et al. propose a model for ovarian tumor progression in which the initial methylation of FANCF is followed by FANCF demethylation and ultimately results in cisplatin resistance [7]. DNA methylation of several genes in ovarian cancers including HSulf-1 [8], ABCG2 [9], EZH2 [10] have been found to be associated with drug resistance. Boettcher et al. analyzed high-definition DNA methylation profiles of 800 CpG islands (CGIs) of selected genes and identified that hyper-methylation in CGIs of BRCA1, CDH1, DNAJC15 and SULF2 as well as hypo-methylation of CGIs for ABCB1, APC and HIC1 genes showed increased doxorubicin tolerance [11]. Chang et al. used global gene expression profiling to analyze cancer cells before and after treatment of DNA methyltransferase inhibitor5aza-29-deoxycytidine, which re-activates methylation silenced gene [12]. They identified several hundred genes that were downregulated in cisplatin resistant cancer cells and reactivated by the DNA methyltransferase inhibitor 5-aza-29-deoxycytidine [12]. Li et al. compared the methylation pattern of A2780 and their derived resistant cell line after several cycles of drug selections using global CGI methylation arrays and mRNA expression microarrays [13].
Taking advantage of the the next-generation sequencing (NGS) technologies, the methylated fraction of genome captured by MeDIP-seq [14], MethylCap-seq [15] and methylcap-seq [16] have been profiled in a greater depth than the array-based platform, revealing many novel regions that are differentially methylated with the biological contents. Here we report the comparison of the global DNA methylation patterns of ciplatin sensitive (A2780) and resistant (A2780CP) ovarian cell lines for the key epigenetic controls responsible for cisplatin resistance. We identified 1224 hyper-methylated and 1216 hypomethylated DMRs (differentially methylated region) in A2780CP compared to A2780. We also found that A2780CP, compared to A2780, has a lower global CpG methylation. Several genes were confirmed to have both promoter methylation and corresponding expression changes including PTK6, PRKCE and BCL2L1, and treatment with the demethylation reagent 5-aza-dC in A2780 cells demethylated the promoters and restored their expression.

Results
Methylomic analysis of the cisplatin resistant A2780CP and its isogenic cisplatin sensitive A2780 cell line Human ovarian carcinoma A2780CP cell line (cisplatinresistant) was derived from A2780 (cisplatin-sensitive) cell line, but show increased resistance to cisplatin [17]. A2780CP and A2780 are isogenic (same genomic DNA). To confirm that the cell line A2780CP we had still maintained cisplatin resistance, we performed MTT assay to evaluate cisplatin drug responses of A2780CP and A2780 cells. We found that the IC50 of A2780CP (5.0 mg/ml) is nearly 3-fold higher than that of A2780 (1.8 mg/ml) (Fig. 1). Our data is consistent with the previous reported cisplatin response profiles of A2780CP [17], indicating that the cell lines we had in the laboratory retains the difference in cisplatin resistance.
We have therefore applied MethylCap-seq to analyze the methylomic profiles of both A2780 and A2780CP cells. We obtained 21,763,035 raw reads for the drug resistant cell line A2780CP and 18,821,061reads for the sensitive cell line A2780. 13,950,202 (64.1%) and 13,671,899 (72.6%) reads were uniquely mapped to human genome for A2780CP and A2780 respectively.
We annotated the reads with respect to the CpG islands in the human genome (http://genome.ucsc.edu/) and found that about 1.4% and 2.5% of the reads locate inside CpG islands respectively for A2780CP and A2780.The average sequence depth for the CpG islands covered is about 16 and 24.5 times. About 68% and 66% of the human CpG islands were covered in our analysis for A2780CP and A2780 respectively ( Table 1). The CpG coverage and depth we obtained in our MethylCap-seq is similar to what previously published [18].
To detect differential methylation regions (DMRs) in the human genome between the sensitive and resistant cells, we used MEDIPS, a recently developed software tool specialized for analyzing immunoprecipitation based methylation analysis (e.g. MeDIP-seq and MethylCap-seq) [19]. We set the criteria for significant DMRs: the length of peaks was set to 500 base pairs and peaks with .20 rpm (reads per million), p-value less than 0.001 and ratio of rpm between two cell line .20. We obtained 1224 DMRs that is hyper-methylated in A2780CP compared to A2780 (Table S1) and 1216 DMRs that is hyper-methylated in A2780 compared to A2780CP (Table S2). These DMRs were annotated with their genomic locations and associated genes within 25 K bp to +5 K bp to the transcription start sites (Table  S3 and 4).
Nearly half of DMRs were located to within 5 k bp of transcript start sites (TSS) of genes. 190 and 270 DMRs were located to upstream of TSS respectively in A2780 and A2780CP ( Table 1). The rest mapped to other regions of the human genome including introns, exons, intergenic regions and repeats (Ensembl human genome annotations) ( Fig. 2A). Hypermethylated regions in promoters of genes are known to decrease gene expression [17], but the effect of hypermethylation regions elsewhere within the gene region remains inconclusive. The hypermethylation of repetitive genomic sequences might prevent chromosomal instability, translocations, and gene disruption caused by reactivation of transposable DNA sequences [17,20].
Li et al. compared isogenic A2780 and their selected resistant cell line using the 60-mer oligo microarray containing 40,000 CpG-rich fragments from 12,000 known gene promoters (Agilent, Santa Clara, CA) [13]. Using a cutoff value of 1.5 fold, they identified 595, 870 and 1,176 hypermethylated genes for the Round1, Round3 and Round5 resistant sublines comparing to the parental (''Round0'') A2780 cells [13]. The technology they used is differential methylation hybridization (DMH) [21], in which methylated DNA was separated from the unmethylated by methylation sensitive enzyme diagestion for probing the oligonucleotide array of the human CpG island regions. DMH is different from MDB-seq, in which methylated DNA was precipitated by recombinant methyl-CpG binding domain of MBD2 protein, and then sequenced. In DMH, isolated DNA was digested with the methylation-insensitive restriction enzyme BfaI (C ' TAG), ligated to linkers, and digested again with the methylation sensitive enzymes HinP1I (G ' CGC) and HpaII (C ' CGG). The digestion products were then amplified using linkers and labeled with Cy3 or Cy5 dyes for comparative hybridizations [21]. Despite different technologies used, we compared our MDB-seq data with Li et al.'s DMH data. We found that 221 (of 1224, 18.1%) and 142 (of 1216, 11.68%) of the hypermethylated and hypomethylated genes in A2780CP in comparison with A2780 were also on the array that Li et al. used (Table S3,S4). We found that 15 regions (corresponding to 13 genes) and 23 regions (21 genes) that are common between the two data sets (Table 2), which are highly significant overlaps with hypergeometric probabilities of 1.11E-5 and 4.22E-20 respectively. The commonly hypermethylated genes in the resistant cells include retinoblastoma binding protein 8 (RBBP8), SRY-box 1 (SOX1), wingless-type MMTV integration site family (WNT9A), general transcription factor IIIA (GTF3), and the commonly hypomethylated genes in the resistant cells include solute carrier family 22 (extraneuronal monoamine transporter) (SLC22A3), aldehyde dehydrogenase 1 family (ALDH1A3), hyaluronan synthase 3 (HAS3) and CUB domain containing protein 1 (CDCP1) ( Table 2).
There are 410 hypermethylated DMRs (Table S3) and 316 hypomethylated (Table S4) in A2780CP in comparison with A2780 that were not on the array that Li et al. used, and would not have been identified by the DMH approach. These new DMRs revealed the power of MDB-seq in identifying novel DMRs without prior knowledge of candidate promoters to be put on arrays to be used in DMH. Furthermore, these new DMRs might be due to differences in the derivative resistant cell lines that were used by Li et al. and us. A lower global CpG methylation in the cisplatin resistant A2780CP cells compared to the cisplatin sensitive A2780 cells To understand the global pattern of CpG methylation in both cell lines, we analyzed the characteristics of hypermethylated CpGs across the entire genome by counting the number of neighboring CpGs within the 500 bp regions surrounding hypermethylated sites and calculated the CpG observed/expected ratio (CpG o/e ) for each CpG using the method described by Ruike et al. [22]. We found that, in general, A2780CP has lower CpG o/e than that of A2780, suggesting a lower global CpG methylation in A2780CP compared to A2780 (Fig. 2B, top left panel). The differences were mostly manifested in the higher CpG o/e in the intron and repeat sequences of A2780 compared to that of A2780CP. Comparing CpG o/e across different genomic region categories (Fig. 2B), the CpG o/e was greater than 0.6 in exons and promoters for both resistant and sensitive cell lines, but less than 0.6 in other genomic regions (repeats, introns and intergenic regions) (Fig. 2B).

Functional annotations and pathway enrichment analysis for the DMR genes
We performed GO analysis for hypermethylated genes in both the resistant (Table S5) and sensitive cell line (Table S6). We found that hypermethylated genes in the sensitive cell line were enriched for GO terms: anti-apoptosis (GO:0006916), negative regulation of cell death (GO:0060548), and regulation of intracellular protein kinase cascade (GO:0010627). The hypermethylated genes in the resistant cell line have several enriched terms related to transcript factor regulation such as GO:0030528 transcription regulator activity and GO:0003677 DNA binding.
We also performed KEGG pathway analysis for the DRM genes and found that genes with hypermethylated promoters were enriched in the following KEGG pathways including 11 genes in hsa04010 (MAPK signaling pathway), 18 genes in hsa05200 (Pathways in cancer), 5 genes in hsa04310 (Wnt signaling pathway), 5 genes in hsa00982 (Drug metabolism), and 5 genes in hsa04115 (p53 signaling pathway) (Table 3).
Interestingly, many ECM-related and membrane channel proteins were hypermethylated. For example, KCNS1 and KCNA1, two potassium voltage-gated channel protein, and SLC15A4 (solute carrier family 15, member 4) are hypermethylated in A2780CP cells while SLC4A11 (solute carrier family 4, sodium borate transporter, member 11) is hypermethylated in A2780 cells. COL18A1 (collagen, type XVIII, alpha 1) is hypermethylated in A2780 cells while KRTAP10-6 (keratin associated protein 10-6) and LRFN5 (leucine rich repeat and fibronectin type III domain containing 5) are hypermethylated in A2780CP cells. We also identified hypermethylation at several microRNA loci including hypermethylation of MI6A2, MIR129-2, MIR124-1, MIR124-3, and MIR10A in A2780CP and hypermethylation of MIR185, MIR548Q, MIR642A and MIR661 in A2780 (Table 4). Validation of genes with DMRs by MS-PCR and bisulfite sequencing The expression of genes with hypermethylated promoter was checked by RT-qPCR, if the expression was significant changed then methylation specific PCR (MS-PCR) was used to validate the methylation status of DMR regions between two cell lines. We identified the following genes BCL2L1 (BCL2-like 1), PPKCE (protein kinase C, epsilon), PTK6 (PTK6 protein tyrosine kinase 6), RAC2 (ras-related C3 botulinum toxin substrate 2), SECTM1 (secreted and transmembrane 1) and DDIT3 (DNA-damageinducible transcript 3) that changed in both promoter methylation and expression. Figure 3A showed the expression changes of these genes, and Figure 3B showed the promoter methylation patterns of these genes. DDIT3 was hypomethylated in A2780 cells compared to A2780CP, while the rest showed the opposite (Fig. 3B). To further confirm DMRs between A2780 and A2780CP, we used bisulfite sequencing to sequence the promoter regions of three randomly picked gene from the 6 genes. Figure 3C showed that all promoter regions of the picked gene PTK6, PRKCE and BCL2L1 were hypomethylated in A2780CP compared with A2780.

Restoration of gene expression of DMR affected genes by treatmentwith 5-aza-dc demethylation reagent
In order to further confirm our data, we used 5-aza-29deoxycytidine (5-aza-dC) to see whether we can re-activate methylation-silenced genes that we identified. 5-aza-29-deoxycytidine (5-aza-dC) is an analogue of cytosine. When incorporated into DNA, it irreversibly binds the methyltransferase enzymes, resulting in passive de-methylationand reactivation of epigenetically silenced genes [23]. We used 5-aza-dC to treat both A2780CP and A2780 cells in different doses from 0 mM to 10 mM. The expression of genes treated with the lowest (0 mM) was compared to that treated with the highest (10 mM) dose of demethylation reagent. We showed that we were able to demethylate the promoters of the hypermethylated PRKCE, BCL2L1, RAC2, PTK6, SECTM1 (Fig. 4A,B), and restored their expression after treatment with 5aza-dC in A2780 cells. Similarly, we were able to demethylate the promoters of the hypermethylated gene DDIT3, and restored its expression after treatment with 5-aza-dC in A2780CP cells (Fig. 4A,B).

Discussion
We describe here for the first time an MDB-seq analysis of a cisplatin response model of ovarian cancer cells. The epigenetic changes between the isogenic pair of the cisplatin sensitive A2780 and its resistant derivative A2780CP cells suggest an epigenetic control mechanism for cisplatin resistance. The global MDB-seq data set generated for this pair of isogenic cells will be a useful resource for the research community interested in cisplatin resistance and epigenetics, and in ovarian cancer in general. Global DNA methylation changes during carcinogenesis and is often a predictor of therapy responses. For example, Anisowicz et al. showed that even the normal part of the lung from a cancer patient has already experienced a loss of global DNA methylation compared to a normal individual [24]. Kim et al. showed that global CpG methylation plays a role in epigenetic control of the radiosensitivity in lung cancer cell lines [25]. In gliomas, LINE-1 methylation, a global DNA methylation marker, is proportional to MGMT promoter methylation and higher LINE-1 methylation is a favorable prognostic factor in primary GBMs, even compared to MGMT promoter methylation [26]. Here we observed a lower global CpG methylation in the cisplatin resistant A2780CP cells compared to the cisplatin sensitive A2780 cells, suggesting that global DNA methylation patterns might be able to predict cisplatin resistance for ovarian cancers. Further experimentation is necessary to evaluate this hypothesis. We obtained 1,224 and 1,216 DMRs that is hyper-methylated or hypo-methylated in A2780CP compared to A2780 cells (Table S1, S2). We found that hypermethylated genes in the sensitive cell line were enriched for GO terms for anti-apoptosis (GO:0006916) and negative regulation of cell death (GO:00 60548), suggesting a possible general mechanism of epigenetic silencing of anti-apoptosis genes, resulting in sensitivity to cisplatin. We also found that the DMRs are enriched in several signaling pathways including the hsa04010 (MAPK signaling pathway), hsa04310 (Wnt signaling pathway), and the hsa04115 (p53 signaling pathway) pathways ( Table 3). The Wnt signaling pathway has been implicated in ovarian cancer progression and chemoresistance [27]. Activation of JNK/p38 MAPK pathways in response to cisplatin could lead to Fas ligand induction and cell death in ovarian carcinoma cells [4]. P53 and its signaling pathway were also been implicated in cisplatin resistance in ovarian cancer cells [28,29]. Our data suggest that epigenetic regulation of these signaling pathways is one of the mechanisms for the involvement of these pathways in cisplatin resistance in ovarian cancer cells. We also identified hypermethylation at several microRNA loci including hypermethylation of MI6A2, MIR129-2, MIR124-1, MIR124-3, and MIR10A in A2780CP and hypermethylation of MIR185, MIR548Q, MIR642A and MIR661 in A2780 (Table 4). MicroRNAs have been implicated in cisplatin resistance in ovarian cancer [30]. For example, Yang et al. showed that many miRNAs are deregulated in human ovarian cancer including miR-214, miR-199a*, miR-200a, miR-100, miR-125b, and let-7 cluster, and that miR-214 induces cell survival and cisplatin resistance by targeting PTEN [30]. Epigenetically silenced microRNAs have been implicated in cancers [31,32]. However, methylation of miRNA locus has not been reported previously for ovarian cancers; neither are reports on the relationship between miRNA methylation and cisplatin resistance. Our data might point to a new direction for the study of methylation of microRNA loci and cisplatin resistance.
We found that BCL2L1 (BCL2-like 1), PPKCE (protein kinase C, epsilon), PTK6 (PTK6 protein tyrosine kinase 6), RAC2 (ras-related C3 botulinum toxin substrate 2), SECTM1 (secreted and transmembrane 1) are hypermethylated in the cisplatin sensitive A2780 cells and DDIT3 (DNA-damage-inducible transcript 3) is hypermethylated in the cisplatin resistant A2780CP cells. Their expression also were also changed accordingly to the hypothesis that hypermethylation would silence gene expression. We confirmed the methylation pattern of these genes by MS-PCR, and also were able to demethylate the promoters of these genes and restore their expression after treatment with 5-aza-dC, an agent that demethylates and reactivates of epigenetically silenced genes [23].
We found BCL2L1 is hypermethylated in the sensitive A2780 cells. BCL2L1 (bcl-xl) encodes a protein belonging to the BCL-2 protein family, whose protein members act as anti-or proapoptotic regulators [33]. Over-expression of the Bcl-xL protein is known to confer resistance to a broad range of potentially apoptotic stimuli in carcinogenesis processes including oncogene activation, hypoxia and matrix detachment [34][35][36][37]. Our result reveals that the hypomethylation of BCL2L1 in the resistant cells (i.e. hypermethylation of BCL2L1 in the sensitive cells) would result in increased BCL2L1 expression, thus conferring resistance to cisplatin. Our data is consistent with the observation by Williams et al. that expression of Bcl-xL in ovarian carcinoma is associated with chemoresistance and recurrent disease [38]. Taking together, this suggest that modulating the expression of BCL2L1 (Bcl-xL) by epigenetic means might be a way to overcome cisplatin resistance in ovarian cancer cells.
We also identified PTK6 (protein tyrosine kinase 6) as hypermethylated in the cisplatin sensitive A2780 cells compared to A2780CP cells. PTK6 directly phosphorylates AKT and promotes AKT activation in response to epidermal growth factor [39]. The relationship of PTK6 with cisplatin has not been previously studied. Further functional analysis of PTK6's role in modulating cisplatin resistance is warranted. PPKCE (protein kinase C, epsilon) is another gene that is hypermethylated in the cisplatin sensitive A2780 cells compared to A2780CP cells. PPKCE is a member of the protein kinase C (PKC) family, whose members phosphorylate a wide variety of protein targets and are involved in diverse cellular signaling pathways [40]. Interestingly, cisplatin was able to induce phosphorylation and translocation of PPKCE from the plasma membrane to the nuclear membrane and to the cytosolic fraction [41]. The hypermethylation of PPKCE in the sensitive cells would reduce its expression, resulting in less translocation to the nuclear membrane or cytosolic fraction upon addition of cisplatin. However, whether reduced expression of PPKCE confers sensitivity to cisplatin in ovarian cancer cells remains to be studied. Finally, we identified DDIT3 (DNA-damage-inducible transcript 3) as hypermethylated in the cisplatin resistant A2780CP cells compared to A2780 cells. Also named GADD153 (Growth arrest and DNA damage-inducible protein 153), DDIT3 encodes a member of the CCAAT/enhancer-binding protein (C/EBP) family of transcription factors, and is activated by endoplasmic reticulum stress, and promotes apoptosis. DDIT3 (GADD153) expression could be induced by cisplatin [42]. We observed that it is hypermethylation in the resistant cells, which would result in reduced expression of DDIT3. It is possible that cisplatin acts through DDIT3 to promote growth arrest and apoptosis, and reduced expression of DDIT3 in the cells would make them more resistant to cisplatin. However, the detailed mechanism remained to be investigated.
In summary, we generated a global dataset for an ovarian cancer cisplatin resistance model and identified several genes that were subjected to epigenetic control to modulate cisplatin resistance. These genes might serves as targets to overcome chemoresistance to cisplatin in ovarian cancers.

Cytotoxicity analysis
Both resistant and sensitive cells were seeded into 96 well plates at a concentration of 4000 cells/well in five replicates with complete culture medium. Then cells were treated with cisplatin in different dose from 0 mg/ml to 16 mg/ml. After 72 hours incubation, both viable and dead cells were counted by using MTT assay, and only the viable cells were included in data analysis.

Bisulfite-modified DNA sequencing and Methylation-Specific PCR
We prepared genomic DNA from cultured cells using the DNeasy Blood & tissue Kit (Qiagen). Approximately 200 ng of DNA was bisulfite-treated with the EZ DNA Methylation-Gold kit (Zymo Research) according to the manufacture's protocol. ZymoTaq Premix was used for Methylation-Specific PCR (MSP) and according to these primers (Table S7). Methylation-specific PCR was run in a total volume of 25 mL by using ZymoTaq Premix (ZYMO research). MSP reactions were subjected to initial incubation at 95uC for 10 minutes, followed by 40 cycles of 95uC for 30 seconds, and annealing at appropriate temperature for 35 seconds and 72uC for 40 seconds. Final extension was done by incubation at 72uC for 7 minutes. MSP products were separated on 2% agarose gels and visualized by ethidium bromide staining.
The primers for bisulfite-modified sequencing (Table S7) and MSP were designed by web tools MethPrimer (http://www. urogene.org/methprimer/index1.html) [43]. The conditions for bisulfite-modified sequencing reaction are same as MSP. The products were purified using MinElute PCR purification Kit, cloned into the pMD 18-T vector (Takara) and sequenced.

RNA isolation and Real time RT-PCR
RNA was extracted using the protocol for Trizol reagent (Invitrogen). 2 mg RNA was firstly reverse-transcribed in 25 mL with an Archive Kit (Applied Biosystems) and 2 mL were amplified by Real Time PCR (Bio-RAD CFX96 Rea-Time System). cDNA samples were amplified with SYBRH Pre-mix Ex Taq TM . The thermal cycling profile consisted of initial denaturation at 95uC for 30 seconds and 40 cycles at 95uC for 5 seconds, 60uC for 15 second, and 72uC for 30 seconds. Each sample was processed in triplicate.

Methylated DNA Binding Protein sequencing (MethylCap-seq)
3 mg of genomic DNA isolated as described above was fragmented by sonication with Bioruptor (Diagenode, Belgium) and end-repaired, A-tailed, and ligated to 2.5 mMol of ''pairedend'' adapters (IDT Inc.) following the manufacturer's recommend protocol (Illumina Inc.). 1.2 mg of the DNA was fractionated on a home-made GST-MBD resin with a buffer of a stepwise increased salt concentration [16]. The high salt fraction was directly amplified by 12-cycle PCR [18]. The 350-450 bp fraction of the PCR products was gel purified and quantified using an Agilent DNA 1000.
The 250-400 bp DNA fractions were excised and purified as described above. The products were assessed and quantified using an Agilent DNA 1000 series II assay and Qubitfluorometer (Invitrogen) respectively. Each library was diluted to 8 nM for sequencing on an Illumina Genome Analyzer II following the manufacture's recommended protocol. Obtained images were analyzed and base called using Illumina provided GA pipeline software OLB 1.6.0 with default setting. The raw reads from MethylCap-seq were submitted to GEO database which accession number isGSE31418.
The 250-400 bp DNA fraction of the amplified fraction was gel-purified. Each library was diluted to 8 nM for sequencing by Illumina Genome Analyzer II following the manufacture's recommended protocol. Obtained images were analyzed and base-called using Illumina provided GA pipeline software OLB 1.6.0 with default setting. The raw reads from MethylCap-seq were submitted to GEO database (accession number GSE31418).

Mapping of Sequence Reads
Human genome sequence and mapping information (Feb. 2009, GRCh37/hg19) was downloaded from ENSEMBL FTP site (http://asia.ensembl.org/info/data/ftp/index.html). Filtered reads were firstly cut into 36 bp to get higher quality reads then mapped to HG19 reference allowing up 2 mismatches by software SOAP2 [44]. The matched results from SOAP2 were converted into bed format by our own script.

Identification and Annotation of DMRs
To identify differential methylation regions (DMRs), the bed format files were used as input for the MEDIPS program [19]. The criteria for significant DMRs were: the length of peaks was set to 500 base pairs and peaks with .20 rpm (reads per million), pvalue less than 0.001 and ratio of rpm between two cell line .20. The length of each DMR was set to 500 bp. Genomic feature data was retrieved from ENSEMBL database (build 55), promoter was defined as the 5 kb upstream of transcription start sites. Each DMR was annotated according to its genomic location by our own script. The DMRs that locate within the proximal promoter (25 kb to +5 kb from their transcription start site) were annotated with GO, KEGG pathway and Reactome Pathway by package ChIPpeakAnno [45] in R and web tool IDConverter (http:// idconverter.bioinfo.cnio.es/) [46]. CpG Obs/Exp ratio was calculated by our own script for each DMR and plotted by the R language. GO Miner [47] was used to analyze GO enrichment and a P value ,0.01 was considered significant. Pathway enrichment was analyzed through KEGG database (http://www.genome.jp/kegg/tool/map_pathway1.html) [48]. To compare to Li et al.'s data, we downloaded the data ftp: //ftp.ncbi.nih.gov/pub/geo/DATA/SeriesMatrix/GSE15373/. When there are multiple probes corresponding to the same genes, the probe that showed the most differentially expressed were used for comparison. For calculating hypergeometric probabilities, the population size is set at 12,000 as there are 12, 000 genes in the array that Li et al. used.