CDKN2A LoF drives BE and EAC evolution, but not EAC initiation

We collected whole-genome sequencing (WGS), whole-exome sequencing (WES) and gene panel sequencing data for 1,032 EACs from the literature6,32,33,34,35,36,37,38 or sequenced de novo by the Esophageal Cancer Clinical and Molecular Stratification (OCCAMS) Consortium (Supplementary Table 1). Our cohort reflected EAC high male prevalence, with almost 9:1 male-to-female incidence ratio39 (Supplementary Table 1). To ensure consistency, we annotated damaging mutations and copy-number alterations in all datasets using the same approach (Methods and Extended Data Fig. 1a–e). Because CDKN2A can be silenced also via epigenetic modifications, we analyzed methylation data for a subset of EACs32,40 (Supplementary Table 1). We then identified the damaged drivers in each sample using a curated list of 54 known (canonical) EAC drivers (Supplementary Table 2). In agreement with previous studies29,32,41, CDKN2A was the second most frequently damaged EAC driver, with LoF in 25% of samples (Fig. 1b). More than 56% of EACs (90% considering also TP53) had damaging alterations in other cell cycle regulators (Fig. 1c and Supplementary Table 2), suggesting that cell cycle disruption is key in EAC evolution but does not always involve CDKN2A.

Next, we measured the frequency of CDKN2A LoF in 257 BEs that progressed to high-grade dysplasia or EAC (P-BEs), again sequenced for this study or gathered from published datasets15,40,42,43,44 (Supplementary Table 1 and Extended Data Fig. 1f–j). CDKN2A LoF occurred significantly more frequently in P-BE than EAC (P = 4 × 10−9, two-sided Fisher’s exact test; Fig. 1d), suggesting that EAC does not always originate from a CDKN2A-damaged BE. To further investigate this, we analyzed 66 matched EAC-BE pairs with CDKN2A LoF in BE or EAC (Supplementary Table 1). Only 15 matched lesions had either identical or clonally related CDKN2A alterations (Fig. 1e), confirming that CDKN2A LoF is not required for precancer to cancer transition. Interestingly, 28 EACs lost CDKN2A independently of the paired BEs (Fig. 1e), suggesting that either EAC developed from a different CDKN2A-damaged BE clone or CDKN2A LoF was acquired after transformation.

Finally, we analyzed 99 BEs that did not progress to high-grade dysplasia or EAC (NP-BEs)15,40,43,44 (Supplementary Table 1 and Extended Data Fig. 1f–j). The frequency of CDKN2A LoF in NP-BE was even higher than P-BE and EAC (P = 3 × 10−3 and P = 3 × 10−13, respectively, two-sided Fisher’s exact test, Fig. 1f). Moreover, although in EAC, the dysregulation of cell cycle could occur through alterations of other genes, CDKN2A was the only gene encoding a cell cycle regulator damaged in BE (Fig. 1d). Therefore, unlike EAC, only CDKN2A LoF is relevant for BE evolution.

As observed previously22,45, P-BEs had significantly more damaged drivers than NP-BEs (P = 7 × 10−6, two-sided Fisher’s exact test; Supplementary Table 2), indicating that EAC initiation requires several driver events, most frequently TP53 complete loss. Given its high recurrence, we used TP53 LoF to assess the role of CDKN2A LoF in EAC initiation calculating the odds of cancer progression based on the mutational status of CDKN2A and TP53 in BE. As expected, the odds of cancer progression in BE cases with TP53 LoF was 1 irrespective of CDKN2A status (Supplementary Table 3), confirming that TP53 is a strong driver of EAC initiation. However, the odds of cancer progression in BEs with CDKN2A LoF and wild-type TP53 was lower than those of BEs with both wild-type genes (0.58 and 0.72, respectively; Supplementary Table 3). This suggested that an early occurrence of CDKN2A LoF in BE may reduce the likelihood of EAC initiation. To test this further, we compared two logistic regression models, one assuming a role in EAC initiation only for TP53 LoF (model 1) and the other for both TP53 and CDKN2A LoFs (model 2; Methods). Model 2 was a significantly better predictor of EAC initiation than model 1 (P = 0.01, ANOVA test), with expected occurrences of P-BEs with any status of TP53 and CDKN2A perfectly matching the observed occurrences (Supplementary Table 3). The negative β coefficient of CDKN2A in model 2 further confirmed that CDKN2A LoF may reduce risk of cancer progression (Methods and Supplementary Table 3).

