Oncogenic enhancers drive esophageal squamous cell carcinogenesis and metastasis

The role of cis-elements and their aberrations remains unclear in esophageal squamous cell carcinoma (ESCC, further abbreviated EC). Here we survey 28 H3K27ac-marked active enhancer profiles and 50 transcriptomes in primary EC, metastatic lymph node cancer (LNC), and adjacent normal (Nor) esophageal tissues. Thousands of gained or lost enhancers and hundreds of altered putative super-enhancers are identified in EC and LNC samples respectively relative to Nor, with a large number of common gained or lost enhancers. Moreover, these differential enhancers contribute to the transcriptomic aberrations in ECs and LNCs. We also reveal putative driver onco-transcription factors, depletion of which diminishes cell proliferation and migration. The administration of chemical inhibitors to suppress the predicted targets of gained super-enhances reveals HSP90AA1 and PDE4B as potential therapeutic targets for ESCC. Thus, our epigenomic profiling reveals a compendium of reprogrammed cis-regulatory elements during ESCC carcinogenesis and metastasis for uncovering promising targets for cancer treatment.

Esophageal cancer has two main subtypes, esophageal squamous cell carcinoma (ESCC, further abbreviated EC) and esophageal adenocarcinoma (EAC)1. It is associated with a poor prognosis and a 5‐year overall survival rate of 18% only2. ESCC accounts for ~90% of esophageal cancer cases worldwide and predominates in eastern Asia and Africa3. Although the 5-year survival rate of patients with esophageal cancer has improved in the past decade, this cancer still has a very poor prognosis because patients generally do not show clinical features before an advanced stage unless via screening. Stage III tumors invade through the muscular layer and involve lymph nodes or other adjacent structures3. The 5-year survival rate after esophagectomy can benefit from lymph node removal, and the depth of invasion and presence of node metastasis are important predictors4. Thus, a better understanding of molecular dependencies and vulnerabilities during esophageal carcinogenesis and its metastasis to lymph nodes is urgently needed to develop promising diagnostic and therapeutic strategies.

An aberrant and reprogrammed transcriptome is a universal hallmark of human cancers that is associated with deregulated cell proliferation, invasion, and metastasis5. The leading causes of aberrant transcriptional patterns are genetic alterations and the deregulated signaling pathways and epigenomic traits, including DNA methylation, histone modification, and other cis-regulatory elements6. Multiple studies have highlighted the importance of histone modifications in the establishment and maintenance of disease and cancer statuses7,8, and super-enhancers (SEs), comprising clustered enhancers, are considered to be the key controllers of different cell types and disease conditions9, associated with disease-specific genetic variants and mutations8,10. We and others recently proposed that key transcription factor (TF)-driven reprogramming of enhancers, which are marked by histone H3 lysine 27 acetylation (H3K27ac) and localized distal to promoters and transcription start sites (TSSs), promotes carcinogenesis and metastasis11,12,13.

Several familial genetic risk mutations in MSR1, ASCC1, CTHRC1, CRTC1, BARX1, FOXP1, TP53, CDKN2A, SMAD4, and NFE2L2, have been identified in Barrett’s esophagus and adenocarcinoma patients14,15,16,17,18. Meanwhile, some familial form of ESCC and a number of susceptibility loci in ESCC patients as well some genetic drivers (e.g., TP53) have been reported16,17,19,20,21,22, while the complexity of carcinogenesis might be more than just genetic alterations. In recent years, the epigenomic profiling of tumor-specific promoter, enhancer, and super-enhancer landscapes during carcinogenesis has been widely used to elucidate the molecular mechanisms underlying the dysregulated local and regional expression of oncogenic genes in several cancers, such as colon cancer, gastric cancer, renal cell carcinoma (RCC), and ependymoma8,23,24,25,26,27. Although global histone modifications have been linked to the evolution of distant metastasis during pancreatic and prostate cancer progression28,29, the genome-wide epigenomic landscape of ESCC tumorigenesis and metastasis is poorly understood.

In this work, to address the role of enhancer dynamics during ESCC progression and metastasis to lymph nodes, we characterize the enhancer and SE landscapes in primary ESCC and lymph node cancer (LNC) tissues as well as normal esophageal epithelium. We identify a large set of commonly altered enhancers and LNC-specific enhancers associated with regional gene expression that can be used for precise subclassification. By analyzing the trans-elements within altered enhancers, we discover several key TFs regulating esophageal cancer cell proliferation and migration. By performing an integrated analysis of SE-associated genes with drug interaction databases, we identify some promising oncotargets that are responsive to pharmacologic inhibition. These findings indicate that the identification of variant enhancers and SEs is of great help for illustrating core transcriptional regulatory circuitries and discovering therapeutic targets for esophageal and LNCs.

H3K27ac profiling defines active regulatory elements of esophageal carcinogenesis and metastasis

To characterize active cis-regulatory elements in esophageal carcinogenesis and metastasis, we generated 28 active histone mark (H3K27ac) profiles from ten paired ESCC patients by chromatin immunoprecipitation and sequencing (ChIP-seq) with freshly dissected samples of adjacent normal esophageal (Nor or Adj) tissues, primary ESCC (further abbreviated EC herein), and LNC tissues. In addition, a total of 50 transcriptomes (Supplementary Data 1) from 18 paired ESCC patients were determined by RNA sequencing (RNA-seq), and four LNC samples with few tumor cells identified in the clinicopathological examination were excluded (Fig. 1a; detailed information for patients is provided in Supplementary Fig. 1a and representative histological images of normal and cancer tissues are presented in Supplementary Fig. 1b). All tumors were stage III or above. Quality control analysis showed that the Q30 score of the sequencing reads was ~90% and that the mapping ratio to the human genome was over 75% for both ChIP-seq and RNA-seq data (Supplementary Fig. 2a–d). Before element analysis, we performed quantile normalization of the H3K27ac signals to comparable levels in 28 ChIP-seq datasets as presented in box plots (Supplementary Fig. 2e, f). In this study, active typical enhancers were defined as significant H3K27ac peaks located more than 2 kb from the nearest TSSs, and regions within ±2 kb of TSSs (TSSs ±2 kb) with significant H3K27ac occupancy were considered as active promoter elements as described in previous studies26,30. The numbers of active promoter elements and enhancers reached saturation after 24 and 19 samples, respectively, suggesting that the majority of cis-regulatory elements could be retrieved with 28 samples from ten ESCC patients (Supplementary Fig. 2g, h).

