Context-dependent effects of CDKN2A and other 9p21 gene losses during the evolution of esophageal cancer

Abstract

CDKN2A is a tumor suppressor located in chromosome 9p21 and frequently lost in Barrett's esophagus (BE) and esophageal adenocarcinoma (EAC). How CDKN2A and other 9p21 gene co-deletions affect EAC evolution remains understudied. We explored the effects of 9p21 loss in EACs and cancer progressor and non-progressor BEs with matched genomic, transcriptomic and clinical data. Despite its cancer driver role, CDKN2A loss in BE prevents EAC initiation by counterselecting subsequent TP53 alterations. 9p21 gene co-deletions predict poor patient survival in EAC but not BE through context-dependent effects on cell cycle, oxidative phosphorylation and interferon response. Immune quantifications using bulk transcriptome, RNAscope and high-dimensional tissue imaging showed that IFNE loss reduces immune infiltration in BE, but not EAC. Mechanistically, CDKN2A loss suppresses the maintenance of squamous epithelium, contributing to a more aggressive phenotype. Our study demonstrates context-dependent roles of cancer genes during disease evolution, with consequences for cancer detection and patient management.

Fig. 1. CDKN2A LoF occurrence in BE and EAC.

a, Gene composition of chromosome 9p21 locus. b, Canonical EAC drivers damaged in at least 5% of EACs (n = 1,032 patients). All cell cycle regulators are reported in bold. c, Alterations in cell cycle regulators in BE and EAC. CDKN2A gene products (p14-ARF and p16-INK4a) regulate the cell cycle through the E2F genes. p14-ARF blocks MDM2 and TP53 degradation, which induces CDKN1A transcription. CDKN1A in turn inhibits the CCNE1/CDK2 complex ultimately blocking cell cycle through E2F1 inhibition. p16-INK4a directly inhibits the CCND/CDK6/CDK4 complex preventing RB1 phosphorylation. Unphosphorylated RB1 can bind E2F1, leading to cell cycle arrest. CDKN2A LoF favors cell cycle progression resulting in uncontrolled cell proliferation. Values within the circle represent the proportion of EACs, P-BEs and NP-BEs with at least one damaged cell cycle regulator (except TP53). d, Canonical EAC drivers damaged in at least 5% of P-BEs (n = 257 patients). e, Paired BE-EACs (n = 66 patients) with CDKN2A LoF. Clonally related alterations refer to either identical CDKN2A alterations in both lesions or CDKN2A alterations in BE that could further evolve in EAC. f, Canonical drivers damaged in at least 5% of NP-BEs (n = 99 patients). Alteration frequency of EAC canonical drivers in b, d and f is indicated in brackets. The alteration frequency of all EAC drivers in the three cohorts is available in Supplementary Table 2. FHCRC, Fred Hutchinson Cancer Research Center; LoF, loss of function; MSKCC, Memorial Sloan Kettering Cancer Center; NP-BE, non-progressor Barrett’s esophagus; EAC, esophageal adenocarcinoma; P-BE, progressor Barrett’s esophagus; SNV, single-nucleotide variant; TCGA, The Cancer Genome Atlas; UoC, University of Cambridge. Source data

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⁻⁴, 8 × 10⁻⁶ and 2 × 10⁻⁷, 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

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⁻⁴) (a), only homozygous deletions (P = 6 × 10⁻³) (b) and only LoF mutations (P = 3 × 10⁻³) (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

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 GSEA 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

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

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

a, Gene regulatory network linking CDKN2A LoF to the downregulation of keratinization genes through TF deregulations. b–d, Distributions of gene expression values of SOX15 in n = 17 patients with P-BE (P = 0.002) (b) and SOX15 (P = 0.04) (c) and TP63 (d) in n = 26 patients with EAC of group 4 and 9p21 wild type (31 P-BE and 184 patients with EAC, respectively). Distributions were compared using two-sided Wilcoxon’s rank sum test. FC and FDR from the differential gene expression analysis with DESeq2 (ref. ) are also shown. 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. e, Overlap between keratinization genes targeted by SOX15 and TP63 in P-BE and EAC. f,g, Preranked GSEA plots using as signature keratinization genes targeted by SOX15 in P-BE (f) and by SOX15 and TP63 in EAC (g). Genes were ranked from the most upregulated to the most downregulated in group 4 compared to 9p21 wild-type samples. For EAC, only the top 2,000 downregulated genes are shown. h–j, Correlation plots between keratinization GSEA NES and the gene expression values of SOX15 in P-BE (h) and SOX15 (i) and TP63 (j) in EAC. Coefficients and associated P values from two-sided Spearman’s correlation test are reported. k–n, Preranked GSEA plots using gene signatures for quiescent basal cells (k), proliferating basal cells (l), early suprabasal cells (m) and late suprabasal cells (n) in P-BE and EAC group 4. P values in (e–g and k–n) were estimated by permutation. ES, enrichment score; TF, transcription factor; GSEA, gene set enrichment analysis. Source data

Extended Data Fig. 1. Distribution of damaged genes in BE and EAC.

Damaged genes per sample in EAC WGS/WES cohorts (n = 779 patients) with any type of damaging alterations (A), homozygous deletions (B), gene amplifications (C), double hits (D) and damaging SNVs and indels (E). Number of damaged genes per sample in P-BE (n = 218 patients) and NP-BE (n = 63 patients) with any type of damaging alterations (F), homozygous deletions (G), gene amplifications (H), double hits (I) and damaging SNVs and indels (J). FHCRC, Fred Hutchinson Cancer Research Center; NP-BE, non-progressor Barrett’s esophagus; EAC, esophageal adenocarcinoma; P-BE, progressor Barrett’s esophagus; SNVs, single nucleotide variants; TCGA, The Cancer Genome Atlas; UoC, University of Cambridge. All boxplots show first and third quartiles, whiskers extend to 1.5X the interquartile lower and upper range and the line indicates the median. Source data

Extended Data Fig. 2. Effect of alterations in cell cycle regulators on EAC survival.

A. Kaplan-Meier survival curves of patients with EAC with wild type CDKN2A compared to those with heterozygous loss of CDKN2A. Heterozygous deletions could be inferred only in n = 779 patients with EAC with WGS or WES data. Kaplan-Meier survival curves of patients with EAC with alterations in CCND1 (p-value = 0.01) (B), TP53 (C), CDKN1A (D), CCNE1 (E) and MDM2 (F) compared to the corresponding wild type samples. Survival curves B-F were done using the whole cohort of 1032 EACs. GoF, gain-of-function; LoF, loss-of-function; EAC, esophageal adenocarcinoma; Het, Heterozygous deletions. Log-rank method was used to estimate the p-values. Source data

Extended Data Fig. 3. Expression of 9p21 genes in normal esophagus.

Transcript Per Million (TPM) expression values of the 26 9p21 genes in esophagus samples from n = 139 healthy individuals. Data were derived from the Genotype-Tissue Expression (GTEx) repository (https://gtexportal.org/). Source data

Extended Data Fig. 4. Workflow of causal gene network analysis.

The workflow to infer the causal gene network linking CDKN2A LoF to the downregulation of keratinization was divided into three steps: (A) identification and filtering of the keratinization-related gene modules and associated transcription factors (TFs) using cMonkey2 and ARACNE-AP (Step 1); (B) prediction of the causal models that link CDKN2A LoF to the dysregulation of keratinization genes through specific TFs using the network edge orienting (NEO) method^87,91 (Step 2); and (C) retention of the causal models with differential expression (FDR < 0.1) of TFs in group 4 compared to 9p21 wild type samples assessed using DESeq2 and significant positive correlation (R > 0.5 and FDR < 0.1) between TF expression and the GSEA NES score of the predicted targets in P-BEs and EACs (Step 3).