TP53 loss reduces proliferation of CDKN2A LoF BE cells

Next, we set out to investigate how CDKN2A LoF in BE could prevent EAC initiation. As the proportion of BEs with both CDKN2A and TP53 LoF was significantly lower than that of BEs with CDKN2A LoF only (P = 0.05, two-sided Fisher’s exact test; Fig. 2a), we hypothesized that negative selection might act on BE cells losing both genes. To test this hypothesis, we compared CDKN2A and TP53 LoF clonality in 580 EACs with WGS or WES data, as clonality informs on when alterations are acquired during cancer evolution. Despite the well-known EAC intratumor heterogeneity14, CDKN2A or TP53 LoFs were clonal in almost 70% of EACs (397/580), confirming that both alterations are early events. However, EACs with fully clonal CDKN2A LoF were significantly fewer than those with fully clonal TP53 LoF (P = 0.001, two-sided Fisher’s exact test; Fig. 2b), suggesting that overallTP53 LoF tends to predate CDKN2A LoF. In support of this, CDKN2A LoF occurred before TP53 LoF in only 6% of the 47 EACs with LoF alterations in both genes as compared to 38% where TP53 LoF occurred before that of CDKN2A (Fig. 2c). This finding confirmed that the subsequent loss of TP53 in the presence of CDKN2A LoF is a rare event, suggesting that it might be selected against.

Fig. 2: Effect of TP53 loss in BE with CDKN2A LoF.