a Diagram of the experimental design for profiling the transcriptome and epigenome of adjacent normal tissues (Nor), esophageal squamous cell cancer (EC) tissues, and lymph node cancer (LNC) tissues. b Identification of the differentially distributed enhancers in Nor (n = 10), EC (n = 10), and LNC (n = 8) tissues. The numbers of altered (gained or lost) promoters (Pro) or enhancers (Enh) upon the comparison of Nor with EC or LNC, or EC with LNC are presented. Rel. relative. c Heatmap of the altered enhancers across Nor (n = 10), EC (n = 10), and LNC (n = 8) samples with H3K27ac enrichment signals. Six groups (G1–G6) of enhancers with H3K27ac enrichment signals with the indicated number of enhancer elements are presented. d Unsupervised hierarchical clustering analysis of 28 Nor, EC, and LNC samples’ differential enhancers shown in c. The number of patients number is also provided. e Principal component analysis of 28 Nor (n = 10), EC (n = 10), and LNC (n = 8) samples using all differentially distributed promoter elements shown in Supplementary Fig. 5aand enhancers shown in c. f Venn diagrams depicting the numbers of shared enhancers (left) and super-enhancers (right) across Nor, EC, and LNC samples. g Identification of super-enhancers in EC and LNC samples in infection plots. The top super-enhancer-associated genes are labeled. The number of shared, EC-specific, and LNC-specific super-enhancers across the three types of samples are displayed in pie charts. h Tracks of H3K27ac ChIP-seq and RNA-seq data at the EHF and CCND1loci in four representatives paired Nor, EC, and LNC samples. The previously identified enhancer79upstream of the CCND1 promoter is indicated. Source data are provided as a Source Data file.

Full size image

Next, we established a differential analysis pipeline to identify altered cis-regulatory elements; meanwhile, differentially expressed genes (DEGs) were identified with DESeq2 (Supplementary Fig. 3). To check the reliability of our predicted cis-elements, we first overlapped our identified H3K27ac peaks (P value <1e-09) with the Epigenomics Roadmap dataset31 and found an overlap rate of ~90% (Supplementary Fig. 4a). By combining the predicted promoter and distal enhancer elements in the 28 samples together (Supplementary Fig. 4b), we recovered 24,823 active promoter elements and 57,675 active enhancers (Supplementary Fig. 5a and Supplementary Data 2). Unsupervised hierarchical clustering of enhancer profiles easily discriminated normal samples from primary EC and LNC samples (Supplementary Fig. 4c), demonstrating a significant epigenomic difference during tumourigenesis. Referring to a previous study26, we defined gained or lost elements as those with a fold change (FC) in H3K27ac reads per kilobase per million reads (RPKM) of ≥2 and an absolute difference of ≥0.5. Considering the large variation between patients and samples, we evaluated the distribution of altered elements in tumor or LNC/normal pairs. Approximately 80% of the altered enhancer regions were significantly different (q < 0.1, paired t-test; Benjamini–Hochberg corrected) at the threshold of ≥6/10 patients (EC vs. Nor) (Supplementary Fig. 4d, e). For the simplicity and consistency of subsequent analyses, we set the same threshold for comparisons between LNC and Nor (≥6/10) or EC (≥6/8) tissues (Supplementary Fig. 4f–i). With these criteria, we obtained a high-confidence and comprehensive set of 2917 gained and 3027 lost promoter elements and 8587 gained and 9642 lost enhancers in EC tissues (relative to Nor), 1471 gained and 1403 lost promoter elements and 4399 gained and 4626 lost enhancers in LNC (relative to Nor), and 864 gained and 715 lost promoter elements and 2733 gained and 2466 lost enhancers in LNC (relative to EC) (Fig. 1b and Supplementary Fig. 5a). These numbers were comparable to those in previous studies in other cancers8,26. The altered elements exhibited significant gain or loss in analyzed patients (Supplementary Fig. 5b and Supplementary Data 3) or in a representative patient (Supplementary Fig. 5c), as indicated in heatmaps.

Next, these differential enhancers were classified into six groups: G1 (EC-specific gain; n = 985), G2 (LNC-specific gain; n = 837), G3 (common gain; n = 1974), G4 (EC-specific loss; n = 871), G5 (LNC-specific loss; n = 548), and G6 (common loss; n = 2715) (Fig. 1c and Supplementary Data 4). These results demonstrate that primary and metastatic cancers have a majority (59%) of common and a minority (41%) of location-specific epigenomic features, a pattern distinct from that of promoter variances, which showed only a small number of commonly altered promoter elements (Supplementary Fig. 5d, e). Unsupervised hierarchical clustering analysis was performed with these differential enhancers, and these three types of enhancer profiles were clearly separated into group-specific traits (Fig. 1d). Similarly, the differential promoter element- and enhancer-associated genes were subjected to principal component analysis (PCA), which successfully separated normal and cancer samples along PC1; in addition, EC and LNC samples were discriminated along PC2 (Fig. 1e). Among the top 100 altered element-associated genes (50 for gain and 50 for loss), genes with positive PC1 scores (such as GFI1, MSI2, RUNX3, CD72, and EXO1) in the promoter analysis and genes with negative PC1 scores (such as CD72, RUNX3, TNFSF8, NCF1C, and CTNNB1) in the enhancer analysis were considered potential oncogenes promoting carcinogenesis (Supplementary Fig. 6a, b). In addition, genes with positive and negative PC2 scores in the promoter analysis partially overlapped with top genes with positive and negative PC2 scores in the enhancer analysis respectively, including functional regulators involved in carcinogenesis and metastasis (Supplementary Fig. 6a, b).

