Identification of DNA methylation changes associated with human gastric cancer

Gastric tissue samples

We obtained three snap-frozen gastric tumors and matched normal gastric tissue from Seoul National University College of Medicine for methylome study. Additionally, twenty-eight matched pairs of normal and tumor stomach tissues were obtained for further confirmation. All samples were obtained by endoscopic resection during examination of the patients who gave informed consent.

Methylated DNA recovery assay (MIRA)

Genomic DNA from 25 mg of gastric tissue was purified by using DNeasy Blood & Tissue Kit (Qiagen, Valencia, CA). Genomic DNA samples from 3 individuals were pooled at the same concentration. MIRA was carried out as previously described [30–32]. Briefly, GST-tagged MBD2b and His-tagged MBD3L1 proteins were prepared as described. 15 ug of genomic DNA was fragmented to 100 ~ 500 bp by sonication and incubated with 28 ug of purified GST-MBD2b protein, 28 ug of His-MBD3L1 protein and 7 ug of JM110 Bacterial DNA for 6 hours. 30 ul of MagneGST beads (Promega, Madison, WI) preblocked with 7 ug of JM110 bacterial RNA were added and incubated at 4°C with rotating for 45 minutes in final 600 ul of MIRA binding reaction mixture. Beads were washed three times with 1 ml of washing buffer, and methylated fragments were eluted by incubation at RT for 5 minutes and then 56°C for 30 minutes with 30 ul of TE containing RNase A (100 ug, Qiagen) and Proteinase K (15 ug, Qiagen). Eluted DNA fragments were further purified by using Qiaquick PCR purification kits (Qiagen).

Illumina Genome Analyzer sequencing

We used 10 ng of eluted DNA for Illumina Genome Analyzer sequencing. Following ligation of a pair of Solexa adaptors, ligation products with the maximum insert size of 200 bp were gel purified on 2% agarose and subjected to PCR amplification. Cluster generation and 36 cycles of sequencing were performed following the manufacturer's instructions. We sequenced 120 ul of adaptor-ligated, size-fractionated DNA (2 ~ 4 pM) on the Illumina Genome Analyzer. Sequence tags were mapped to the human genome (UCSC hg18 database based on NCBI Build 36.1 assembly) using the Solexa Analysis Pipeline (version 0.3.0). Sequenced reads of 34 bp (excluding the first and last nucleotide) that passed quality control filters were used.

Data processing and MES calculation

We extended the 3' end of the 34-bp reads by 200 bp to cover DNA fragments bound by the MBD proteins. The readout was converted to browser extensible data (BED) files for visualization in the UCSC genome browser http://genome.ucsc.edu/. We counted overlapping sequence tags at 50 bp resolution. To find enriched genomic regions, the number of mapped reads in a sliding window of 1 kb was compared to the total number of reads or the background number of reads in the genome. As such, MES was calculated in two ways; one is as the log2 of (target read count/target size)/(total read count/genome size) and floored to zero, the other is as the log2 of (target read count/target size)/(background read count/background size) and floored to zero. To adjust for background sequencing bias, MESbg was calculated in the same manner for input sequencing without affinity purification and subtracted from MES.

Genomic positions of CGIs, promoters, transcript bodies, CDSs, and repetitive elements

All genomic positions of CGIs, transcripts and repeat elements were downloaded from the UCSC genome browser. A total of 27,639 CGIs (except randomly located CGIs) were predicted by the following criteria: GC content of 50% or greater, length greater than 200 bp, and ratio greater than 0.6 of observed number of CpG dinucleotides to the expected number [33]. The NCBI mRNA reference sequences collection (RefSeq from release version 46; March 11, 2011) was used for identifying transcription units with the defined transcription start, end sites and CDS start, end sites. For promoters, we used the region 500 bp upstream ~ 500 bp downstream of the transcription start site. We obtained ~ 5 million repeat locations that had been determined by the RepeatMasker program based on the RepBase library of repeats.

Methylation level of genomic elements

The methylation level of a CGI, promoter, gene-body, and repeat element was estimated by means of MES overlapping each element. MES = 0 was used to define unmethylated elements. To measure hypermethylation or hypomethylation in cancer, we calculated the differential MESs as (Cancer MES - Normal MES). Differential MES > 1.0 was used as a threshold. To understand the functions of selected genes, we used the ontology classification of genes through the DAVID Functional Annotation Clustering tool http://david.abcc.ncifcrf.gov/.

Gene expression analysis

