Molecular subtypes of oropharyngeal cancer show distinct immune microenvironment related with immune checkpoint blockade response

Abstract

Background: Oropharyngeal cancer (OPC) exhibits diverse immunological properties; however, their implications for immunotherapy are unknown.

Methods: We analysed 37 surgically resected and nine recurrent or metastatic anti-programmed death-1 (PD-1)/programmed death-ligand 1 (PD-L1)-treated OPC tumours. OPCs were classified into immune-rich (IR), mesenchymal (MS) and xenobiotic (XB) subtypes based on RNA-sequencing data.

Results: All IR type tumours were human papillomavirus (HPV) positive, most XB types were HPV negative, and MS types showed mixed HPV status. The IR type showed an enriched T cell exhaustion signature with PD-1⁺ CD8⁺ T cells and type I macrophages infiltrating the tumour nest on multiplex immunohistochemistry. The MS type showed an exclusion of CD8⁺ T cells from the tumour nest and high MS and tumour growth factor-β signatures. The XB type showed scant CD8⁺ T cell infiltration and focal CD73 expression. The IR type was associated with a favourable response signature during anti-PD-1/PD-L1 therapy and showed a high APOBEC mutation signature, whereas the MS and XB types showed resistance signature upregulation. Among anti-PD-1/PD-L1-treated OPC patients, the IR type showed a favourable clinical response (3/4 patients), whereas the XB type showed early progression (3/3 patients).

Conclusion: Our analysis classified OPCs into three subtypes with distinct immune microenvironments that are potentially related to the response to anti-PD-1/PD-L1 therapy.

Fig. 1. Three molecular subtypes of OPC identified by RNA-seq analysis of tumours show different survival outcomes.

a The two-dimensional t-SNE plot of RNA-seq data from 37 surgically resected OPC tumours (YOPC) with p16 IHC status. The k-means clustering analysis classifies the tumours into three subtypes: immune-rich (IR), mesenchymal (MS) and classic (CL) type. The patient number (YOPC) corresponds to the plots. b Comparison of smoking doses of patients among three subtypes. The smoking dose in pack-years (PYR) was compared using Kruskal–Wallis test with post hoc analysis using Dunn’s multiple comparison test. The difference between XB and IR type OPCs was statistically significant (p < 0.05) in the post hoc analysis and indicated with the asterisk (*p < 0.05). c The t-SNE plotting for OPC tumours of the combined cohort, which consists of 37 OPC tumours from patients at the Yonsei Cancer Center (YOPC) and 33 HNSC TCGA tumours of oropharyngeal origin (the soft palate, the base of tongue, the tonsils and the side and back wall of the throat). The subtype of TCGA OPC tumours was determined by k-means clustering. d The Kaplan–Meier curves for overall survival of OPC patients in the combined cohort. The patient survival was compared using log-rank test. The difference between XB and IR type OPCs was statistically significant (p < 0.167) after the Bonferroni post hoc correction and indicated with the asterisk (*). e The heatmap of significantly altered genes among three subtypes (log₂ fold change ≥2 and p value < 0.05). The CL-, IR- and MS-specific gene signatures were selected according to their expression in each subtype. f The three-dimensional plotting and k-means clustering analysis using gene set enrichment score of each subtype signature derived from e calculated by GSVA algorithm in 518 HNSC TCGA patients (upper panel). The Kaplan–Meier curves for overall survival of 518 HNSC TCGA patients according to subtype classification (lower panel). The difference between XB type vs. IR type OPCs was statistically significant (p < 0.167) after the Bonferroni post hoc correction. Throughout the figures, data are mean and SEM.

Fig. 2. Characterisation of subtype-specific gene signatures.

a The differentially expressed genes among the subtypes (log fold change ≥1 and p value < 0.05) were categorised into six clusters: xenobiotic metabolism (C1), adaptive immune response (C2), cell cycle-related (C3), smooth muscle contraction and cell adhesion (C4), keratinisation (C5) and membrane potential regulation (C6) genes in 37 surgically resected OPC tumours (YOPC). b Subtype-specific expression of six cluster genes. The average z-score expression of each gene is shown by a grey line, and the average z-score expression of all genes is shown by a red line. c The network analysis of differentially expressed genes via pair-wise comparison of each subtype using Enrichment map software in 37 YOPC tumours. d The heatmap showing enrichment score of hallmark gene signatures (H) in MSigDB on GSVA. Red arrows indicate the gene signatures related to metabolism in 37 YOPC tumours. e Tumour-to-liver uptake ratio (TLR) from F-18 FDG PET-CT scan of OPC tumours according to their subtype. The pre-operative PET-CT scan was available in 31 out of 37 OPC tumours, and the representative PET/CT images of subtypes are shown. The TLR values were compared using one-way ANOVA with Bonferroni post hoc test. The difference between XB type and MS type OPCs was statistically insignificant in the post hoc analysis. n.s., Not significant.

