DNA methylation screening identifies driver epigenetic events of cancer cell survival

Abstract

Cancer cells typically exhibit aberrant DNA methylation patterns that can drive malignant transformation. Whether cancer cells are dependent on these abnormal epigenetic modifications remains elusive. We used experimental and bioinformatic approaches to unveil genomic regions that require DNA methylation for survival of cancer cells. First, we surveyed the residual DNA methylation profiles in cancer cells with highly impaired DNA methyltransferases. Then, we clustered these profiles according to their DNA methylation status in primary normal and tumor tissues. Finally, we used gene expression meta-analysis to identify regions that are dependent on DNA methylation-mediated gene silencing. We further showed experimentally that these genes must be silenced by DNA methylation for cancer cell survival, suggesting these are key epigenetic events associated with tumorigenesis.

Figure 1. Clustering of DNMT deficient cells identifies three classes of putative driver genes marked by DNA methylation. See also Figure S1 and Table S1

(A) One-dimensional hierarchical clustering using Euclidean distance and average linkage was performed with the ~24,000 Infinium DNA methylation probes located outside of repeats or known SNPs in HCT116 wild type, DKO8 and DKO1 cell lines. Each row represents a probe; each column represents a sample. The beta value (level of DNA methylation) for each probe is represented with a color scale as shown in the key. (B) K-means (K=4) clustering of the 566 Infinium DNA methylation probes that maintain DNA methylation in DKO1 sample (a beta value of at least 0.6 and a difference between HCT116 and DKO1 smaller than 0.2) in (A) for 10 TCGA samples (n=4 normal colon and n=6 primary colon adenocarcinoma). (C) Heatmap of 566 infinium DNA methylation probes in 32 normal tissues retaining the probe order from (B). Primary normal bladder (n=4), sperm (n=1), and primary normal TCGA kidney (n=15), lung (n=4) and ovary (n=8). Whole Genome Amplified DNA (WGA) was used as a negative control for DNA methylation. (D) Bisulfite sequencing validation of Infinium DNA methylation data from two regions (LDHAL6B and ADAM2) from the somatic-specific DNA methylation cluster and two regions (ARMCX1 and MEOX2) from the cancer-specific DNA methylation cluster. Arrow indicates the position of the Infinium probe. Empty and filled circles denote unmethylated and methylated CpG sites, respectively. Each horizontal row represents one sequenced DNA clone. The number on the right represents the mean DNA methylation score of each region and the number in the parentheses represents the mean DNA methylation score of the specific Infinium CpG site.

Figure 2. Residual methylation in DKO1 is not caused by an inherent susceptibility to DNA methylation

(A) Validation of H3K27me3 status in ESC as a predictive method for DNA methylation in HCT116 cells. Methylation status of ~27,000 CpG sites was determined by Infinium. T-test with Mann Whitney post-test. Data represent the mean ± SEM. (B) Frequency of probes marked by H3K27me3 in ESC in the cohort of DNA methylated probes (Beta value >0.6) in HCT116, DKO8 and DKO1 cells. (C) Validation of the predictive method based on genomic architecture (Estecio et al., 2010) in HCT116 cells. Methylation status of ~27,000 CpG sites was determined by Infinium. One way ANOVA with Kruskal-Wallis test. Data represent the mean ± SEM. (D) Frequency of methylation-prone genes in the cohort of DNA methylated genes (Beta value >0.6) in HCT116, DKO8 and DKO1 cells.

Figure 3. Validation of CpG sites identified with cancer-specific DNA methylation using independent datasets and association with gene repression

