9p21 loss confers a cold tumor immune microenvironment and primary resistance to immune checkpoint therapy

Abstract

Immune checkpoint therapy (ICT) provides substantial clinical benefits to cancer patients, but a large proportion of cancers do not respond to ICT. To date, the genomic underpinnings of primary resistance to ICT remain elusive. Here, we performed immunogenomic analysis of data from TCGA and clinical trials of anti-PD-1/PD-L1 therapy, with a particular focus on homozygous deletion of 9p21.3 (9p21 loss), one of the most frequent genomic defects occurring in ~13% of all cancers. We demonstrate that 9p21 loss confers "cold" tumor-immune phenotypes, characterized by reduced abundance of tumor-infiltrating leukocytes (TILs), particularly, T/B/NK cells, altered spatial TILs patterns, diminished immune cell trafficking/activation, decreased rate of PD-L1 positivity, along with activation of immunosuppressive signaling. Notably, patients with 9p21 loss exhibited significantly lower response rates to ICT and worse outcomes, which were corroborated in eight ICT trials of >1,000 patients. Further, 9p21 loss synergizes with PD-L1/TMB for patient stratification. A "response score" was derived by incorporating 9p21 loss, PD-L1 expression and TMB levels in pre-treatment tumors, which outperforms PD-L1, TMB, and their combination in identifying patients with high likelihood of achieving sustained response from otherwise non-responders. Moreover, we describe potential druggable targets in 9p21-loss tumors, which could be exploited to design rational therapeutic interventions.

Fig. 1. 9p21 loss is frequently observed in human cancer and associated with significantly shortened survival.

a Schematic view of the chromosomal region 9p21.3 showing genes mapped to this focal region, their relative genomic locations, and frequency of 9p21 homozygous deletion (HD) observed in human cancer, based on data from the TCGA studies. b Pie charts showing the relative proportions of different types of somatic copy number variations (SCNAs) identified in MTAP and CDKN2A, respectively. The genomic data of 10,435 tumors from the TCGA program were analyzed. WT wildtype and diploid, LOH loss of heterozygosity (hemizygous deletion), HD homozygous deletion, Gain copy number gain or amplification. c The mRNA expression levels of MTAP (left) and CDKN2A (right) were markedly reduced in tumors with homozygous deletion of the genes. The Numbers of biologically independent samples were labeled on the violinplots. P values were calculated by two-sided Wilcoxon rank-sum test and adjusted for multiple testing. Box, median ± interquartile range; whiskers, 1.5× interquartile range. ***P value < 0.001. Exact P values were P < 2 × 10⁻¹⁶ for all comparisons. d (left) The relationship of different types of SCNAs between MTAP and CDKN2A and their relative frequencies (right). Mut mutation. e The landscape of 9p21 SCNAs across TCGA cohorts. The colors are the same as shown in the panel d. (bottom) Histogram showing the fraction of different types of 9p21 SCNAs (as defined in panel d) across TCGA cancer types (see Supplementary Data 1 for a complete list). (top) Line plot showing the fraction of MTAP and CDKN2A specific events and co-deletions. f Representative tumor types demonstrating great variation in the frequencies of 9p21 loss across previously defined molecular subtypes (see Supplementary Data 2 for the abbreviations of disease codes). P values were calculated by two-tailed Fisher’s exact tests. g The prognostic significance of 9p21 loss at pan-cancer level in TCGA cohorts. A total of 10,283 patients with available survival data were included in survival analysis. The line colors are the same as shown in the panel d. Log-Rank P values and the median overall survival time (in months) are shown. mo, months. h Univariate Cox regression analysis of 9p21 loss for overall survival across 12 TCGA cohorts with frequent 9p21 loss (>10%, see Supplementary Data 3). Numbers within the parentheses indicate the sample size. P values were calculated by Cox proportional hazards (PH) regression model. Error bars indicate the estimated 95% confidence interval of the hazard ratio. i Representative examples showing that 9p21 loss is associated with significantly shortened overall survival in individual cancer cohorts. The cancer type, molecular subtype, sample size, and Log-Rank P values are labeled on each plot. P values were calculated by two-sided Log-rank test.