a, Frequency of CDKN2A LoF in 356 BEs (n = 257 P-BE and n = 99 NP-BE individuals, respectively) with or without TP53 LoF. Statistical significance was assessed using a two-sided Fisher’s exact test (P = 0.05). b, Frequency of EACs with clonal LoF alterations in CDKN2A and TP53 genes. For this analysis, n = 580/779 patients with EAC with WGS or WES data and LoF in these genes were considered. Statistical significance was assessed using a two-sided Fisher’s exact test (P = 0.001). c, Frequency of EACs with clonal and subclonal LoF alterations in CDKN2A and TP53 genes in n = 47 patients with WGS or WES data and damaging alterations in both genes. d, CDKN2A and TP53 gene expression levels quantified by RT–qPCR of RNA from TP53 wild-type CP-A cells (CP-A_TP53wt), three TP53 KO clones (CP-A_2c8, CP-A_3d2, CP-A_5f4) and control RNA relativized to ACTB expression. One biological replicate was performed with three technical replicates. e, TP53 gene structure in CP-A_TP53wt, CP-A_2c8, CP-A_3d2 and CP-A_5f4 cells. Exon-intron arrangement was derived from the UCSC genome browser (https://genome.ucsc.edu/) based on NM_000546 mRNA sequence (chr17:7,668,421-7,687,490, hg38 assembly). Dotted lines represent edited regions. f, Growth curves of CP-A_TP53wt, CP-A_2c8, CP-A_3d2 and CP-A_5f4 cells. Proliferation was assessed every 24 h and normalized to time zero. Mean values at 72 h were compared by two-tailed Student’s t-test (P = 1 × 10−4, 8 × 10−6 and 2 × 10−7, respectively). Error bars show standard deviation. Three biological replicates were performed, each in two to four technical replicates. ctrl, control; KO, knockout; P-BE, progressor Barrett’s esophagus; RLU, relative light unit; RT–qPCR, real-time quantitative PCR; RQ, relative quantification; UTR, untranslated region; wt, wild type.

Source data

Interestingly, BAR-T cells, derived from BE with constitutive loss of CDKN2A, increase cell doubling times upon TP53 knockdown46, supporting the hypothesis that the additional loss of TP53 reduces cell growth rate. To test this experimentally, we induced TP53 knockout (KO) in metaplastic BE CP-A cells derived from a male individual with CDKN2A LoF and wild-type TP53 (ref. 47). First, we confirmed that CP-A cells expressed TP53 but did not express CDKN2A (Fig. 2d). We then used CRISPR-Cas9 to edit TP53 (Supplementary Table 4) and performed single cell cloning to expand cell colonies. To control for off target effects and clonal differences, we selected three clones with a partial deletion of TP53 exons 5 and 6 (Fig. 2e), as assessed via amplicon sequencing (Supplementary Table 4). We confirmed that these clones did not express CDKN2A nor TP53 (Fig. 2d). The fact that we could isolate clones losing both genes implied that BE cells with CDKN2A LoF can survive subsequent TP53 loss. However, compared to TP53 wild-type CP-A cells, all three TP53 KO CP-A clones showed significantly slower growth rate that was already visible after 72 h (two-sided t-test test, Fig. 2f).

This finding was in line with the reported increase in cell doubling times of TP53 knockdown BAR-T cells46 and supported the tumor-preventive role of early CDKN2A inactivation due to the reduced fitness, defined as proliferative capacity, of cells additionally losing TP53.

LoF of 9p21 genes predicts poor survival in EAC, but not in BE

Because CDKN2A LoF has been associated with poor patient survival5,6,7, we investigated the survival effect of CDKN2A and other 9p21 gene LoF in our extended BE and EAC cohorts. Patients with EAC and CDKN2A LoF showed significantly worse survival than those with the wild-type gene (Fig. 3a). This difference held true even when patients with CDKN2A homozygous deletions (Fig. 3b) or damaging mutations (Fig. 3c) were considered separately. However, we did not observe lower survival in patients with CDKN2A heterozygous deletions only (Extended Data Fig. 2a), suggesting that CDKN2A complete loss is required to affect prognosis. Damaging alterations in TP53 or other cell cycle regulators had no effect on survival (Extended Data Fig. 2b–f) despite their frequent EAC alterations (Fig. 1c). Therefore, the survival effect of CDKN2A LoF does not depend on its function as cell cycle regulator. Moreover, CDKN2A LoF was not a predictor of worse survival in P-BE (Fig. 3d), again suggesting context-dependent consequences of its loss.

Fig. 3: Effect of the LoF of CDKN2A and other 9p21 genes on survival.

a–c, Kaplan–Meier survival curves of n = 1,032 patients with EAC with wild-type CDKN2A compared to those with all types of LoFs (P = 2 × 10−4) (a), only homozygous deletions (P = 6 × 10−3) (b) and only LoF mutations (P = 3 × 10−3) (c). d, Kaplan–Meier survival curves of n = 129 patients with P-BE with and without CDKN2A LoF. Log-rank method was used to estimate the P values. ns, not significant. e, Approach to test the effect of the co-damage in 9p21 genes on patient survival. Only n = 779 patients with EAC with WGS or WES data were used for the survival analysis, whereas n = 337 patients with RNA-seq data were used to measure 9p21 gene expression. Letters correspond to the 26 genes according to their order in the chromosomal locus. f, LoF frequency of 9p21 genes in n = 779 patients with EAC. g, Distribution of normalized expression values in the 9p21 genes in n = 337 patients with EAC. Boxplot shows first and third quartiles, whiskers extend to the lowest and highest value within the 1.5× interquartile range and the line indicates the median. h, Kaplan–Meier survival analysis of patients with EAC with co-alterations in the ten expressed 9p21 genes and n = 413 patients with EAC with a wild-type locus. Only groups with significantly poor survival (FDR < 0.1) are shown and genes of interest are outlined in black. All groups used in the analysis are listed in Supplementary Table 5. The minimum and maximum number and percent of damaged EACs in f and h are reported in the corresponding heatmap. HD, homozygous deletion; WES, whole-exome sequencing; WGS, whole-genome sequencing. Cartoon in (e) was created with BioRender.com.

Source data

We then investigated whether the co-occurring loss of other 9p21 genes could also contribute to poor survival, restricting the analysis to 779 EACs with WGS or WES data (Fig. 3e). Although CDKN2A was the most frequently occurring alteration in the locus, confirming that it is the event under positive selection, the other 25 genes were frequently co-lost with it (Fig. 3f). However, only ten 9p21 genes were expressed in EAC (Fig. 3g) or normal esophagus (Extended Data Fig. 3), suggesting that the loss of the remaining 16 genes likely had no functional consequences. We therefore tested the potential impact on survival of the ten 9p21 expressed genes by dividing patients with EAC in nine groups. Each of these groups represented at least 5% of the cohort and was composed of patients with the same 9p21 mutation and copy-number profile (Supplementary Table 5). Patients in all nine groups had worse survival than 413 patients with EAC with a wild-type 9p21 locus (FDR < 0.1; Fig. 3h and Supplementary Table 5). All patients lost KLHL9, IFNE, MTAP, CDKN2A, CDKN2B and DMRTA1 (Fig. 3h), suggesting that alterations in these genes may contribute to poor prognosis.

LoF of 9p21 genes has distinct consequences in BE and EAC

Our results suggested that the LoFs of CDKN2A and other 9p21 genes have functional and survival consequences that depend on time and context. Disentangling these variable effects is challenging because 9p21 genes are often co-damaged (Fig. 3f). To tease out the contribution of individual 9p21 genes, we divided 22 NP-BEs, 108 P-BEs and 337 EACs with matched genomic and transcriptomic data (Supplementary Table 1) into four groups (Fig. 4a). Each group had the same LoF profile of the six genes whose loss impacted survival (KLHL9, IFNE, MTAP, CDKN2A, CDKN2B and DMRTA1; Fig. 3h). Group 1 included all samples with CDKN2A LoF independently of the status of the other genes (Fig. 4b), closely resembling the cohorts tested in the survival analysis (Fig. 3a,d). The other three groups were subsets of group 1 with variable LoF frequency in the six genes (Fig. 4b).

Fig. 4: Functional consequences of 9p21 gene LoF in BE and EAC.

a, Frequency of damaged 9p21 genes in the four groups estimated over n = 22 patients with NP-BE, n = 108 patients with P-BE and n = 337 patients with EAC with matched genomic and transcriptomic data. b, Proportions of samples with LoF in KLHL9 (N), IFNE (T), MTAP (U), and CDKN2B (W) and DMRTA1 (X) over samples with CDKN2A LoF (V) in each group of NP-BEs, P-BEs and EACs. The number of samples in each group and condition is reported. c, Relative proportion of dysregulated pathways in NP-BE, P-BE and EAC cohorts mapping to cell cycle regulation, metabolism, signal transduction, immune response and development. Numbers in brackets represent the number of unique pathways. d–f, Results of pre-ranked GSEA48 showing the normalized enrichment score (NES), FDR and gene ratio (number of leading-edge genes over the total expressed genes) of pathways dysregulated in each group of NP-BEs (d), P-BEs (e) and EACs (f). NES > 0 indicates pathway upregulation, whereas NES < 0 indicates downregulation. P values were estimated by permutation and corrected for multiple testing using the Benjamini–Hochberg method. g,h, Fold change of expression and correlation plot of the shared leading-edge (LE) genes of interferon gamma (g) and alpha (h) response pathways enriched in P-BE and EAC group 2 as compared to 9p21 wild-type samples. Coefficients and associated P values from two-sided Pearson’s correlation test are reported for both pathways. i, Overlap of leading-edge genes between interferon gamma and alpha response pathways enriched in P-BE and EAC group 2. The 19 shared genes are listed. FC, fold change.

Source data

We identified the dysregulated biological processes in each group as compared to the corresponding 9p21 wild-type samples by performing a pre-ranked gene set enrichment analysis (GSEA)48 in NP-BEs, P-BEs and EAC separately. Overall, we detected 72, 62 and 28 unique pathways significantly dysregulated (FDR ≤ 0.01) in NP-BE, P-BE and EAC, respectively (Supplementary Table 6). Almost 80% of these pathways mapped to only five biological processes, namely cell cycle regulation, metabolism, immune response, signal transduction, and development. Overall NP-BE and P-BE showed a higher fraction of dysregulated pathways than EAC (Fig. 4c), suggesting that 9p21 LoF had higher impact in premalignant conditions.

As expected, given CDKN2A, CDKN2B and KLHL9 role in cell cycle regulation role, we found cell cycle dysregulation across groups and conditions except group 4 (CDKN2A LoF only; Fig. 4d–f and Supplementary Table 6), suggesting that the co-deletion of KLHL9, CDKN2A and CDKN2B maximizes the effect.

CDKN2A LoF alone might not be sufficient also to trigger metabolic or immune dysregulation (Fig. 4d–f and Supplementary Table 6). In this case MTAP and IFNE LoF could play a role given their functions in metabolic reprogramming49,50 and activation of immune response through metabolic regulation51, respectively. Interestingly, oxidative phosphorylation was consistently downregulated in NP-BE, upregulated in P-BE, and showed no difference in EAC (Fig. 4d–f and Supplementary Table 6). This once again suggested that the same genetic alterations may trigger different functional responses depending on the context. Similarly, the disruption of immune pathways differed between BE and EAC (Fig. 4d–f and Supplementary Table 6). Although interferon alpha and gamma responses were consistently downregulated in NP-BE and P-BE, both were upregulated in EAC, particularly in group 2 (Fig. 4b). Consistently, we observed a significant inverse correlation between expression fold changes of interferon gamma (Fig. 4g) and alpha (Fig. 4h) genes in BE and EAC groups 2 compared to 9p21 wild-type samples. Moreover, there was substantial overlap between altered genes in the two pathways (Fig. 4i), suggesting a comprehensive transcriptional reprogramming of interferon response. The most likely candidates for this reprogramming were again MTAP, given its recently reported ability to regulate the TME11, and IFNE, a type-1 interferon expressed in adult epithelia. Since the effect was most visible in group 2, which had LoF in both genes, and not in group 3, which had MTAP LoF and IFNE wild-type (Fig. 4a,b), the effect on interferon response might be due to IFNE loss.

CDKN2A LoF alone might instead be enough for the pervasive downregulation of keratinization genes given that these pathways were consistently dysregulated also in group 4 (Supplementary Table 6 and Fig. 4d–f).

Loss of IFNE reduces immune infiltration in BE, but not in EAC

To further investigate the opposite effect of IFNE on interferon alpha and gamma response in BE and EAC (Fig. 4g,h), we quantified the infiltration of 18 immune cell populations in NP-BEs, P-BEs and EACs from their bulk transcriptomic data. We then compared the abundance of immune infiltrates between each of the four 9p21 LoF groups (Fig. 4a) and the corresponding 9p21 wild-type samples.

Immune infiltrates were depleted in NP-BE groups 1 to 3 (Fig. 5a and Supplementary Table 7) and P-BE groups 1 and 2 as compared to 9p21 wild-type samples (Fig. 5b and Supplementary Table 7), where the impact of IFNE LoF was more appreciable. This again suggested that the immune depletion is a consequence of IFNE loss consistent with recent observations of a cold TME when IFNE10 or the whole IFN locus9 are lost in melanoma ovarian, or pancreatic cancers (Supplementary Table 8). However, the same studies also reported an increased infiltration of Treg cells, MDSCs and B cells (Supplementary Table 9) that we did not observe (Fig. 5a,b). The TME of group 4 (CDKN2A LoF only) was not significantly different to that of 9p21 wild-type samples in both NP-BE and P-BE, confirming that CDKN2A LoF does not directly interfere with the immune system.

Fig. 5: Impact of 9p21 gene loss on immune infiltration in BE and EAC.

a–c, Comparison of NESs of 18 immune populations between 9p21 LoF and wild-type samples in n = 22 patients with NP-BE (a), n = 108 patients with P-BE (b) and n = 337 patients with EAC (c), respectively. NES distributions were compared using a two-sided Wilcoxon’s rank sum test and corrected for multiple testing using the Benjamini–Hochberg method. Numbers of samples are reported in brackets. Immune populations with significant differences (FDR < 0.1) are outlined in red. d–f, Representative IMC images from group 2 (n = 4 patients, d), group 4 (n = 3 patients, e) and 9p21 wild-type (n = 3 patients, f) EACs showing the expression of 9p21 targeted proteins and mRNAs. Cadherin-1 and pan-keratin denote tumor. Arrows indicate examples of epithelial staining. Scale bar: 200 μm. g, Relative abundance of immune cells over all cells in 9p21 LoF and wild-type EACs. Samples in groups 2 and 4 were pooled together to form group 1 (n = 7 patients). Distributions were compared using a two-sided Wilcoxon rank sum test. h, Relative abundance of CD4+ cells over all CD3+ cells in 9p21 LoF and wild-type EACs. Distributions were compared using a two-sided Wilcoxon rank sum test. i, Median marker intensity across the T cell clusters at a clustering resolution of 0.5. j, UMAP map of 9750 T cells in n = 10 patients with EAC. Cells were grouped in 12 clusters based on the expression of six markers and colored according to the mean intensities of CD3 and CD4. The cluster enriched in group 1 is circled. Boxplots in g and h show first and third quartiles, whiskers extend to the lowest and highest value within the 1.5X interquartile range and the line indicates the median. Samples in groups 2 (n = 4 patients) and 4 (n = 3 patients) were pooled together to form group 1 (n = 7 patients). For 9p21 wt groups n = 3 patients are shown for all populations, except NK and dendritic cells where samples with no staining were removed. DCs, dendritic cells; TAMs, tumour-associated macrophages.

Source data

Unlike other cancer types (Supplementary Table 8) and BE (Fig. 5a,b), we did not observe any significant TME difference between 9p21 LoF and wild-type EACs (Fig. 5c and Supplementary Table 7). To investigate this at higher resolution, we performed high-dimensional imaging mass cytometry (IMC) on tissue sections representative of group 1, group 2, group 4 and 9p21 wild-type EACs (Supplementary Table 9). We used a panel of 26 antibodies targeting structural, immune and 9p21-encoded proteins as well as RNAscope probes against IFNE and IFNB1 mRNAs to increase the detection signal (Supplementary Table 10). We confirmed that group 2 lost the expression of all 9p21-encoded proteins in the tumor, whereas group 4 lost CDKN2A only compared to 9p21 wild-type EACs (Fig. 5d–f). Moreover, IFNE was the only interferon clearly expressed in EAC epithelium (Fig. 5d–f).

We performed single-cell segmentation of the IMC images to quantify T cells, NK cells, macrophages, dendritic cells, monocytic (M) and granulocytic (G) MDSCs, and neutrophils (Methods). We then compared the relative abundance of each immune population over all cells in each slide across EAC groups. We confirmed no significant difference in immune infiltration between 9p21 LoF and wild-type EACs, except for a borderline significant enrichment in dendritic cells in groups 1 and 2 (Fig. 5g). We further applied unsupervised clustering to T cells and macrophages, for which we had multiple markers (Supplementary Table 10), to test whether there was any difference in specific subpopulations. Again, we detected no major differences in any subpopulations of macrophages or T cells, except a borderline significant depletion of CD4+ T cells in groups 1 and 2 compared to 9p21 wild-type EAC (Fig. 5h–j). These results confirmed that, unlike BE, the loss of IFNE or any other 9p21 genes does not lead to any major difference in the TME of EAC.

CDKN2A LoF favors squamous to columnar epithelium transition

We observed a pervasive downregulation of processes responsible for terminal differentiation of keratinocytes, such as keratinization and formation of the cornified envelope, across all 9p21 LoF groups (Fig. 4d–f). In particular, P-BE and EAC groups 4 were associated with the downregulation of keratinization, suggesting that CD2KNA LoF alone was sufficient for triggering this process. To gain further mechanistic insights, we rebuilt the gene regulatory network linking CD2KNA LoF to keratinization in P-BE and EAC group 4 (Fig. 6a).

Fig. 6: Impact of CDKN2A LoF on epithelium differentiation in P-BE and EAC.