Considering that the epigenetic discrepancy among the three groups in distal regions, the epigenetic features of which were more dynamic among multiple patients, was much more stable than that in proximal promoter regions (Supplementary Fig. 5), we focused on investigating enhancer reprogramming in subsequent analyses. For altered enhancer-linked genes, genes with positive PC1 scores were associated mainly with the regulation of RNA polymerase II activities and epidermal development involving the Rap1 and Ras signaling pathways, while genes with negative PC1 scores were related to cell migration, the immune response, and the T-cell receptor signaling (Supplementary Fig. 6c, d). In the positive PC2 dimension with enrichment in LNC samples, the enhancer-linked genes were involved mainly in immune cell functions as well as cell shape regulation, while genes enriched in the negative PC2 direction were associated with cell adhesion, cellular organization, and pathway regulation in cancer (Supplementary Fig. 6e, f). These data validate the accuracy of sample collection and demonstrate that genome-wide alterations in distal enhancer-linked genes are potential regulators of esophageal carcinogenesis and metastasis.

SEs, a subtype of enhancers with extended physical proximity, modulate the expression of master regulators of cell identity and disease states9,32. Accumulating evidence suggest that SEs are enriched in cancer signatures as a hallmark in multiple cancers23,25,26. To examine the role of SEs in ESCC, we identified 1042 nonredundant predicted SEs in total (Supplementary Data 5) using ROSE33 and found that the three cohorts shared 35.7% of SE regions, higher than the percentage of common typical enhancers (Fig. 1f, g; Supplementary Fig. 6g), indicating the evolutionary conservation of SE regions across patients and cell types. It is worth noticing that the vast majority of SEs were altered, of which 436 and 427 SEs were gained and 481 and 310 SEs were lost in EC and LNC, respectively (Supplementary Fig. 6g), validating the heterogeneity of predicted SEs23. When PCA was performed on all of these SEs, the three groups of cancer profiles were roughly separated, with stringent overlap (Supplementary Fig. 6h).

Relative to the normal esophageal epithelium, 392 (46.3%) and 494 (52.1%) cancer-specific SEs were recovered in the EC and LNC H3K27ac profiles, respectively, and these SEs were enriched with cancer-associated genes or oncogenes, such as CCND125,26, CWH43, ANO1, and CXCR4 (Fig. 1f, g and Supplementary Fig. 7). For example, neighboring regions of CCND1, CLUAP1, and MTMR8 showed a significant gain of promoter elements, enhancers, and SEs in cancer tissues, and loss of promoter elements, enhancers, and SEs was observed in TFAP2B, MAU2, and EHF, accompanied by transcriptional activation (CCND1) or repression (EHF) (Fig. 1h and Supplementary Figs. 6i, 7). Our consistent retrieval of previously identified distal enhancers (such as CCND1 at the rs7105934 locus) from other cancers25,26,34 in our unbiased profiling of H3K27ac enrichment further supports the reliability of our data.

Enhancer alterations correlate with aberrant transcriptional programs in ESCC

To explore relationships between epigenomic enhancer alterations and gene expression, we analyzed the RNA-seq data from the same cohorts as well as 22 transcriptomes from eight additional ESCC patients. As previously reported9, we assigned predicted enhancers to nearby genes with the nearest TSSs, and then analyzed the expression of differential enhancer-linked genes in the six groups (Fig. 1c). As expected, the expression of genes with nearby differential enhancers was generally positively correlated with H3K27ac enrichment signals, especially for the common gained (G3) and lost (G6) enhancers (Fig. 2a). In the groups with more altered enhancers, the numbers of upregulated genes in G3 and downregulated genes in G6 were much larger than those in the other groups (Fig. 2b). Genes linked to gained enhancers (G1–3) were generally expressed at higher levels in EC or LNC samples, while genes linked to lost enhancers (G4–6) were expressed at lower levels in cancer samples than in normal epithelium (Fig. 2c). For differential enhancer-linked genes with detectable expression levels, the expression of a relatively high percentage of DEGs followed the pattern of H3K27ac alterations, especially in the G3 (n = 413; 78.8%) and G6 (n = 442; 70.8%) groups, although the expression of some genes was not consistent with the patterns of their nearby H3K27ac marks (Supplementary Fig. 8a, b and Supplementary Data 6). Moreover, the total number of DEGs was larger than the number of altered enhancer-linked DEGs (Supplementary Fig. 8c and Supplementary Data 7), possibly because DESeq2 was used for DEG analysis (n = 18) without considering the pairing of patients and the threshold of ≥6/10 patients used for differential enhancer identification (Supplementary Fig. 3). These results indicate that enhancer acquisition and loss is a reliable and stable signature of primary and metastatic ESCC and is even more distinctive and predictive than aberrant gene expression (Supplementary Fig. 8c)35.

