Multi-faceted epigenetic dysregulation of gene expression promotes esophageal squamous cell carcinoma

Abstract

Epigenetic landscapes can shape physiologic and disease phenotypes. We used integrative, high resolution multi-omics methods to delineate the methylome landscape and characterize the oncogenic drivers of esophageal squamous cell carcinoma (ESCC). We found 98% of CpGs are hypomethylated across the ESCC genome. Hypo-methylated regions are enriched in areas with heterochromatin binding markers (H3K9me3, H3K27me3), while hyper-methylated regions are enriched in polycomb repressive complex (EZH2/SUZ12) recognizing regions. Altered methylation in promoters, enhancers, and gene bodies, as well as in polycomb repressive complex occupancy and CTCF binding sites are associated with cancer-specific gene dysregulation. Epigenetic-mediated activation of non-canonical WNT/β-catenin/MMP signaling and a YY1/lncRNA ESCCAL-1/ribosomal protein network are uncovered and validated as potential novel ESCC driver alterations. This study advances our understanding of how epigenetic landscapes shape cancer pathogenesis and provides a resource for biomarker and target discovery.

Conflict of interest statement

E.C. is consultant at Takeda, Merck, Loxo, and Pear Diagnostics, reports receiving commercial research grants from AstraZeneca, Ferro Therapeutics, Senti Biosciences, Merck KgA and Bayerand stock ownership of Tatara Therapeutics, Clara Health, BloodQ Guardant Health, Illumina, Pacific Biosciences and Exact Biosciences. T.G.B. is an advisor to Revolution Medicine, Novartis, Astrazeneca, Takeda, Springworks, Jazz, and Array Biopharma, and receives research funding from Revolution Medicine and Novartis. M.S. is Cofounder and scientific advisory board member of Personalis, SensOmics, Mirvie, Qbio, January, Filtricine, and Genome Heart. He serves on the scientific advisory board of these companies and Genapsys and Jupiter. Other authors declare no competing interests.

Figures

Fig. 1. Epigenetic landscape and heterogeneity in esophageal squamous cell carcinoma (ESCC).

a Ten pairs of ESCC and adjacent normal tissues were performed whole-genome bisulfite sequencing (WGBS). The asymmetric density distribution of all CpG methylation statuses in the normal esophageal tissues versus ESCC. ESCCs lose methylation which leaves most CpGs partially methylated. Normal = blue, tumor = red. b Circos plot of >5 million differentially methylated CpGs (DMCs) between ESCC tumor and adjacent normal tissue. DMCs are substantially hypomethylated in ESCC (97.3%). Only 2.7% are hypermethylated in ESCC. c Principal component analysis (PCA) shows that characteristic CpGs discriminate tumor samples from normal samples. d t-Distributed Stochastic Neighbor Embedding (t-SNE) showed CpG methylation profiling of TCGA-esophageal cancer from human methylation 450K analysis clustered into either normal tissue (n = 15, green circles) or ESCC (n = 97, red circles) or esophageal adenocarcinoma (n = 89, black circles) subtypes. e Entropy analysis of all CpGs showed variations per CpG in normal esophageal tissues (blue bars) and ESCC (red bars). The entropy of CpGs in ESCC was higher than in normal samples. f Multivariate cox proportional hazard analysis demonstrated TCGA-ESCC patients (n = 92) with lower variance of CpG methylation in tumors showed better survival time than those with higher variance. Median variance was used to discriminate high versus low variance groups. Variance of DNA methylation changes were normalized for age, gender, and alcohol consumption in patients. Statistical significance was assessed by two-sided Wald test, p = 0.002.

Fig. 2. Differentially methylated regions (DMRs) and their functional impacts on the ESCC genome.