Fig. 2. 9p21 loss is associated with ‘cold’ tumor-immune phenotypes.

a Schema showing the patterns of spatial distribution of TILs defined by a previous TCGA study by Saltz et al. b Gradient changes in the spatial TILs patterns among 9p21-WT tumors, 9p21-LOH tumors, and 9p21-loss tumors were observed, which corresponded to progressive copy number loss of 9p21. The plots of three representative cancer types are shown (see more details in Supplementary Fig. 7). The FDR q-values did not reach significance level at 0.05. c 9p21 loss in shaping the immune cell abundance and cell composition in tumor microenvironment. Immune deconvolution was performed by applying MCP-counter to the bulk RNA-seq data, similarly as described in our recent studies^,. The data is shown for 12 TCGA cohorts (14 molecular subtypes) with frequent 9p21 loss (>10%, see Supplementary Data 3). The bubble plot is drawn using computed log2-transformed fold change (9p21-Loss vs. 9p21-WT) and adjusted p-values (FDR q-value). The size of the bubble indicates statistical difference, the bigger the more significant. The color of the bubble indicates change in the immune cell abundance in 9p21-loss tumors (vs. 9p21-WT), with blue denotes depletion and red denotes enrichment. d Box plots of representative examples selected from the panel c (see more details in Supplementary Fig. 8). P values were calculated by two-sided Wilcoxon rank-sum test. Number of samples: B cells in HNSC_HPV-: WT (n = 57); Loss (n = 135); T cells in SKCM: WT (n = 71); Loss (n = 112); CD8 T cells in HNSC_HPV-: WT (n = 57); Loss (n = 135); CD8 T cells in PAAD: WT (n = 40); Loss (n = 44); CTLs in PAAD: WT (n = 42); Loss (n = 44); CTLs in STAD: WT (n = 163); Loss (n = 46). Box, median ± interquartile range; whiskers, 1.5× interquartile range. e The richness and diversity of T-cell receptor (TCR) repertoire was decreased in tumors with 9p21 loss in multiple TCGA cohorts. The diversity of TCR repertoire is indicated by the Shannon entropy. The color of the bars indicates the significance level of changes in 9p21-loss tumors (vs. 9p21-WT). f Changes in immunomodulatory gene expression in 9p21-loss tumors in comparison with 9p21-WT tumors. A list of 28 immunomodulatory genes (see a full list in Supplementary Data 9) were analyzed and the most significant ones are shown. The color of the bubble corresponds to Log2 fold change in gene expression levels in 9p21-loss tumors (vs. 9p21-WT), with blue denotes decrease and red denotes increase in 9p21-loss tumors.

Fig. 3. 9p21 loss is associated with immune resistance to anti-PD-1/L1 monotherapy in solid tumors.

a–c the MDA (MD Anderson Cancer Center) solid tumor cohort (n = 94 patients). a Schematic view of the information collection and analysis flow. b 9p21 loss is associated with lack of response to anti-PD-1/L1 monotherapy in the ICT “responsive” tumor cohort. The response rates (percentages of CR/PR) were compared between the two groups. c 9p21 loss is associated with disease progression following anti-PD-1/L1 monotherapy in the ICT “refractory” tumor cohort. The progression rates (percentages of PD) were compared between the two groups. P values were calculated by two-tailed Fisher’s exact tests. d, e the high-risk resectable melanoma cohort from Helmink et al. hi_hi, tumors with mRNA expression levels of both CDKN2A and MTAP above the group median and lo_lo, tumors with expression levels of both genes below the group median. d Comparison of the response rates (percentages of CR/PR) to ICT between the hi_hi and lo_lo groups. P values were calculated with two-sided Fisher-exact test. e Waterfall plot showing the RECIST response calculated based on the percentage of change in tumor volume relative to baseline. P value was calculated using the two-sided Mann–Whitney U test. f The combined melanoma cohort from 4 studies (see Table 1 and Supplementary Data 10 for details). The response rates were compared between the hi_hi and lo_lo groups for all patients together (left), and in individual patient subpopulations receiving nivolumab (middle) and pembrolizumab (right), respectively. P values were calculated using the two-tailed Fisher’s Exact tests.