The microarray product used in this study was Codelink Human Whole Genome 55 K chip (GE Healthcare, USA). All experimental procedures including cRNA target preparation, hybridization, post-hybridization dye coupling were performed using vendor recommended protocols. The result files were imported into GeneSpring GX 7.3 (Agilent Technologies, USA) for filtering and basic statistical analysis. Among 55 K genes on the microarray, only the genes with present flags in at least 50% of samples were selected for subsequent analysis. The microarray data were deposited at the GEO http://www.ncbi.nlm.nih.gov/geo/ (accession number GSE33651).

MIRA and real-time qPCR

MIRA was performed on four additional individual samples. DNA was purified from the supernatant and monitored by real-time qPCR using Roche 480 machine. The sequences of used primers are presented in Additional file 1: Table S1.

Bisulfite treatment, methylation-specific PCR and pyrosequencing

We isolated the genomic DNA from individual sample by using a Qiagen DNeasy Tissue Kit (Qiagen). Bisulfite treatment was carried out using the EZ DNA methylation gold kit (Zymo research) according to the manufacturer's instructions. Bisulfite-treated DNA was stored at -80°C until further use. The primers used for MSP were designed using Methprimer [34], and are shown in Additional file 1: Table S1. PCR was performed with HotStarTaq Polymerase (Qiagen) and included an initial incubation at 95°C for 15 min, followed by 40 cycles of 95°C for 1 min, 59°C for 1 min and 72°C for 40 sec, followed by one cycle of 72°C for 10 minutes. MSP products were separated on 2% agarose gels and visualized by ETBR staining. The pyrosequencing reactions were automatically performed with a PSQ 96 system (Pyrosequencing AB) according to the manufacturer's instructions. Briefly, the biotinylated PCR product (50 ul) was purified by using streptavidin-sepharose beads (Amersham Biosciences). The purified product was loaded into the reagent cartridge with the enzyme, substrate and dNTP included in the PSQ96 SNP Reagent Kit (Pyrosequencing AB). The sequencing primers for pyrosequencing are shown in Additional file 1: Table S1.

Processing of MIRA-seq methylome data

We purified the methylated DNA enriched through MIRA (methylated CpG island recovery assay) and sequenced the DNA using next-generation sequencing. DNA methylation levels were determined using sequencing read counts of the corresponding regions, at 50 bp intervals, as described under Methods. We created DNA methylation maps for both normal and cancerous gastric tissues. For each sample, we obtained about 10 million sequence reads (Additional file 1: Table S2). Each methylome contained ~140 million CpG reads, covering ~48% of all genomic CpG sites excluding centromeres (Additional file 1: Table S3). The average coverage of CpG reads in each methylome was 4.5X. In support of the high sensitivity of MIRA, genomic segments containing only one CpG had higher read counts than those with no CpG (p value = 0), suggesting that single CpG changes could be resolved using MIRA. The average sequence reads increased in proportion to the number of CpGs within a 50-bp interval, and in fact, MIRA coverage was not low, even for regions of low CpG density (Additional file 2: Figure S1). Taken together, these results show that MIRA was successful in recovering a sufficient fraction of methylated regions. As for the accuracy of MIRA, ~99% of MIRA-captured fragments had at least one CpG site within their sequence, indicating a low false detection rate.

To measure enrichment of local methylation signals, we calculated methylation enrichment scores (MESs) by obtaining a read count in a given region and then performing normalization to control for the total read count (MESt) in the sample (global normalization) or the local read count (MESl) in a user-defined surrounding region (local normalization) (see Methods). This enables a direct comparison of independent samples with different read density. We then carried out a logarithmic transformation of the derived score. Along with having other mathematical merits, this provides the benefit of variance stabilization, particularly for high read counts, which are often coupled with high technical variations that may introduce significant bias in the data.

We assessed the statistical significance of the MES in two ways. Randomized MESs were generated numerically by permuting the genomic positions of our sequence reads. The background MES (MESbg) was experimentally obtained by sequencing the normal genome without affinity purification. As expected, the real data yielded markedly higher enrichment scores (Additional file 2: Figure S2). Notably, MESbg was higher than MES from randomized genomes, an indication that background sequences alone can create enrichment, probably due to chromatin accessibility and amplification bias. Consistent with recent reports [35], this illustrates the need for a proper calibration for inherent sequencing bias. Therefore, we normalized our MES with MESbg.

To find the optimal condition for normalization, we compared the statistical fitness of various normalization methods. Tag distribution along the genome can be modeled by the Poisson distribution [36, 37]. The goodness of fit was tested using the Kolmogorov-Smirnov test. In this test, a low D statistic indicates a good fit. While the Poisson model outperformed the Gaussian overall, the MES showed a better fit than raw read counts (Additional file 2: Figure S3), illustrating the rare event nature of the log-scaled read count measure. The normalized MESl calibrated by control sequencing (MESbg) yielded even better results than normalized MESt calibrated by control sequencing (MESbg).