a A DMR identification algorithm from DMC was developed using two criteria: (1) two flanking DMCs should be close (150 base pairs, bp) given the minimum size of CpG island (CGI) 150 bp and (2) the methylation pattern should be consistent, either hypomethylated or hypermethylated within a DMR. Our algorithm revealed the distribution of DMR size and CpG density within DMRs. Both distributions of DMR size and CpG density are asymmetric and have long tails as DMR size increases in length and CpG density is more compact. The peak of DMR size is 200–300 bp and the peak of CpG density is approximately 4% (also seen in Supplementary Fig. 8). b Methylation level of CpGs within 15,000 bp upstream and downstream relative to a transcription start site (TSS) was assessed in ESCC and normal esophagus separately. Overall methylation conversions between normal esophagus tissue and ESCC were observed. CpGs in ESCC tend to be hypermethylated between 2990 bp upstream and 6990 bp downstream of a given TSS. The arrow indicates the p value for the specific region of significant methylation changes. c A representative genomic region at chr8:55,360,000–55,400,000 with hypermethylation in CpG island and hypomethylation in the CpG shore. d SOX17 expression is decreased in ESCC tumors (T, 10 cases) relative to normal tissue (N, 10 cases). Y-axis is the normalized gene expression levels (transcript per million reads, TPM) from RNAseq data. Box and whisker plot: center line, median; box limits, upper and lower quartiles; and whiskers, maximum and minimum values, circle dots, outliers. Statistical significance was assessed by two-sided t-test, p = 0.67, FDR = 0.8. e Genomic regional enrichment analysis of ESCC-associated DMRs. In all, 289,973 hypo-DMRs and 5322 hyper-DMRs were mapped to LOLA core database with ENCODE, DNase, CODEX, and UCSC genomic annotations. The significant overlapping genomic regions were selected with p value <0.05 and odds ratio >2 from Fisher’s exact test. Each bar with pseudocolor gradience represents each dataset in LOLA core databases (see Supplementary Data File 3).

Fig. 3. Integrative analysis of WGBS and RNAseq uncovered methylation-mediated diverse gene regulation.

a A total of 694 genes were selected for the methylome–transcriptome association analysis. The selected genes have statistically significant DMRs (FDR ≤ 0.001) in promoters, defined as 4500 base pair (bp) upstream and 500 bp downstream relative to transcription start sites, and are statistically significant DEGs (differentially expressed genes) (|log2(fold change)| > 1, FDR ≤ 0.05). b, c The association of promoter methylation and expression of the 694 genes were identified. There are four clusters: C1 (n = 48): genes that are hypermethylated in the promoter with low expression level in ESCC; C2 (n = 389): genes that are hypomethylated in the promoter with high expression in ESCC; C3 (n = 67): genes that are hypermethylated in the promoter with high expression in ESCC; C4 (n = 190): genes that are hypomethylated in the promoter with low expression in ESCC. Genes in C1 and C2 fit the canonical model of regulation, while genes in C3 and C4 are not well explained by current understanding. Representative genes are listed in each cluster. See Supplementary Data File 4. d The quantification of CpG methylation in gene promoters and gene bodies in C3 (n = 67) is significantly higher than in C1 (n = 48), p < 0.001, while no significant difference (NS) for enhancers. e CTCF-binding sites are significantly higher in C3, indicating hypermethylation of inhibitors leads to de-repression to promote gene expression. The sample sizes for C1, C2, C3, and C4 in d and e are the same as defined in b and c. Box and whisker plot: center line, median or mean; box limits, upper and lower quartiles; and whiskers, maximum and minimum values, circle dots, outliers. Statistic significance was assessed by one-way ANOVA.

Fig. 4. DNA methylation at regulatory consensus protein-binding sites and impacts on gene expression.

a Functional annotation for the four distinct methylation-transcriptome clusters. Lists of genes in C2 (n = 389), C3 (n = 67), and C4 (n = 190) were subjected to gene set analysis using hypergeometric statistics for gene sets collected from multiple databases (ENCODE, CHEA, KEGG, WikiPathways, Reactome, GO molecular function, Panther, BIOGRID, etc.). The significance of the hypergeometric analysis is indicated as –Log10 (p value) in the form of a horizontal histogram where bar heights represent level of significance. Bars are color coded based on their inclusion in each cluster. Gene pathways or GO terms in different clusters, including polycomb repression complex 2 (PRC2) subunit, EZH2 (Ester of Zinc Finger Homolog 2), and SUZ12 (Polycomb Repressive Complex 2 Subunit)-binding sites significantly enriched in C3, suggesting hypermethylation of PRC2 de-represses gene expression. Genes in C1 have no significant gene set enrichment. b, c In silico analysis of EZH2 or SUZ12-binding sites within gene promoters in each cluster (C1–C4) and show more significant binding scores in C3 than in other groups, Box and whisker plot: center line, median (red) or mean (green); box limits, upper and lower quartiles; and whiskers, maximum and minimum values; red dots, outliers. Statistic significance was assessed by one-way ANOVA. d The probability of the top 20 transcription factor consensus-binding sites in C3 showed EZH2 has the highest binding scores in a subset of genes in C3, including WNT2. The heatmaps for C1, C2, and C4 are in Supplementary Fig. 18. e In non-canonical gene cluster C3 (n = 67), WNT2 is significantly hypermethylated in the promoter (one-way ANOVA, FDR = 6.6005e−03) and highly expressed in ESCC (one-way ANOVA, FDR = 0.0039). WNT2 also shows the highest fold change in gene expression in ESCC relative to adjacent normal tissues. Each dots represent individual gene; colored dots reflect different categories genes.