When considering the difference among the three cohorts, we showed that the expression of genes with gained enhancers was significantly upregulated and that of genes associated with lost enhancers was markedly downregulated, compared to that of control genes with unaltered enhancers, in the three groups of comparisons (Fig. 2d), further verifying the control of gene dysregulation by enhancer alterations. In addition, although a single altered enhancer was sufficient to induce aberrant gene expression, an increased number of gained enhancers was positively correlated with gene activation, while an increased number of lost enhancers was negatively correlated with gene repression in a quantitative manner (Fig. 2e). It further demonstrates that enhancer alterations and reprogramming contribute to the aberrant transcriptional program in primary and metastatic tumors. To identify biomarkers for EC and LNC, we selected G1-, G2-, G3-, and G6-associated genes as input, and analyzed the correlation between their mean expression levels and H3K27ac signals (exemplified by EHD3, ANO1, and EN1); for example, the H3K27ac level in an EHD3 enhancer was highly correlated with its expression (r = 0.86, Spearman correlation) (Supplementary Fig. 9a, band Supplementary Data 8). The top 100 genes with a positive correlation were ranked by their expression FCs, and the top 10 or 15 genes were assigned to subgroup g1 (FAM19A5, LAMC2, SNAI2, PMEPA1, etc.), g2 (TCF7, KLHL6, SCIMP, CXCR4, etc.), g3 (FCGR2A, HOXD11, ZNF469, CLDN3, etc.), and g6 (such as TFAP2B, MAL, EMP1, EHD3, etc.) with region-specific enhancer and transcript enrichment (Supplementary Fig. 9c, d and Supplementary Data 9). The combinations of these four groups of genes were defined as biomarkers for Nor (g1− g2−g3− g6+), EC (g1+ g2− g3+ g6−), and LNC (g1− g2+ g3+ g6−) tissues (Supplementary Fig. 9e).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses showed that gained enhancers in EC and LNC were enriched in essential cancer processes, such as cell proliferation, motility, and migration, and were associated with genes involved in proteoglycans, chemokine signaling, transcriptional misregulation, Rap1 signaling, and other pathways in cancer (Supplementary Fig. 10a, b). Importantly, consistent with the PCA results (Supplementary Fig. 6), gained enhancers in LNC were associated with cell migration as well as immune system processes regarding the T-cell receptor signaling, which might be elicited by the mixture of minimal immune cells from lymph nodes (Supplementary Fig. 10a, b). However, we also observed that relative to those in EC samples, the gained enhancers in LNC samples were enriched in cell adhesion, cell migration, and other metastatic traits, with the activation of cancer-related pathways (Supplementary Fig. 10c, d). Interestingly, some cancer-specific features, including PI3K-AKT signaling, HIF1 signaling, and Hippo signaling36, emerged as EC-specific features but were diminished in LNC samples (Supplementary Fig. 10d). Intriguingly, we found that some genes with nearby active enhancers were not transcriptionally activated in specific states. It has been proposed that the epigenetic patterning and information of enhancers may be established and embedded within enhancer elements before cellular processing37. Therefore, we selected the genes for which the nearby enhancers were pre-activated (such as SORL1, BCL2, KLHL6, and CLDN3) or pre-silenced (such as MAPK6, ETS2, LCN2, and LGALSL) before transcriptional alterations in EC samples, and considered that these enhancers or their associated genes might be the priming factors for the initiation of cancer metastasis-related processes (Supplementary Fig. 11 and Supplementary Data 10). Collectively, we propose that enhancer reprogramming confers aberrant gene expression in esophageal carcinogenesis and metastasis.

Super-enhancer signatures and heterogeneity in primary and lymph node tumors

A much higher proportion of SEs than typical enhancers were recurrently gained or lost in esophageal or LNC (Supplementary Figs. 5a and 6g). However, the active typical enhancers were relatively consistent in the same cohort of profiles among multiple patients and were more stable than in multiple cancer cell lines23, while SEs varied with a modest recurrent tendency (estimated recurrence rates of 6–16% in ≥5/10 or 5/8 patients) among multiple patients (Fig. 3a), indicating the heterogeneity of SEs in primary and lymph node tumors.

a Percentage of predicted typical enhancers and super-enhancers (SE) across Nor (n = 10), EC (n = 10), and LNC (n = 8) tissues showing H3K27ac enrichment above that of randomly selected regions (99%) across an increasing number of patients. b Gene Ontology (GO) analysis of the top 2000 predicted typical enhancer-linked genes and SE-associated genes; the top significantly associated biological processes are presented. c–e A total of 1317 super-enhancers are ranked by their differential H3K27ac intensity between EC and Nor, LNC and Nor, as well as LNC and EC samples. Genes associated with the top gained and lost SEs are listed; oncogenes are highlighted in red, and tumor suppressors are highlighted in blue. f, g Fold changes in H3K27ac enrichment signals (f) and expression levels for differential SE-associated genes (g) (EC/Nor, LNC/Nor, or LNC/EC) and for differential SE (gained and lost)-associated genes (EC vs. Nor, LNC vs. Nor, or LNC vs. EC). Unaltered SE-associated genes were used as controls. Data points represent means of Nor (n = 10), EC (n = 10), and LNC (n = 8) samples. Boxes correspond to interquartile ranges (IQR), thick black lines indicate the median values, and the whiskers extend to the lowest or highest data point that are still within 1.5 IQR of the bottom or top quartile, respectively. P value was calculated by the two-sided Wilcoxon test (no adjustments). *P < 10e-10, **P < 10e-20, ***P < 10e-30. h Cancer hallmark analysis using differentially predicted SEs showing recurrently gained, recurrently lost, and unaltered H3K27ac signals in EC or LNC relative to Nor. The log (adjust P value) obtained from the hypergeometric test is shown. i, j Survival analysis comparing groups of patients with high or low expression of the top common lost (i) or gained (j) SE-associated genes in ESCC TCGA data using Kaplan–Meier plotter. A poor prognosis observed for ESCC patients with tumors possessing a high expression signature of gained SE-associated genes and a low expression signature of lost SE-associated genes. Survival data are presented every 20 months. A Log-rank test was performed for survival data. Source data are provided as a Source Data file.

Full size image

GO analysis showed that the SEs in the EC and LNC cohorts were enriched in genes involved in the regulation of metabolic processes, signal transduction, responses to cytokines, cell proliferation, cell migration, and cell mobility, and the enrichment was more significant than that of the GO terms predicted from the top 2000 typical enhancers with the strongest average H3K27ac signals (Fig. 3b). Because the vast majority of SEs were recurrently gained or lost in EC or LNC samples, we propose that consistent with our previous study11, enhancer reprogramming (especially of SEs) may predominantly mediate cellular responses to cytokines, such as TGFβ signaling, to promote tumourigenesis and metastasis.