Fig. 3. Distinct tumour immune microenvironment of OPC subtypes revealed by RNA-seq.

a Significantly altered immune-related genes (log₂ fold change ≥1 and p value < 0.05 in GO term immune response: GO:0006955) among three subtypes. The representative genes related to T cell adaptive immune response and macrophage/granulocyte activation are indicated. b GSEA of expression profile comparing IR type vs. XB and MS types using T cell exhaustion (LAYN and Mel75) and myeloid cell activation (GO:0002274: myeloid leucocyte activation and GO:0043030: regulation of macrophage activation) gene signatures. c Immune cell composition analysis using the CIBERSORT deconvolution method. The heatmap showing the relative fraction of various types of immune cells in OPC tumours (the upper panel). The dot plot showing ratio of CD8⁺ T cells, M1 macrophages and M2 macrophages in OPC tumours (the lower panels). The proportion of the cells was compared using one-way ANOVA with Bonferroni post hoc test. The statistically significant results per post hoc analysis are indicated with asterisks (*p < 0.05). d The waterfall plot of TIDE prediction score of 37 OPC tumours. The molecular subtype is indicated by a different colour. A low TIDE score means a high probability of response to immune checkpoint blockade therapy. e The heatmap of the gene set enrichment score of 19 IPRES gene sets calculated by GSVA in 37 OPC patients.

Fig. 4. Distinct immune microenvironment of OPC subtypes revealed by multiplex IHC.

a The representative multiplex IHC image of OPC tumour. The representative tumour nest and stroma sites were selected, and the density of CD8⁺ T cells and CD68⁺ macrophages, as well as the proportion of CD73⁺ tumour cells were counted both in tumour nests and stroma. b The representative multiplex IHC images of subtypes of OPC tumours. c The dot plots of immune cell densities of OPC tumours, total CD8⁺ T cells, PD-1⁺CD8⁺ T cells and CD68⁺ macrophages both in tumour nests and stroma. d The ratio of CD8⁺ T cell density in tumour nest to CD8⁺ T cell density in stroma in the three subtypes (left panel). The f-TBRS score of the three subtypes in 37 YOPC tumours (right panel). e The representative CD73 IHC image of two XB type tumours (YOPC8 and YOPC32, upper panel). The dot plot of the proportion of CD73⁺ tumour cells in the subtypes (lower panel). The cell densities (c) and f-TBRS (d) were compared using one-way ANOVA with Bonferroni post hoc test. The statistically significant results per post hoc analysis are indicated with asterisks (*p < 0.05).

Fig. 5. Analysis of subtype-specific genetic properties by next-generation sequencing of OPC tumour DNA.

a The oncomap showing genetic alterations of OPC tumours. The tumour-blood-matched targeted panel sequencing was performed in 13 out 37 OPC patients with available samples. The tumour subtype is indicated at the bottom of the map and the tumour mutation burden (TMB, non-synonymous mutations/Mb) was plotted. b The plot showing correlation between tumour mutation burden and TIDE score in OPC tumours. c Pie charts showing the proportion of COSMIC mutational signatures in XB type, MS type and IR type of OPC tumours. b The dot plot comparing the expression level of APOBEC family genes, APOBEC3B, APOBEC3D and APOBEC3G, in the three subtypes. The normalised expression level on RNA-seq was compared by unpaired t test. The expression values were compared using one-way ANOVA with Bonferroni post hoc test. The statistically significant results per post hoc analysis are indicated with asterisks (*p < 0.05).

Fig. 6. The molecular subtype predicts response to anti-PD-1/PD-L1 therapy in OPC patients.

a The patient list and clinical characteristics of nine recurrent or metastatic OPC patients receiving anti-PD-1/PD-L1 therapy (YOPD). *The best response was determined by RECIST 1.1 criteria. **Patient still undergoing anti-PD-1/PD-L1 therapy. b The t-SNE plot for nine YOPD tumours combined with 37 YOPC tumours. The subtype of YOPD tumours was determined by distance from 37 YOPC tumours on the t-SNE plot. c Spider plot showing changes in target lesion diameters on anti-PD-1/PD-L1 therapy in nine OPC tumours according to subtype. d Representative CT scan images of IR type and XB type tumours on anti-PD-1/PD-L1 therapy.