Fig. 5. Hypermethylation in the WNT2 promoter leads to high WNT2 expression in ESCC.

a Chromatin immunoprecipitation by high-throughput DNA sequencing (ChIP-seq) showed EZH2 preferentially binding to the hypomethylated WNT2 promoter region in normal cells (Het-1A) relative to the hypermethylated WNT2 promoter region in ESCC cells (EC109). b The WNT2 promoter region is hypermethylated in ESCC cell lines by methylation-specific PCR (MS-PCR) analysis. M methylation detection, U unmethylation detection, PC positive control, NC negative control. Het-1A cells are an immortalized normal esophageal epithelial cell line. EC9706, EC109, and EC1 are patient-derived ESCC cell lines. c Western blot analysis showed WNT2 protein is overexpressed in three ESCC cancer cell lines compared to normal cell line Het-1A. d EC9706 cancer cells were treated with or without the DNA methyltransferase inhibitor, 5-Azacytidine (5-AzaC), for 24 h. The methylation change status in promoter region of WNT2 was detected by MS-PCR. e ChIP-PCR detection of EZH2 pull-down in EC9706 cancer cells were treated with or without 5-AzaC for 24 h. f RT-PCR detection of WNT2 expression in EC9706 cancer cells were treated with or without 5-AzaC for 24 h. g WNT2 is overexpressed in ESCC tumor (n = 10) and adjacent normal samples (n = 10). Left panel: protein expression from western blots, right panel: RNA abundance (FPKM) from RNAseq. h WNT2, MMP3, and MMP9 expression in the tumor and normal samples. Representative blots from 10 paired normal and ESCC tumor samples. i WNT2 depletion with two independent siRNAs inhibited MMP3 and MMP9 expression in the ESCC cell line EC9706. j Two independent siRNAs to silence WNT2 expression reduced the invasion and migration of ESCC cells relative to siRNA controls. Scale bar:100 µm. k Xenograft tumor growth curve and tumor weight in EC9706 cells expressing either shRNA control (sh-NC, n = 5) or sh-Wnt2 (n = 5). l Schematic representation of the mechanism of EZH2/PRC2-WNT2-MMP signaling upregulation in ESCC. Functional study of WNT2 in EC109 cell line is presented in Supplementary Fig. 27. Representative results from three independent experiments are presented in b–e, j; triplicates in each condition are shown in f. The bar plots are plotted as mean ± s.d. Statistic significance was measured using two-sided t-test. *p < 0.05, **p < 0.01. ***p < 0.001, ****p < 0.0001. Source data are provided as a Source Data file.

Fig. 6. Hypomethylation-mediated upregulation of long non-coding RNAs (lncRNAs) in ESCC.