Fig. 4. 9p21 loss is associated with immune resistance to anti-PD-1/L1 monotherapy in large metastatic urothelial cancer (mUC) and advanced non-small-cell lung cancer (NSCLC) cohorts.

a–d the MDA mUC cohort. A total of 86 mUC patients who received either pembrolizumab or atezolizumab monotherapy were included and 80 patients with available response and MTAP IHC data were taken into subsequent analyses. Samples were collected prior to ICT. a Schematic view of the information collection and analysis flow. b Decreased trend of PD-L1 stain positivity in MTAP-negative tumors. Colors in this plot indicates the four categories of PD-L1 IHC staining results. c MTAP loss is associated with primary resistance to ICT and disease progression following pembrolizumab or atezolizumab monotherapy. P values were calculated using the two-tailed Fisher’s Exact tests by comparing the rates of disease progression (percentages of PD) between two groups. d MTAP loss is associated with worse progression-free survival (PFS) and disease-specific survival (DSS) in mUC patients received pembrolizumab or atezolizumab monotherapy. e–l the MSK NSCLC cohort from Rizvi et al.. e Schematic view of the information collection and analysis flow. A total of 151 LUAD patients received PD-1/L1 as monotherapy with available genomic and response data were included in subsequent analyses. f Decreased trend of PD-L1 positivity in tumors with 9p21 loss. Colors in this plot indicates the categorized PD-L1 IHC staining results. g 9p21 loss is associated with a lower rate of DCB (response defined and shorter PFS (h). DCB, durable clinical benefit, defined as complete/partial response or stable disease that lasted >6 months by the original study (the detailed classification of CR, PR, SD, PD were not available). NDB no durable benefit. Integration of 9p21 status with TMB (i, j) or PD-L1 expression (k, l) in patient stratification for response and PFS. TMB tumor mutation burden. PD-L1 expression was measured by immunohistochemistry staining by the original study. P value in panel H was calculated with two-sided Log-rank test. P values in panels g, i and k were calculated by two-tailed Fisher’s exact tests.

Fig. 5. Validation of the translational impact of 9p21 loss on ICT in large-scale metastatic urothelial cancer (mUC) cohort.

Patients were from the IMvigor210 trial investigating Atezolizumab (anti-PD-L1) blockade in the mUC cohort (n = 298 patients). Pre-treatment samples were collected for bulk RNA-seq and immune profiling and the data was downloaded from a published study from Mariathasan et al.. a The proportions of immune (left) and tumor (middle) cells that were positive for PD-L1 staining (by SP142 immunohistochemistry) were significantly lower in tumors with decreased CDKN2A expression (top and bottom quantiles, Q4 vs. Q1), and the fraction of “inflamed” immune phenotype (right) was also significantly lower in low (Q1) than those with high (Q4) CDKN2A expression. The immune phenotypes were defined by CD8 IHC staining by the original study. P values were calculated by two-tailed Fisher’s exact tests. b (left, middle) The CDKN2A and MTAP expression levels, respectively, were stratified into decreasing quantiles, and the response rates (percentages of CR/PR) decreased significantly with decreasing MTAP/CDKN2A expression. (right) CDKN2A and MTAP co-expression patterns can better stratify patients for response (hi_hi, tumors with mRNA expression levels of both CDKN2A and MTAP above the group median, n = 124, and lo_lo, tumors with expression levels of both genes below the group median, n = 127), and overall survival (c) following PD-L1 blockade by Atezolizumab. P values in panel b were calculated by two-tailed Fisher’s exact tests. P value in panels c was calculated with two-sided Log-rank test. d Multivariable Cox regression analysis showing that 9p21 loss was a strong prognosticator of short survival, independent of other variables listed. Cox proportional hazards (PH) regression model was used to calculate the Hazard Ratio (HR), the 95% confidence interval (95%CI) and P values. Error bars indicate the estimated 95%CI of the HR. Number of samples: CDKN2A_MTAP, lo_lo (n = 127), hi_hi (n = 124); Tumor-cell PD-L1, High (n = 42), Low (n = 255); Immune cell cell PD-L1, High (n = 102), Low (n = 195); TMB, High (n = 120), Low (n = 114); Sex, Male (n = 233), Female (n = 65); Tobacco use, Smoker (n = 32), and non-smoker (n = 266). e 9p21 status can compensate PD-L1 expression and TMB in identifying the responders and non-responders to Atezolizumab and showed significant correlates with survival (f). The cut off of PD-L1 expression was 5% as suggested by the original study, and the median value of TMB was used to split patients into TMB-high and TMB-low groups. Patients without PD-L1 IHC data (n = 1) and those without TMB data (n = 64) were excluded from corresponding analysis. Response scores were calculated by incorporating three factors (9p21, PD-L1 expression on immune cells, TMB), which stratified patients into three groups, with high^,, intermediate^,, and low(0) response score. Log-Rank P values and the median overall survival time (in months) are shown. mo months. P values in panel e, f were calculated by two-tailed Fisher’s exact tests.

Fig. 6. Therapeutic vulnerability and potential immunotherapy targets in tumors with 9p21 loss.

a–c Identification of potential immunotherapy targets in the mUC cohort from Mariathasan et al.. a Differentially expressed immune-related genes in the CNKN2A_MTAP: lo_lo tumors. A curated list of ~500 genes (including known and emerging viable immunomodulatory targets and other druggable targets of cancer and cytokines, see Supplementary Data 17 for the complete list) were analyzed the most significant genes that upregulated in the lo_lo group (except CD274 which was downregulated) were labeled on the plot. Two vertical lines indicate gene expression fold change (lo_lo vs. hi_hi) >1.2 and <−1.2, respectively, and the horizontal line indicates the adjusted P value (FDR q-value) of 0.05. P values were calculated by two-sided Wilcoxon rank-sum test. The color of the dot represents the FDR (q-value) levels. b Spearman correlation analysis identified potential immunotherapy targets that were reversely correlated with CDKN2A/MTAP expression, i.e. upregulated in tumors with low CDKN2A/MTAP expression. The Spearman correlation ecoefficiency is shown on the x axis and the bars are color coded by FDR q-value. Two vertical lines indicate Spearman’s ρ < −0.2 and <−0.4, respectively. The color of the bar represents the FDR (q-value) levels. c Box plots showing representative genes displayed in panels a and b. The expression levels were compared in the pre-treatment tumors between the lo_lo and hi_hi groups and stratified by patient’s response status (SD and PD). Sample size: SD, hi_hi (n = 24), lo_lo (n = 18); PD, hi_hi (n = 57), and lo_lo (n = 67). P values were calculated by two-sided Wilcoxon rank-sum test. Box, median ± interquartile range; whiskers, 1.5× interquartile range. d–f Identification of potential immunotherapy targets in the TCGA cohorts. d Spearman correlation of gene expression with CDKN2A/MTAP across 12 TCGA cohorts (14 molecular subtypes) with frequent 9p21 loss (>10%, see Supplementary Data 3). The size of the bubble represents the correlation levels. The color of the bubble represents the FDR levels. Red: positive correlation. Blue: negative correlation. e Scatter plots showing representative genes displayed in the panel d. The cancer type and molecular subtype, Spearman correlation ecoefficiency and FDR q-value are labeled on each plot. Error bands indicate the estimated interval of correlation level. f Box plot showing VTCN1 (B7-H4) expression between the 9p21-loss and 9p21-WT groups. Sample size: STAD_All, WT (n = 171), Loss (n = 50); STAD_CIN, WT (n = 45), Loss (n = 39); ESCA_All, WT (n = 32), Loss (n = 64); ESCA_ESCC, WT (n = 12), Loss (n = 52); LUAD_All, WT (n = 162), and Loss (n = 87). Box, median ± interquartile range; whiskers, 1.5× interquartile range. g Schema summaries the immunological modulation of 9p21 to the TME and potential immunotherapy targets identified in this study.