Using the previously reported strategy23,26 to map gained or lost SEs (Supplementary Fig. 6g), we identified some putative targets of the top altered enhancers (Supplementary Data 11), including some notably gained SE-linked oncogenes associated with cell activation, communication, adhesion, and motility (such as CCND1, GNAS, FGF3, HOXA10, and CEBPB), as well as many lost SE-linked genes (such as EHF, TFAP2B, NFIX, MAL, and MEIS1) involved in epithelial development, morphogenesis, and developmental processes (Fig. 3c, d and Supplementary Fig. 12a). Interestingly, the putative targets of lost SEs highly overlapped with those reported in RCC, and these genes were expressed at low levels in tumor tissues but at high levels in normal tissues (Fig. 3c, d)26, indicating the conserved mechanisms underlying the loss of epithelial characteristics. Moreover, the SEs gained in LNC relative to EC were enriched mainly in the immune system and oncogenic features (Supplementary Fig. 12b), as represented by FOXA1, CXCR4, SOX2, CCND2, etc. (Fig. 3e). Compared with the trends in transcriptional changes associated with enhancer alterations (Fig. 2d), the gained or lost SEs-associated genes were deregulated with increased FCs in gene upregulation or downregulation, respectively (Fig. 3f, g). We next performed cancer hallmark analysis9 to further explore the biological significance of the predicted SE heterogeneity with a gain or loss status. In line with the GO analysis results (Supplementary Fig. 12a, b), somatic gained SEs were enriched in genes related to angiogenesis, invasion, and evasion of growth suppressors, and this enrichment was more significant than that for the somatic loss and unaltered categories (Fig. 3h). These data indicate that gained somatic SEs in EC and LNC can be used to predict the progression and aggressiveness of ESCC and LNC, and highlight the importance of SEs in carcinogenesis.

The heterogeneity and variance of SEs in clinical sample profiling (Supplementary Fig. 12c) may lead to some top-ranked differential SEs exhibiting a large absolute difference value but a small FC. Therefore, we established another pipeline to rank the differential predicted SEs considering sample pairing for each patient with a threshold of ≥5/10 or 5/8 (for LNC analysis) patients (the SEs identified in individual samples and recurrently identified in ≥3 patients are listed in Supplementary Data 12 and Supplementary Fig. 12d). The numbers of gained and lost SEs obtained in each group were slightly greater than those calculated from the mean signals, with a large proportion of overlapping altered SEs (Supplementary Fig. 12e). To avoid missing essential SEs, we combined the differential SEs identified through the two strategies and ranked these SEs by their FC (Supplementary Fig. 12f and Supplementary Data 12). The SEs with the highest FC (FC >2 for gain and FC <0.7 for loss; because there were much fewer lost SEs than gained SEs) overlapped with those with the highest absolute difference values (>2 for gain and < −2 for loss), and dozens of high-confidence altered SEs were identified (Supplementary Fig. 12f). As shown in Supplementary Figs. 13a and 14(Supplementary Data 13), H3K27ac signal gain or loss was positively correlated with the upregulation or downregulation patterns of nearby genes, respectively, in a stringent manner, and numerous top representative genes with gained (such as FGF3, ANO1, HOXA9, and CCND1) or lost (such as NFIX, TFAP2B, and HOXA3) SEs overlapped with the results shown in Fig. 3c–e. These gained SE-linked genes in EC or LNC were also associated with tumor invasion and inflammation processes (Supplementary Fig. 13b). Taking together, we identified four groups of deregulated genes associated with the differential SEs (n = 456) (Supplementary Data 14), the majority of which were commonly gained or lost in EC and LNC, along with a minor group of genes with LNC-specific gained or lost SEs (exemplified by MEIS1 and CXCR4) (Supplementary Figs. 13c and 14), consistent with altered typical enhancer-linked genes (Supplementary Fig. 8b).

To explore the potential roles of the putative targets of the predicted SEs in other cancers, we analyzed the expression of differential SE-linked genes (Supplementary Fig. 13c) in primary and metastatic pancreatic adenocarcinoma (PDA) and colorectal cancer (CA) datasets from The Cancer Genome Atlas (TCGA). Surprisingly, we identified 78 and 212 DEGs that showed either common gain/loss or metastasis-specific gain/loss patterns in PDA and CA, respectively, among which the common upregulated (gain) and downregulated (loss) genes accounted for the majority of the DEGs (Supplementary Fig. 15a, b and Supplementary Data 15). The examined DEGs from ESCC, PDA, and CA overlapped, demonstrating that many DEGs are coincident among these three cancers, especially those in the groups with common gain or loss between ESCC and CA (Supplementary Fig. 15c); 12 genes, including ANO1, CTTN, PYGB, SOX4, SOX9, PMEPA1, etc., were commonly upregulated in all three types of cancers (Supplementary Fig. 15d and Supplementary Data 16), an effect possibly mediated by SE activation. We also identified six genes (KHDRBS2, IQSEC1, ARHGEF1, ARID5A, IRX5, and TMC6) that were activated in metastatic cancers in two of the three cancer types (Supplementary Fig. 15d) and might be associated with metastatic traits. To elucidate the clinical relationship between altered SEs and ESCC patient survival, we identified differential SE-linked genes that displayed a high correlation with the expression of their target genes and performed survival analysis via Kaplan–Meier Plotter with TCGA data. Higher expression of common lost SE-associated genes predicted a higher probability of patient survival (Fig. 3i), while higher expression of common gained SE-associated genes predicted a lower probability of survival in EC, PDA, and CA patients (Fig. 3j and Supplementary Fig. 15e–g). It demonstrates the close relationship between cancer-associated SEs and clinical events during ESCC carcinogenesis.

Transcription factor circuitries of ESCC and LNC