Global and chromosomal views of DNA methylation

After confirming the best method for scoring genome-wide methylation levels at 50-bp intervals, we first examined the chromosomal methylation patterns of normal samples. The average MES calculated for each chromosome suggested that CpG-rich and gene-rich chromosomes tend to be highly methylated (Figure 1A). The methylation levels of chromosomes with large amounts of long interspersed nuclear elements (LINEs) were relatively low (e.g., chromosome 4). Interestingly, the quantity of short interspersed nuclear elements (SINEs) was proportional to the chromosome methylation pattern. This is likely caused by the fact that SINEs are typically clustered in gene-rich regions. Sex chromosomes were globally hypomethylated with lower CpG density and higher repeat content than autosomes. Since we used the tissues taken from a male in this experiment, the global hypomethylation of the X chromosome observed is not associated with X inactivation. Chromosome-wide views recapitulated high CpG density and high methylation around gene-rich (see black bars at the bottom) and CGI-rich regions (see blue bars at the top) (Figure 1B). In contrast, low CpG density and low methylation were observed around gene-poor regions that were rich in long-range repeats (> 1 kb) (see red bars at the top). The average MES suggests that the methylation level of CGIs is considerably higher than that of genic regions or repeats (Figure 1B).

Generally, CGIs tend to remain methylation free in normal tissue. To analyze the high methylation patterns of CGIs, we checked the average MES distribution and found a slightly bimodal pattern (Figure 1C). About 66% (11,376/17,284) of CGIs in the left peak overlapped with a promoter (1 kb by our definition). In contrast, 13% (1,386/10,357) of CGIs in the right peak overlapped with a promoter, suggesting that most promoter-associated CGIs are unmethylated. In contrast to promoter-related CGIs, promoter-independent CGIs were heavily methylated (Figure 1C). Although most CGI-positive promoters were not methylated, CGI-negative promoters showed relatively high methylation levels (Figure 1D). We also checked the methylation level of promoters by CpG density as previously defined [38] (Additional file 3: Table S4). The methylation pattern of promoters was inversely related to CpG density (Additional file 2: Figure S4). On the other hand, CGI-containing gene bodies had higher methylation levels than those without CGIs (Figure 1D).

Changes in DNA methylation patterns associated with gastric cancer

When the average chromosomal MES of the cancer methylome was compared to that of control tissue, we found that all chromosomes in the cancer tissue tended to be hypomethylated (Additional file 2: Figure S5). With chromosome-wide views, CGI-rich regions were found to be specifically hypermethylated, while repeat-rich regions were widely hypomethylated (Figure 2A; Additional file 2: Figure S6). To analyze methylation changes in genomic elements, we aligned each element at the start and end sites and then obtained the average MES at each respective position. Strikingly, we detected hypermethylation in the upstream region, particularly from 500 bp upstream to the transcription start site (Figure 2B). This is in accordance with the hypermethylation of promoter regions frequently observed in cancer.

The region centered at the transcription start site showed completely different patterns depending on the presence of a CGI, reflecting the low methylation status of CGI-containing promoters (Figure 2C). We also found that, in cancerous tissue, remarkable hypermethylation of CGI-containing promoters occurs and that the density of CpGs is crucial for the increase in DNA methylation (Figure 2C). To further analyze whether 5' regions of genes were hypermethylated similarly to gene promoters, we checked the methylation pattern of the first exons. Interestingly, we found that the first exon was hypermethylated only when it was the 5'-end of a coding exon, but not when it was a 5' UTR exon (Figure 2D). These regions also contained high CpG density. Therefore, CGIs at the upstream regions of genes, the promoter, and the coding start appear to be the major targets of DNA hypermethylation in cancer.

Methylation pattern of CpG islands

To explore the correlation between the location of CGIs and DNA methylation, we subgrouped CGIs according to their position within the genome. Specifically, they were categorized as 5' (located between 1 kb upstream and the coding start site of a gene), intragenic (intragenic CGIs outside the 5' end), and intergenic (located in non-genic region) (Additional file 1: Table S5). Although CpG density was similar among the three groups, non-5' CGIs (intragenic and intergenic CGIs) were significantly more methylated than 5' CGIs (Additional file 2: Figure S7). We further compared the average MES of subgrouped CGIs and found the methylation of all CGIs was generally increased. However, the relative differential MES suggested that the change in methylation in 5' CGIs was significantly greater than that for other CGIs (Figure 3A), reflecting the important roles of 5' CGIs in cancer. The extent of 5' CGI hypermethylation significantly correlated with the overlap of the transcription start site (Figure 3B).

DNA methylation of repetitive elements

