We used five cohorts, namely, the Morris cohort (Samstein et al., 2019), the Allen cohort (Miao et al., 2018), the Ho cohort (Morris et al., 2017), the Berger cohort (Zehir et al., 2017), and TCGA cohort. The inclusion criteria for r/m HNSCC is shown in Figure 1. The Morris cohort finally included 120 patients, and the Allen cohort included 12 patients (Supplemental Table S1). All patients were diagnosed with r/m HNSCC, not including nasopharyngeal carcinoma. Those 132 patients had undergone immunotherapy with anti-programmed cell death protein 1 (PD1)/PD-L1 agents alone or combined with CTLA-4 inhibitors. Only patients for whom whole-exome sequencing (WES) and survival information were available were analyzed for this study. The Ho cohort included 115 r/m HNSCC patients, and the Berger cohort included 49 r/m HNSCC patients who did not receive immunotherapy (Supplemental Table S2, S3). TCGA cohort included 488 patients with HNSCC with few r/m HNSCC patients included (Supplemental Table S4). No patient from TCGA cohort had received immunotherapy. All relevant information from the five cohorts was downloaded from cbioportal.org, including mutation landscape data, WES, and RNA sequencing (RNA-seq) data, and survival information.

TMB was calculated as total variants divided by the total callable bases covered by the MSK-IMPACT or WES panel in megabases (Mb) for both the intronic and exonic sequences (total variants/Mb) or for exonic and splice site sequence (nonsynonymous mutations/Mb). TMB score of the Morris cohort, the Allen cohort, the Ho cohort, and the Berger cohort was downloaded from cbioportal.org. TMB score of TCGA cohort was downloaded from the literature (Thorsson et al., 2018). The mean coverage of the WES was 38 Mb. Variant allele frequencies cutoffs can vary from 0.5 to 10%, with lower thresholds increasing the risk of including false positives arising from sequencing artifacts such as C-to-T transitions introduced by formalin fixation (Chan et al., 2019). Thus, we selected variants with an allelic fraction of at least 10%.