The core circuitry of cell type- or status-specific transcriptional regulation is often dominated by several key TFs to establish context dependence38. To identify principal TFs that promote ESCC tumourigenesis and metastasis, we predicted trans-acting factors associated with common gained or LNC-specific gained enhancers using HOMER39. Binding motifs for many essential TFs, such as RXRA, NFE2L2, ESRRA, IRF2, RELA, ESRRA, SOX2, and SMAD2/3 (Fig. 4a, b), key TFs that mediate TGFβ signaling in epithelial-mesenchymal transition (EMT) and cancer metastasis11,40, were enriched in common gained or LNC-specific gained enhancers, with some overlapping TFs. Core regulatory circuitry analysis25,41further showed that many TFs predicted by gained enhancers, including SMAD3, RXRA, IRF2, ETS1, etc., were displayed among the top rank-ordered TFs interacting with other predicted TFs (Supplementary Data 17); the identification of FOSL2, SOX2, and RXRA (Fig. 4c), which were also previously revealed in ependymoma25, supports the reliability of our data and strategy. Notably, the majority of the top enriched TFs were highly expressed in the EC and/or LNC groups compared to the normal esophageal epithelium group (Fig. 4d and Supplementary Data 18). The trans-acting factors predicted by commonly lost enhancers, such as TBX21, NFIA, MEIS1, and ELF3, were involved in epithelial development (Supplementary Fig. 16a), while the functions of some LNC-specific lost enhancer predicted TFs remained largely unclear in LNC (Supplementary Fig. 16b). Interestingly, many TFs predicted by commonly lost enhancers were negatively correlated with TFs predicted by common gained or LNC-specific gained enhancers (Supplementary Fig. 16c and Supplementary Data 19), indicating the distinctive roles of these TFs in normal and cancer tissues.

Two independent short hairpin RNAs (shRNAs) for each target were used to investigate the functional roles of altered enhancer-predicted top-ranked TFs by efficiently knocking down the corresponding TFs (RXRA, NFE2L2, ZNF519, ESRRA, RELA, TCF4, SOX2, ETS1, SOX12, NFATC2, IRF2, and ZEB1), with individual low efficient shRNAs (i.e., ETS1 shRNA1 and ZNF519 shRNA2) (Fig. 4e). Then, cell growth was assessed in TE1 esophageal cancer cells expressing the control shRNA or shRNAs targeting these 12 TFs, showing that the depletion of IRF2, RXRA, RELA, SOX12, and NFATC2 resulted in a significant decrease in cell proliferation and that the knockdown of ETS1 (with shRNA2), ZEB1, TCF3, and NFE2L2 led to a modest decrease in cell growth (Fig. 4f and Supplementary Fig. 17a). Wound healing assays further demonstrated that depletion of NFE2L2, ESRRA, TCF4, SOX2, ETS1 (with shRNA2), SOX12 (weak but significant), and IRF2 resulted in impaired cell migration (Fig. 4g and Supplementary Fig. 17b). Consistent with these cellular phenotypes, the knockdown of some key TFs, such as IRF2, RELA, RXRA, and NFE2L2, induced upregulation of the cell cycle inhibitor p21 and downregulation of the oncogene C-Myc and the EMT inducers SNAI1 and FN1 (Supplementary Fig. 17c).

Moreover, the binding sites of eight TFs in cancer cell lines (data available from ENCODE; the data source is provided in Supplementary Data 20) largely overlapped with predicted enhancers (Supplementary Fig. 18a); dozens to hundreds of common or LNC-specific gained enhancers that were associated with important oncogenesis- or metastasis-related genes, such as ZXDC, CTNNB1, CXCR5, CCND2, STAT3, and SNAI1, were bound by IRF2, RXRA, RELA, or ZEB1 (Supplementary Fig. 18b and Supplementary Data 20). Thus, we propose that these TFs mediate the enhancer gain process. In summary, we successfully identified several important TFs that regulate ESCC carcinogenesis and metastasis through the prediction of trans-TFs by gained enhancers.

Gained SE maps identify candidate drugs against ESCC carcinogenesis and metastasis

To screen candidate drugs for clinical use in ESCC and lymph node metastasis, we performed an integrated analysis of our tumor-specific SE-linked genes with the Washington University Drug Gene interaction database42, informing 69 common gained SE-linked genes and 43 LNC-specific gained SE-linked genes, which were redundantly categorized into 21 and 18 classes, respectively, including druggable genome, kinase, clinically actionable, cell surface, TF binding, etc. (Fig. 5a, b and Supplementary Data 21). Among the 49 and 25 putative druggable genome targets, our analysis revealed 236 drugs interacting with 36 druggable genes associated with commonly gained SEs and 170 drugs interacting with 18 druggable genes associated with LNC-specific gained SEs (Supplementary Fig. 19a and Supplementary Data 21). Because many genes can be targeted by several drugs, we ranked these genes according to the numbers of interactions with drugs, and the top eight druggable targets with the most interactions (HSP90AA1, CCND1, ANO1, and CCR3 from the commonly gained SEs; and BCL2, PDE4B, ROCK1, and CXCR4 from the LNC-specific gained SEs), which were identified as candidate genes responsive to small molecule inhibitors, were selected for subsequent functional investigations using commercial chemical inhibitors (Supplementary Fig. 19b). Except for the ROCK1 and CXCR4 inhibitors, six inhibitors, particularly the HSP90AA1 and PDE4B inhibitors, markedly suppressed TE1 cell growth in a dose-dependent manner, and these two inhibitors exhibited synergistic effects on cell proliferation (Fig. 5c and Supplementary Fig. 19c, d). Moreover, a high rate of cell apoptosis was induced by HSP90AA1 or PDE4B inhibitors or dual inhibitors (Supplementary Fig. 19e, f). Notably, the inhibitory effects on cell proliferation (Supplementary Fig. 20a–d) and proapoptotic effects (Supplementary Fig. 20e, f) of HSP90AA1 inhibitor and PDE4B inhibitor were perfectly validated in another two ESCC lines (KYSE30 and KYSE150). Thus, we demonstrate that the inhibition of HSP90AA1 or PDE4B activity substantially suppresses cell proliferation and induces cell apoptosis.