(A) Volcano plot of the CpG loci identified as cancer-specifically methylated in colon adenocarcinoma (Normal n=16, Cancer n=168) from the TCGA Data Portal. The beta value difference in DNA methylation between the tumor samples and the correspondent normal samples is plotted on the x-axis, and the p value for a FDR-corrected Wilcoxon rank sum test of differences between the tumor and correspondent normal samples (-1* log10 scale) is plotted on the y-axis. Probes that are significantly hypermethylated (FDR adjusted p<0.05) in tumors are shown in red. (B) Volcano plot gene expression data of cancer-specific DNA methylated genes. Gene expression data was obtained from GEO (GSE 8671) from primary normal colon (n=32) and primary colon cancer (n=25). For the volcano plots, gene expression fold change between the normal tissues and the tumor tissues is plotted on the x-axis, and the p value for a FDR-corrected t-test of differences between the normal and the tumor tissues (-1* log10 scale) is plotted on the y-axis. Probes that are significantly (p<0.05) down-regulated in tumor tissues are shown in red. (C) Volcano plot of the CpG loci identified as cancer-specifically methylated in lung adenocarcinoma (Normal n=4, Cancer n=19) from the TCGA Data Portal. The beta value difference in DNA methylation between the tumor samples and the correspondent normal samples is plotted on the x-axis, and the p value for a FDR-corrected Wilcoxon rank sum test of differences between the tumor and correspondent normal samples (-1* log10 scale) is plotted on the y-axis. Probes that are significantly hypermethylated (FDR adjusted p<0.05) in tumors are shown in red. (D) Volcano plot gene expression data of cancer-specific DNA methylated genes. Gene expression data was obtained from GEO (GSE7670) from primary lung adenocarcinoma (n=27) and primary lung (n=30). For the volcano plots, gene expression fold change between the normal tissues and the tumor tissues is plotted on the x-axis, and the p value for a FDR-corrected t-test of differences between the normal and the tumor tissues (-1* log10 scale) is plotted on the y-axis. Probes that are significantly (p<0.05) down-regulated in tumor tissues are shown in red. (E) Venn diagram showing the overlap between the genes statistically hypermethylated in colon adenocarcinoma (n=50, FDR adjusted p<0.05) and lung agenocarcinoma (n=33, FDR adjusted p<0.05). (F) Venn diagram showing the overlap between the genes statistically repressed in colon adenocarcinoma (n=44, FDR adjusted p<0.05) and lung agenocarcinoma (n=25, FDR adjusted p<0.05).

Figure 4. Apoptosis analysis of HCT116 and DKO1 cells. See also Figure S2

(A) HCT116 wild type and HCT116 DKO1 cells were stained with annexin V-FITC and Propidium Iodide (PI) and analyzed by FACS, showing an increased level of basal apoptotic cell death in the HCT116 DKO1 cell line compared to HCT116 wild type. HCT116 DKO1 cells were then sorted in viable (annexin V and PI negative) and early apoptosis (annexin V positive and PI negative). (B) The morphology of viable DKO1 and early apoptosis is clearly distinct. The apoptotic cells (blue) show a characteristic phenotype of higher SSC and lower FSC than the viable cells (red). (C) Bisulfite sequencing analysis of CpG methylation status of four regions, from cancer-specific methylated cluster (EYA4 and IRAK3) and from somatic tissue-specific DNA methylation cluster (SYCP3 and ADAM2). The mean percent methylation at each CpG site is derived from clones showed on Figure S2A. The capped line represents the region analyzed.

Figure 5. Functional validations. See also Figure S3

(A) Overexpression of nine candidate genes from the cancer cluster (P2RY14, IRAK3, CDO1, ESX1, ARMCX1, BCHE), somatic cluster (SYCP3 and ADAM2), cell culture cluster (STEAP4). Shown is the fraction of Empty-Vector at the indicated times, normalized to the day 0 values. NOX4 was used as a control gene, since it is hypermethylated in HCT116 cells and completely demethylated in DKO1 cells. (B) Overexpression of the same nine candidate genes reduces viability of RKO cancer cells. Shown is the fraction of Empty-Vector at the indicated times, normalized to the day 0 values. NOX4 was used as a control gene. (C) Meta-analysis using the oncomine (www.oncomine.org) for IRAK3 expression. Box plots showing decreased expression of IRAK3 during tumorigenesis on datasets performed in colon adenocarcinoma (Kaiser et al., 2007); lung adenocarcinoma (Su et al., 2007); prostate carcinoma (Welsh et al., 2001) and cutaneous melanoma (Talantov et al., 2005). The y-axis represents log₂ median-centered intensity (normalized expression). Shaded boxes represent the interquartile range (25th–75th percentile). Whiskers represent the 10th–90th percentile. The bars denote the median. (D) Overexpression of IRAK3 in HCT116 cells induces a reduction in the Survivin levels. Western-Blot analyzes of IRAK3 and Survivin after lentiviral infection with pLJM1 empty vector (E/V) or pLJM1-IRAK3. Histone H3 was used as a loading control. (E) IRAK3 expression induces cell death of cancer cells. HCT116 infected with pLJM1 empty vector or pLJM1 IRAK3 were stained with annexin V-FITC and Propidium Iodide (PI) and analyzed by FACS, showing and increased level of cell death in the cell overexpressing IRAK3 (upper panel). Re-expression of IRAK3 in HCT116 wild type cells showed a reduced cell number in culture than HCT116 empty vector (lower panel). *P<0.05, ***P<0.0001. Data represent the mean ± SEM. (F) IRAK3 knock-down induces colony formation in a non-tumorigenic cell. UROTSA infected with a shRNA against IRAK3 presented a higher colony formation activity than UROTSA infected with a scrambled shRNA. Western-Blot analysis of IRAK3 after lentiviral infection with shRNA against IRAK3 or a scrambled shRNA. Actin was used as a loading control.