We used the CIBERSORT web portal (https://cibersort.stanford.edu/) to evaluate the infiltration status of 22 immune cell types for TCGA cohort. Then, we analyzed the differences in immune cell infiltration between the samples in different groups.

RNA-seq data were downloaded from UCSC Xena (https://xenabrowser.net/datapages). RNA-seq count data normalized to FPKM (fragments per kilobase per million) values were used as input for single-sample gene set enrichment analysis (ssGSEA). ssGSEA was performed for select signatures using the gene set variation analysis package with default settings. The association of each signature’s ssGSEA scores with TP-PR was tested using Spearman’s correlation test.

Univariate and multivariate Cox regression analyses were performed to analyze the gene mutations significantly associated with OS. TMB and the gene mutations significantly affecting OS were used to construct TP-PR score. The TP-PR score of each sample was calculated by multiplying the values of certain indexes by their weight in the Cox model and then adding them together. The prognostic index TP-PR (formula: TP-PR = TMB * 1.0803 + PIK3CA * 0.654 − TP53 * 0.5581 − ROS1 * 1.2342) was then constructed according to the expression of those three genes. It was scored 1 if TMB ≥10 mut/Mb or if the gene was mutated and was scored 0 if TMB <10 mut/Mb or if the gene was wide-type. The final score was calculated as the TP-PR score. Regression coefficients were generated in the process of multifactor Cox regression analysis with R software (version 4.1.1). The prognostic power of the TP-PR score was evaluated by Kaplan–Meier (K-M) survival curves with log-rank tests with the ICI cohort.

The K-M estimator was used to describe the distribution of survival time. The log-rank Mantel–Cox test was used to compare survival between different strata. Statistical analyses were performed using GraphPad software, version 8.0 (GraphPad Software Inc., La Jolla, CA, USA). All data points and statistical analyses represent individual patients.

We collected the data of 132 r/m HNSCC patients receiving immunotherapy from cBioPortal (https://www.cbioportal.org/study/summary?id=tmb_mskcc_2018; https://www.cbioportal.org/study/summary?id=mixed_allen_2018). The inclusion criteria for r/m HNSCC are shown in Figure 1. Detailed characteristics including age, gender, sample type, and administered immunotherapy drugs are shown in Table 1. We characterized patients as having TMB-H or TMB-L on the basis of the threshold of TMB ≥10 mut/Mb. Our analysis showed that the median OS of patients with TMB-H was 2.5 times as long as that of the patients with TMB-L (25 vs. 10 months), with an area under the curve (AUC) value of 0.56 (95% CI 0.68–0.44) (Figures 2A,B). In other words, immunotherapy was more effective in r/m HNSCC patients with TMB-H (TMB ≥10 mut/Mb). In addition, we collected data from two other cohorts that included r/m HNSCC patients receiving no immunotherapy (Figures 2C,D). Results showed that the difference in median OS between the patients receiving no immunotherapy with TMB-H and those with TMB-L was significant in the Ho cohort (Figure 2C, p = 0.0043) and not significant in the Berger cohort (Figure 2D, p = 0.55). It reminded us that the significance of TMB stratification in non-ICI cohorts of HNSCC was still uncertain. Thus, we need to look for a comprehensive prognostic index other than TMB alone.

The mutational landscape showed that the commonly mutated genes in r/m HNSCC patients from ICI-treated cohorts were TP53, TERT, PIK3CA, NOTCH1, FAT1, KMT2D, and CDKN2A (Supplementary Figure S1A). The flowchart of gene selection is shown in Supplementary Figure S1B. Twenty-seven genes with mutation frequency of no less than 5% were obtained from the results of mutational landscape. The correlation between those gene mutations and the OS was analyzed by K-M univariate analysis (Figures 3A–C; Supplementary Figure S2A–X). K-M univariate analysis showed that TP53, PIK3CA, ARID1A, ROS1, ATR, and ASXL1 might be related to the OS (p < 0.1) for r/m HNSCC patients who received immunotherapy. Then the six mutation genes and TMB were included for multivariate Cox regression analysis. Finally, the three genes, PIK3CA, TP53, and ROS1, were found to be associated with the OS (p < 0.05) (Figure 4). Among them, TMB-H or PIK3CA mutation was positively correlated with OS, while TP53 or ROS1 mutation was negatively correlated with OS. The prognostic index TP-PR (formula: TP-PR = TMB * 1.0803 + PIK3CA * 0.654 − TP53 * 0.5581 − ROS1 * 1.2342) was then constructed according to the expression of those three genes. It was scored 1 if TMB ≥10 mut/Mb or if the gene was mutated and was scored 0 if TMB <10 mut/Mb or if the gene was wide-type. The final score was calculated as the TP-PR score.

The data sets of r/m HNSCC patients who had received immunotherapy were divided into two groups: high TP-PR (TP-PR >0) and low TP-PR (TP-PR ≤0). We found that the high TP-PR group had better median OS (25 vs. 8 months) than the low TP-PR group, with a higher AUC value (0.77, 95% CI 0.86–0.68) than TMB stratification (Figures 5A,B). Moreover, the high TP-PR group screened out more people than the TMB-H group. In other words, TP-PR screened out more people who might benefit from immunotherapy than TMB. We further contrasted the mutational landscape of TP-PR high population and TMB-H population (Figure 6). We found that PIK3CA mutation rate was higher in the TP-PR high population, whereas TP53 and especially ROS1 (0 vs. 14%) mutation rates were higher in the TMB-H population. That was in accordance with our conclusion that PIK3CA mutation was positively correlated with OS, while TP53 or ROS1 mutation was negatively correlated with OS. And it might be the answer why TP-PR screened out more people who might benefit from immunotherapy than TMB alone. Gene mutation analysis was also performed between the high and low TP-PR groups. It was confirmed that the expressions of PIK3CA, TP53, and ROS1 were different in the high and low TP-PR groups; especially, those of PIK3CA and ROS1 were significantly different (data not shown). It supported that PIK3CA, TP53, and ROS1 were related to the prognosis of r/m HNSCC patients.

TCGA cohort included 488 patients with HNSCC who had not received immunotherapy. RNA-seq analysis of TCGA cohort showed that CD8⁺ T cells and activated memory CD4⁺ T-cell infiltration in the tissues of the patients with high TP-PR were higher than those in the patients with low TP-PR, and the differences were statistically significant (Figure 7A). Infiltration of other immune cells was not significantly different between the two groups. It was suggested that high TP-PR might enhance CD8⁺ T cells and activated memory CD4⁺ T-cell infiltration, without affecting the mRNA expression of immunotherapy-related targets like LAG3, CTLA4, IDO1, and HAVCR2 (Figure 7B). But the PD-L1 (p = 0.047) and TIGIT (p = 0.026) mRNA expressions were higher in the high TP-PR group than in the low TP-PR group, indicating that high TP-PR might promote the mRNA expression of PD-L1 and TIGIT.

Then we tried to find biological processes (BPs) and signaling pathways that might be relevant by gene ontology and Kyoto Encyclopedia of Genes and Genomes analyses. Figure 8 shows the 15 most relevant BPs or signaling pathways. It is suggested that the differential gene expression in the high and low TP-PR groups in TCGA cohort mainly influenced BP related to development and metabolism, and metabolism-related signaling pathways, such as the pathways for fatty acid metabolism, glycerolipid metabolism, and nitrogen metabolism.

RNA-seq analysis also showed that CD8⁺ T cells and resting memory CD4⁺ T-cell infiltration in patients with TMB-H were higher than those in patients with TMB-L, with a statistically significant difference (Supplementary Figure S3A). But there was no difference in the mRNA expression of PD-L1, TIGIT, LAG3, CTLA4, IDO1, and HAVCR2 between the TMB-H and TMB-L groups (Supplementary Figure S3B). Compared with analysis stratified by TMB alone, patients with TP-PR stratification showed a difference in the mRNA expression of PD-L1 and TIGIT, which might influence the immune microenvironment. The differential gene expression in TMB-H and TMB-L groups in TCGA cohort also mainly influenced metabolism-related signaling pathways (Supplementary Figure S4).