a In canonical gene cluster C2 (n = 389), ESCCAL-1 is significantly hypomethylated in the promoter (one-way ANOVA, FDR = 1.7386e−04) and highly expressed in ESCC (one-way ANOVA, FDR = 0.01). ESCCAL-1 also shows the most significant and substantial methylation difference between normal esophageal tissues and ESCC among the lncRNAs in C2. Each dots represent individual gene; colored dots reflect different categories genes. b Loss of CpG methylation at the ESCCAL-1 promoter region in ESCC (chr8:76,135,639–76,236,976 of GRCh37/hg19). WGS (whole-genome sequencing) of three ESCC patients shows no mutation or copy number variations detected at the above indicated region. This is validated by TCGA ESCA data (n = 186), where no mutation or copy number variations were observed (see Supplementary Fig. 24d). DMRs around transcription start sites (TSSs) of the two isoforms showed extensive differentiation between ESCC and normal samples. c ESCCAL-1 was significantly differentially expressed and highly abundant in ESCC samples (n = 10) relative to normal samples (n = 10). Statistic significance was assessed by one-way ANOVA, p value = 0.0013, log2 (fold change) >1, FDR = 0.0044. d Methylation status of the ESCCAL-1 promoter region was verified on an independent matched normal (n = 32) and ESCC tumor (n = 32) samples using a methylation-specific PCR (MS-PCR) assay; four representative PCR results are shown. M: PCR with methylation primers, U: PCR with unmethylation primers. e Quantification of MS-PCR results in the 32 paired normal and tumor samples. Chi square analysis tested for significance between groups. P value < 0.01. f ESCCAL-1 expression is significantly higher in ESCC tumors in an independent cohort of 73 matched normal and tumor samples. Error bars are median values with 95% confidence intervals. Significance for comparison between the two cohorts was measured using an unpaired two-sided Student’s t-test, p = 0.002.

Fig. 7. Oncogenic functions of ESCCAL-1 in ESCC.

a Hypomethylation at ESCCAL-1 promoter regions was confirmed in three different ESCC cell lines using methylation-specific PCR. M: PCR with methylation primers, U: PCR with unmethylation primers. Het-1A: an immortalized esophageal epithelial cell line. EC1, EC109, and EC9706 are patient-derived ESCC cell lines. b ESCCAL-1 expression was significantly higher in three ESCC cell lines (EC1, EC109, and EC9706) relative to a normal esophageal epithelial cell line (Het-1A). Triplicates in each cell lines, mean ± s.d., unpaired two-sided t-test, **p < 0.01, ****p < 0.0001. c, d ESCC patients with higher expression of ESCCAL-1 exhibit worse overall (OS) and progression-free survival time (PFS), risk ratio = 2.56, log-rank test, p value = 0.003. e, f shRNA knockdown of ESCCAL-1 inhibited tumor growth in a tumor xenograft mouse model (n = 6 in each group, unpaired two-sided t-test, p value = 0.001). ESCCAL-1 remains greater than 50% lower expression in xenograft tumors of ESCCAL-1 knockdown group relative to the control group. Bar plots (mean ± s.d.) indicate triplicates in each condition for RT-PCR. **p < 0.01 was measured by unpaired two-sided t-test. g Chromatin immunoprecipitation by YY1 transcription factor protein-directed antibody followed by a standard polymerase chain reaction (PCR) assay in ESCC cancer cells (EC109) and normal esophageal cells (Het-1A). IgG was used as negative control. H3 was used as a positive control. Representative results from three independent experiments. h Representative western blot assay showed decreased YY1 protein expression following 72 h post-transfection of three independent siRNAs targeting various transcript regions of YY1 in ESCC cell line (EC109). ESCCAL-1 expression was measured in YY1 knockdown cells using RT-PCR. Bar plots (mean ± s.d.) showed quantification of YY1 protein expression from two independent experiments or ESCCAL-1 expression from triplicates in each condition. Statistical significance was assessed by unpaired two-sided t-test, *p < 0.05, **p < 0.01. i RNAseq was performed in duplicate in the normal esophageal cell line Het-1A, ESCC cancer cells EC9706 with control shRNA, or EC9706 with an shRNA against ESCCAL-1. Unsupervised hierarchical clustering of differential gene clusters between the three conditions is shown. Differentially expressed genes were selected based on an iterative clustering approach selecting for genes with the top 5% of the most variable and differential gene expression. Two hundred and ten genes were identified and subjected to functional annotation using a hypergeometric test in multiple databases. j Gene Set Enrichment Analysis (GSEA) for RNA-seq data from either ESCCAL- 1 knockdown or Myc knockdown in EC9706 showed ribosomal genes (pathways) enriched in EC9706 cells. k Diagram illustrating YY-1 binding to hypomethylated promoter regions of ESCCAL-1, driving its overexpression and leading to dysregulation of ribosomal genes and ESCC progression. Source data are provided as a Source Data file.