To assess the functions of these gained SE-linked putative targets in cancer metastasis, we performed wound healing assays upon inhibitors treatment and found that nearly all of these applied inhibitors showed strong (HSP90AA1, PDE4B, and BCL2) or weak (ANO1, CCR3, CCND1, and CXCR4) suppressive effects on cell migration (Fig. 5d and Supplementary Fig. 19g). Furthermore, treatment of immunodeficient mice bearing esophageal cancer xenografts (TE1 cell-derived) with alvespimycin hydrochloride (AH, an HSP90AA1 inhibitor) or ML-030 (a PDE4B inhibitor) remarkably suppressed tumor growth (Fig. 5e–g). Coincidently, the expression of HSP90AA1 and PDE4B as well as their nearby SEs was more significantly enriched in the EC and LNC groups than in the normal control group (Fig. 5hand Supplementary Fig. 21a, b), validating their potential as clinically achievable antitumor targets. Considering the high expression of HSP90AA1 relative to PDE4B in EC and LNC samples, we prioritized HSP90AA1 as a potential clinical diagnostic and therapeutic target. We divided 18 EC patients recruited in the present study into two groups with high or low HSP90AA1 expression in EC samples (Supplementary Fig. 21c) and found that the 2-year survival rate of HSP90AA1-low patients (66.7%, n = 6/9) was twofold higher than that of HSP90AA1-high patients (33.3%, n = 3/9) (Fig. 5i). Although the survival curve showed no significant difference, the HSP90AA1-low group displayed a better survival trend in our analysis of a limited number of patients (Supplementary Fig. 21d). To support our hypothesis, we analyzed the survival probability of multiple types of cancer patients in the TCGA database, showing that lower HSP90AA1 expression predicted better survival and prognosis in ESCC, EAC, head-neck squamous cell carcinoma, and liver cancer (Fig. 5j and Supplementary Fig. 21e–g). Considering the specific gain of PDE4B SEs, we tested the effect of PDE4B inhibition in cancer metastasis. As we expected, ML-030 treatment remarkably reduced the metastasis of KYSE150 cell-derived primary tumors in an intramuscular injection model to mimic primary EC (Supplementary Fig. 21h). These data suggest that the SE landscape can inform therapeutic targets, such as HSP90AA1 and PDE4B, to combat ESCC carcinogenesis and metastasis. Given that these gained SEs might be mediated by key TFs (Fig. 3), we analyzed the chromatin architectures around representative gained SEs near the ANO1, HSP90AA1, and PMEPA1 genes, showing that enhancer gain or loss occurred mainly within the same topologically associating domains (TADs), which were analyzed from a set of Hi-C data in fetal lung fibroblasts (IMR90)11, with some predicted TFs binding sites (Supplementary Fig. 22a). Therefore, we examined the expression of three putative targets of gained SEs (ANO1, HSP90AA1, and PMEPA1) upon TF depletion, demonstrating that the knockdown of some key TFs, which may mediate SE gain, resulted in the downregulation of ANO1, HSP90AA1, or PMEPA1 expression (Supplementary Fig. 22b). Collectively, our results identify primary and metastatic tumor-specific enhancers and SEs as well as their driven TFs, reveal a molecular basis for the regulation of oncogenic transcriptional programs, and provide a strategy for discovering therapeutic targets involved in ESCC carcinogenesis and metastasis.

Esophageal cancer is a severe disease with high mortality, accompanied by poor prognosis and frequent metastasis to lymph nodes; in particular, its survival rate is inversely correlated with its metastasis and recurrence. In addition to surgery and chemotherapy, very few chemical drugs targeting biological molecules are clinically approved1,2,3,4. Therefore, for this kind of clinically heterogeneous disease, understanding its epigenomic deregulation signatures has emerged as a pivotal strategy to identify biomarkers, vulnerabilities, and therapeutic targets for this clinically heterogeneous disease23. By profiling histone H3K27ac modifications in paired normal, primary, and lymph node tumor pairs, we generated a comprehensive compendium of promoter and enhancer alterations in ESCC.

Our study characterizes a comprehensive epigenetic profile of primary and metastatic ESCC tumors beyond genomic alterations17 and provides a valuable resource of epigenomic information for investigating cancer metastasis. Similar to a recent report regarding DNA methylation alterations in ESCC43, the significance of our work was partially limited by the limited sample size; a larger number of patients may help to dissect more clinical and therapeutic highlights, such as clinical subgrouping and personalized treatment, as reported previously in EAC44. Second, we show that the epigenomic and transcriptomic profiles of LNC are quite similar to those of primary tumors, with some LNC-specific characteristics. It suggests that targeting the common gained enhancer-linked oncogenic targets in primary and metastatic tumors, such as HSP90AA1, is a potential strategy for ESCC therapy. In addition, the blockade of LNC-specific enhancer-linked genes may specifically inhibit cancer cell metastasis. It demonstrates that epigenomic landscapes, which might be more stable than transcriptional differences35, are able to dissect the molecular differences between similar tumor entities. Third, we demonstrate that variant enhancer and super-enhancer reprogramming contributes to transcriptional remodeling, which is further emphasized by the reverse engineering of core transcriptional regulatory circuitries functionally involved in tumor cell proliferation and migration (Fig. 4). Finally, integrated analysis of somatically gained SEs with drug interaction databases enables us to identify and validate potent master regulators and cancer dependencies that are responsive to pharmacologic drugs or inhibitors (Fig. 5), highlighting the reliability of the enhancer landscape to inform precision therapies obtained from drug databases.

In recent years, epigenomic studies regarding promoter alterations24,45 and distal enhancer dynamics have illustrated a fundamental role for epigenomic elements in controlling the cell status in human diseases and cancers8,9,23,25,26,46, which motivate us to investigate a consistent epigenomic feature of carcinogenesis and metastasis in esophageal cancer. Because of the limitations of obtaining an appropriate number of cells with high purity from donors for ChIP-seq analysis, here we mapped only the H3K27ac profiles (Fig. 1), similar to the approach used in a recent report regarding ependymoma;25 however, although enhancer-specific histone H3 lysine 4 mono-methylation (H3K4me1) profiling will be quite helpful for precisely separating active promoter- and enhancer-enriched H3K27ac signals in promoter regions26. Therefore, we focused mainly on the distal enhancers outside promoter regions, which account for the majority of enhancer-linked cis-regulatory elements47. The identification of promoters herein using H3K27ac alone was not accurate, so we defined TSS ±2 kb regions with significant H3K27ac occupancy as active promoter elements (including promoters and a few proximal enhancers) but not as active promoters. Considering the neighboring occupancy between enhancer-marked H3K27ac and H3K4me3 (tri-methylation) around TSSs, the number of our identified active promoter elements was comparable to that of H3K27ac + /H3K4me3 + -defined active promoters within TSS ±2 kb regions reported in another work26, indicating similarities between active promoter elements and active promoters.

