Genome-wide methylated CpG island profiles of melanoma cells reveal a melanoma coregulation network

Abstract

Metastatic melanoma is a malignant cancer with generally poor prognosis, with no targeted chemotherapy. To identify epigenetic changes related to melanoma, we have determined genome-wide methylated CpG island distributions by next-generation sequencing. Melanoma chromosomes tend to be differentially methylated over short CpG island tracts. CpG islands in the upstream regulatory regions of many coding and noncoding RNA genes, including, for example, TERC, which encodes the telomerase RNA, exhibit extensive hypermethylation, whereas several repeated elements, such as LINE 2, and several LTR elements, are hypomethylated in advanced stage melanoma cell lines. By using CpG island demethylation profiles, and by integrating these data with RNA-seq data obtained from melanoma cells, we have identified a co-expression network of differentially methylated genes with significance for cancer related functions. Focused assays of melanoma patient tissue samples for CpG island methylation near the noncoding RNA gene SNORD-10 demonstrated high specificity.

Figure 1. Characterization of genome-wide methylation patterns in melanoma cell lines.

(A) The overall workflow, and results of next-generation sequencing and analysis of MBD2-pulldown DNA. (B) Frequency distributions of sizes of the MBD2-pulldown DNA across five cell lines. The methylated regions were defined by MACS peak calling of the MBD2-pulldown DNA against input (non-enriched) DNA. (C) & (D) Venn diagram of the number of MBD2-pulldown regions shared or unique among the various cell lines. (E) Venn diagram describing the number of MBD2-pulldown DNA regions partitioned among various cell lines, which defines a unique set of stage IV specific MBD2-pulldown DNA regions. (F) Genome-wide distribution of differential methylated regions (as defined as the MBD2-pulldown DNA). The outermost circle displays the human chromosomes. The inner circles represent the genome-wide distribution of MBD2-pulldown DNA. The height of the histogram bins indicates the density of methylated regions. From the inner most circle, which represents the distribution of MBD2-pulldown DNA in 5AzadC treated WM1552C (black), the circles from inside going outside represent HEM-1 (orange), WM793B (blue), WM1552C (green), A375 (red) and SKEML2 (red).

Figure 2. DNA fragments pulled-down with MBD2 are enriched for methylated CpG islands.

Results of bisulfite sequencing of MBD2-pulldown DNA near five genes are shown. Methylated CpG sequences are shown as closed boxes (unmethylated sequences are open boxes). The number at the right of each row corresponds to the number of CpG in the region that was methylated, with the total number of CpG represented as the total number of boxes in each row. Each row is an independent experiment.

Figure 3. Functional significance of methylated CpG islands in melanoma cells.

(A) The functional class distribution of CpG methylation in the five cell lines. Gene annotation was taken from Refseq at UCSC genome browser. Promoter region was defined as 3000 bp upstream of the transcription starting site (TSS), and the downstream region is defined as 3000 bp downstream of translation stop sequence. (B) Heatmap of methylation fold change in ± 1000 bp windows centered at TSS for cancer related genes in the five cell lines. The rows represent 3,688 cancer-related genes downloaded from the CancerGene database by Higgin et al., which were filtered for significant methyl CpG levels in melanoma lines relative to HEM-1. Color scale on the right represents MBD2-enrichment, i.e., methyl CpG, fold change. (C) Functional enrichment analysis of diseases and disorders of the gene set in (B).

Figure 4. CpG methylation status near SNOR-10 gene in different patient samples.

(A) depicts the results obtained with next-generation sequencing of MBD2-enriched DNA in the five cell lines, plus those in WM1552C cells treated with 5AzadC and those in the input DNA. (B) A region between −249 bp and −31 bp of the TSS of SNORD-10 is shown magnified. Methylation status of this region in the DNA obtained from various metastatic melanoma patient samples or in non-melanoma nevi was determined by bisulfite sequencing (see Methods). Each box is a CpG; open boxes are unmethylated, filled boxes are methyl CpG sequences. Numbers at the right represent the number of methyl CpG detected in the region corresponding to each row (independent sequence analysis).

Figure 5. MBD2-enriched DNA from WM1552C cells defines a large co-expression network with cancer and metastasis-related genes.

A co-expression network formed by the genes that exhibit (a) enrichment for MBD2-pulldown relative to HEM-l, (b) differential RNA expression levels measured by RNA-seq relative to HEM-l, and (c) included in cancer-related genes by Higgins et al. (A)–(E) Proteins of the connected component in the co-repression network (Fig. S4) that are known to interact as members of known protein complexes. (F)–(J) Protein complexes constructed by extending the core protein interaction networks in (A–E) as seeds by incorporating known first-degree interaction neighbors of the member proteins. The protein modules are enriched for RNA-dependent DNA replication (F), ligand-dependent transcription factor activity (G), DNA damage response and cell death (H), basal lamina/substrate-dependent cell movement (I), and actin filament sliding/cell movement (J) (see text for P-values).

Figure 6. Patterns of methyl CpG at repetitive elements enriched in the five cell lines.

The relative levels of MBD2-enrichment at repeated DNA elements compared to input (un-enriched) DNA are plotted in a color coding scheme in which the color represents the enrichment significance P-values.