DNA hypomethylation of repetitive elements is a common feature of cancer. The methylation status of individual repetitive elements such as Alu, LINE-1, and Sat2 is a major target of global DNA methylation studies. However, a comprehensive view of the methylation status of all repetitive elements is currently unavailable. We first checked the methylation status of various repetitive elements to verify the changes in each repetitive element in cancer. When we compared the DNA methylation levels in normal and cancerous tissue, striking changes were observed in the distribution of repeat methylation in the cancerous tissue (Figure 4A).

In addition to the hypomethylation patterns of the SINE and LINE, a significant reduction in methylation was found in SVA (SINE-VNTR-Alus), satellites, and LTR. This is consistent with previous reports about cancer-specific hypomethylation [41]. Alu elements are the most abundant class of repetitive elements in the human genome, with these elements having over one million copies and spanning over 30 lineages. AluS and AluY elements, which are younger subfamilies, were significantly hypomethylated when compared with older subfamilies, as previously reported [42] (Figure 4B). SVA elements, which have been extensively mobilized in the human genome, consist of a combination of sequences derived from other retroelements [43].

To understand if a correlation may exist between the degree of methylation changes and genomic location, we analyzed the methylation of repetitive elements based on their genomic location (i.e., promoter, gene body, or intergenic), even though it is unclear if repetitive elements participate in regulating gene expression (Figure 4C). We found that methylation changes were higher in gene-associated repetitive elements than in gene-independent repetitive elements. This suggests that a correlation exists between methylation changes in repetitive elements and the expression of adjacent genes.

Long-range epigenetic silencing (LRES) in gastric cancer

Large chromosomal regions can usually be suppressed in cancer cells, as seen by hypermethylation of neighboring CpG islands and downregulation of most genes within the region. To determine whether LRES occurs in gastric cancer, we identified large-scale genomic regions where methylation enrichment encompasses multiple genes. Interestingly, we found a high degree of hypermethylation around the 12q14 site, even though there were few CGIs and genes with abundant repeats (Figure 5A).

The LRES region around the 12q14 site spanned about 2.7 Mb and harbored about 21 genes. Among these genes was MDM2, which encodes a protein that is considered to be a negative regulator of p53 and a major regulator of cancer development. Promoters of genes in this LRES region displayed tumor-specific hypermethylation (p = 0.03, t-test) (Figure 5B). To determine whether genes in this LRES region show concordant gene silencing, we re-analyzed publicly available expression datasets (GSE27342) for differential gene expression in clinical samples. Consistent with previously reported expression patterns of LRES regions, the public data for this LRES region showed common gene suppression (p = 0.05, t-test) (Figure 5C).

To determine if the hypermethylation pattern of genes within this LRES region is commonly present in gastric cancer, we examined the methylation enrichment frequency of selected genes (Figure 5D). The methylation enrichment levels of cancer samples were over 2-fold higher than that of normal samples. However, one of the patients showed a low methylation level at several target sites (Additional file 2: Figure S9). Because the amplification of this region frequently occurs in many cancers, we examined the amplification frequency of several genes within these regions using real-time PCR. Intriguingly, we detected amplification of MDM2 (Additional file 2: Figure S10), suggesting that an interaction exists between DNA methylation and gene amplification. To examine the generality of MDM2 methylation, we used pyrosequencing to analyze the methylation level of a specific locus in an upstream region of MDM2. We analyzed MDM2 methylation in normal and cancer tissue samples from other 28 gastric cancer patients. Out of 28 independent samples, most patients showed higher MDM2 methylation (Figure 5E), while four samples showed decreased levels of methylation, along with gene amplification (Additional file 2: Figure S11). Therefore, the MDM2-containing region appears to be hypermethylated in cancer in general, but our results suggest that gene amplification in this region interferes with methylation. Therefore, except for cases with gene amplification, LRES across 12q14 appears to be a distinct epigenetic pattern associated with gastric cancer.

Methylation patterns in metastatic lymph nodes

In gastric cancer, lymph node metastasis is one of the major prognostic factors and an important indicator of tumor aggressiveness. Many studies have been conducted to analyze the expression profiles for and epigenetic changes in specific genes, but comprehensive data for whole-genome alterations in DNA methylation remain scarce. Similar to cancerous tissue, metastatic lymph nodes showed hypermethylation at the region centered at the transcription start site only in CGI-containing promoters (Figure 6). GO analysis revealed that several clustered groups of genes with highly differential CGI MES (differential MES > 1) in 5' CGIs overlapped with cancer genes, such as HOX family and histone family genes, even if there was little difference in ranks (Additional file 1: Table S8). After considering all these data, we conclude that the methylation pattern of metastatic lymph nodes is considerably similar to that of the cancer of origin.