Although a large number of nearest deregulated genes were associated with altered enhancer elements (Fig. 2), there were still hundreds of DEGs that were not positively correlated with nearest H3K27ac alterations (Supplementary Fig. 8). As reported previously, approximately half of the predicted enhancer/gene or SE/gene interactions directly regulate the nearest proximal genes23,48. To build a perfect distal interaction model capable of identifying enhancers and paired target genes within a TAD, it will be meaningful to generate ESCC-specific Hi-C data; alternatively, integrating the altered enhancers or SEs with data from interaction datasets, such as PreSTIGE30, GREAT31, and RNAPII ChIA-PET, is another feasible strategy23 to accurately map the functional distal enhancers. Additionally, we identified a subgroup of genes, with delayed expression changes after H3K27ac alterations (Supplementary Fig. 11), which indicates that these enhancers may prime the transcriptional regulation of their targets. This possibility may also account for part of the inconsistency between gene expression and nearby H3K27ac alterations. Intriguingly, it has been proposed that the adjacent “normal” tissues may contain accumulated prevalent genomic alterations49,50,51,52. Investigating the link between epigenomic alterations and genomic mutations in the adjacent normal tissues will be quite meaningful.

Compared with those in matched nonmalignant esophageal tissues, recurrent predicted SEs in EC and LNC tissues were identified at known oncogenes, such as CCND1, CEBPB, and RXRA25,26, and at essential oncogenic promoters. We also observed a wide range of SE variations between different patients, with a heterogeneity markedly exceeding that of typical enhancers (Fig. 3). Although it has been proposed that over 60% of SEs are tissue- or cancer type-specific23, the top differential SE-associated genes were quite similar to those in other types of cancer20,23,25,26, highlighting the importance of these top gained SEs in carcinogenesis (Fig. 3b). We attempted to combine the absolute difference and FC data to successfully examine the highly reliable most altered SEs (Supplementary Figs. 12, 13), which might be a strategy worth recommending for other studies. The high overlap of gained and lost enhancers/SEs between ESCC and LNC (Figs. 1 and 3) validate the origin of metastatic LNC cells from primary ESCC and suggests the relative consistency of epigenomic features between primary and metastatic tumors, as well as the molecular consistency with other types of cancer (Supplementary Fig. 15).

Meanwhile, we recovered a series of potential tumor suppressors (EHF, MAL, and TFAP2B)26that were associated with lost enhancers and could be uncovered only in normal tissues (Fig. 3c–e). We also identified a series of LNC-specific enhancers and SEs; however, we could not exclude the minor possibility of specific gain of a few enhancers reflected by a minimal number of immune cells. Actually, many gained enhancers in LNC tissues relative to EC, which were associated with oncogenes, were also gained in EC relative to normal tissues (Supplementary Fig. 13). Common gained/lost enhancers might be more reliable than LNC-specific altered enhancers. In support of this notion, (Fig. 3i, j), our survival analysis indicated that somatic gained SEs contribute to aberrant gene expression and predict a poor clinical prognosis.

Our epigenetic landscape contains an appealing pool of potential targets that can contribute to ESCC tumorigenesis and metastasis. We performed trans-element prediction using gained enhancers (Fig. 4) and integrated SEs with drug interaction databases to predict potential targets responsive to inhibitors (Fig. 5). By elucidating core transcriptional regulatory circuitries from the set of gained enhancers, we identified several key TFs, including RXRA, NFE2L2, IRF2, RELA, etc., involved in cell growth and migration (Fig. 4f, g), some of which were also revealed to have tumor-promoting functions in ependymoma and RCC25,26. Clarification of the potential roles of these TFs in establishing cancer-specific enhancer profiles will be of great importance. Consistent with this presumption, among TF circuitries predicted from gained enhancers, SOX2 has been previously identified as a factor of the core regulatory circuitry for controlling epigenetic and transcription patterns in ESCC cell lines53.

Our notable finding is the identification of HSP90AA1 and PDE4B as potent therapeutic targets, both of which have been previously proposed as potential biomarkers for ESCC54,55,56, further supporting our findings here. Importantly, the pharmacologic inhibition of these two molecules significantly suppressed cell proliferation, migration, and xenograft tumor growth (Fig. 5). PDE4B is a protective cyclic AMP-phosphodiesterase involved in heart failure and colon cancer57,58. Although the clinical relevance between PDE4B and ESCC is not as strong as that between HSP90AA1 and ESCC, PDE4B inhibitor exhibits very strong antitumor effects. Therefore, the mechanisms underlying the tumorigenic role of PDE4B remain to be elucidated in the future. HSP90AA1, a heat shock protein critical for the stability of its target proteins, is important for autophagy and drug resistance59,60. Here we validated its positive correlation with tumor progression and cancer mortality, and emphasized its role in ESCC carcinogenesis and metastasis.

In summary, we demonstrate that enhancer and SE deregulation attributes primarily to transcriptional reprogramming during carcinogenesis and metastasis and that LNC enhancer signatures are quite similar to that of primary ESCC but possess metastasis-specific features. By trans-element analysis from the set of altered enhancers, we successfully identified IRF2, RELA, NFATC2, etc., as essential TFs involved in cancer cell growth and migration, possibly via establishing oncogenic enhancer programs. SE-linked HSP90AA1 and PDE4B, emerge as crucial oncogenes and potential therapeutic targets for esophageal carcinogenesis.