Introduction

Hypopharyngeal squamous cell carcinoma (HPSCC) is a clinical challenge. It shows the poorest outcome of all head and neck cancer sites, partly due to frequently being diagnosed in late stages1. The total number of cases/year is around 80,000 worldwide2 and the main risk factors are smoking and alcohol, although, a small subset of HPSCC is caused by human papillomavirus (HPV)3,4,5,6,7. HPSCC caused by HPV, analogous with oropharyngeal squamous cell carcinoma, seems to have a much better clinical outcome than HPV-negative HPSCC3,4.

Patients with early stage as well as locoregionally advanced but curable HPSCC are treated with radiotherapy (RT), chemoradiotherapy (CRT), local surgery or combinations of these1. Additionally, there are several current clinical trials evaluating immunotherapy with immune checkpoint inhibitors as an optional therapy in HPSCC, however the majority of those are in the recurrent disease/metastatic setting. Several studies, including a recent meta-analysis, have demonstrated better survival for patients undergoing surgery together with RT, as compared to RT/CRT alone, but larger randomised control trials are needed to conclude this1,8. Hypopharyngeal surgery is demanding for the patients, who often are in a poor condition, and surgery frequently results in severe side-effects. Since HPSCC often is diagnosed in advanced stages, where surgery, when performed, is extensive, there has been a shift towards oncological treatment during the last decades1,8. At the Karolinska University Hospital, oncological treatment has dominated HPSCC treatment since the year 2000. However, if survival could increase with more patients recommended for surgery, there is a need to assess/identify which patients would benefit and there is a need for markers or gene signatures identifying tumors who do not respond satisfactory to oncological treatment alone. Predictive and prognostic biomarkers are thus needed to optimize treatment for patients with limited as well as locoregionally advanced tumors, but also in the metastatic setting where the objective treatment response rate is approximately 20%, it is important to avoid unnecessary toxicity in patients who will not respond.

The aim of this investigation was therefore to find prognostic/predictive biomarkers enabling identification of patients with poor prognosis after oncological treatment. For this purpose, mRNA expression in HPSCC biopsies was evaluated in relation to patient survival using the PanCancer IO 360 panel. This was followed by gene set enrichment analysis, including studying the expression of an 18-gene signature, earlier identified as related to response to PD-L1 inhibitors for several cancer types9,10,11,12,13. Furthermore, the expression of two putative biomarkers identified during the initial analysis, S100A12 and HLA class I, was validated to correlate with survival through immunohistochemistry (IHC).

Results

Gene expression analysis on the PanCancer IO 360 panel

The expression levels of 750 mRNA transcripts in 23 microdissected HPSCC samples were successfully quantified on the PanCancer IO 360 panel. The included tumors, 3 HPV-positive and 20 HPV-negative, were selected based on completed curative treatment and small tumor size. 10 of the HPV-negative samples were from patients having 3-year relapse-free survival and 10 from patients having a relapse or death by HPSCC within 3 years. To identify gene expression related to patient survival when only receiving radiotherapy or chemoradiotherapy, tumors from patients receiving local surgery were excluded, as were patients who died of causes other than HPSCC within 3 years (Table 1).

Table 1 Patient and tumor characteristics.
Full size table

Hierarchical clustering separates samples by prognosis

Upon initial hierarchical clustering, samples could be divided into four clusters; HPV-positive samples formed one cluster, while HPV-negative samples divided into three groups with high, intermediate and low proportions of survivors, respectively (Fig. 1A).

Figure 1
figure 1

(A) Heatmap of all samples. Genes are rows and samples columns. (NED = No evidence of disease, DOD = dead of disease). (B) Volcano plot of transcripts differing between HPV-positive and HPV-negative samples. Transcripts higher in HPV-positive samples are marked red. The horizontal line indicates adjusted p value of 0.05. Blue dots represent genes in the tumor inflammation signature. Analyses of mRNA expression were performed in R 3.6.2. Heatmaps were plotted with the R package ComplexHeatmap40. (https://www.bioconductor.org/packages/release/bioc/html/ComplexHeatmap.html).

Full size image

Differences between HPV-positive and negative HPSCC are related to the cell cycle

Comparison of the gene expression between HPV-positive and HPV-negative samples presented a  > tenfold overexpression of CDKNA (P16INK4A) in HPV-positive samples, confirming the p16 protein expression as evaluated by IHC. In addition, p53 and p21 showed a strong overexpression while CDK6, CDH11 and CCND1 were among the transcripts down-regulated in HPV positive samples (Fig. 1B, Table 2, Supplementary Table S1).

Table 2 Transcripts with > twofold and significant difference in expression between HPV positive and negative HPSCC.
Full size table

Immune genes linked to survival in HPV-negative HPSCC

When HPV-negative samples were evaluated separately, they clustered into the same three groups described above. Although surviving patients were clearly enriched in cluster 1, and underrepresented in cluster 3, this difference was not significant (p = 0.07, Fig. 2A). After adjusting for multiple testing, no transcripts remained significantly different between the patients who responded to treatment (no evidence of disease, NED) and those who did not (dead of disease, DOD) (Fig. 2B, Table 3, Supplementary Table S1). Genes found to have large differences in expression between the NED and the DOD group, as well as those with significant unadjusted p values were used to construct heatmaps shown in Supplementary Fig. S1. Genes with higher expression in the NED group included several calgranulins (S100A8, A9 and A12) showing similar expression patterns, and HLA class I and II genes (HLA-A, HLA-DQA1 and HLA-DQB1), while those with lower expression in the NED group included SPP1, PC and FANCA.

Figure 2
figure 2

(A) Heatmap of all HPV-negative samples. Genes are rows and samples columns. (B) Volcano plot of transcripts differing between surviving and non-surviving patients with HPV-negative HPSCC. Transcripts higher in survivors are marked red. Unadjusted p values are presented and the horizontal line indicates an unadjusted p value of 0.05. Blue dots represent genes in the tumor inflammation signature. Analyses of mRNA expression were performed in R 3.6.2. Heatmaps were plotted with the R package ComplexHeatmap40. (https://www.bioconductor.org/packages/release/bioc/html/ComplexHeatmap.html).

Full size image
Table 3 Transcripts with > twofold and/or significant† differences in expression between the NED and DOD group in HPV negative HPSCC.
Full size table

When comparing clusters 1 and 3 against each other, a large number of genes were found to be differentially expressed (Supplementary Fig. S2, Supplementary Tables S1 and S2). The majority of transcripts up-regulated in cluster 1, consisting of patients with good prognosis, were immune-related and included IDO1, NKG7, CD27, CCL5, CXCL9 and CXCL10.

Immune signaling pathways are enriched in HPV-negative patients with favorable prognosis

Gene set enrichment analysis (GSEA) was performed to identify cell signaling pathways correlating with clinical parameters (Supplementary Fig. S3, Supplementary Table S3). No pathways were significantly enriched in HPV-positive samples after adjusting for multiple testing. However, without adjusting the p values, cell cycle checkpoints were enriched in HPV-positive samples, while cytokine and extracellular matrix signaling were comparatively enriched in HPV-negative samples. When comparing the NED and DOD groups, as well as cluster 1 vs. cluster 3, immune and interferon signaling were significantly enriched for the NED group and cluster 1, with both groups containing tumors responding to treatment. Tyrosine kinase signaling was enriched among non-responders.

Gene expression in relation to tumor inflammation signature for HPV negative HPSCC

An 18-gene signature (CCL5, CD27, CD274, CD276, CD8A, CMKLR1, CXCL9, CXCR6, HLA-DQA1, HLA-DRB1, HLA-E, IDO1, LAG3, NKG7, PDCD1LG2, PSMB10, STAT1, TIGIT) has recently been proposed to be predictive for a clinically favorable response to treatment with PD-L1 inhibitors9,10,11. These genes were all overexpressed in tumors responding to treatment (Fig. 1B), and were markedly enriched in the cluster of HPV-negative patients with favorable prognosis (Supplementary Fig. S2). Confirming these observations, GSEA showed this gene signature to be significantly enriched in tumors responding to treatment (Supplementary Fig. S3, Supplementary Table S3). Clustering all HPV-negative samples based on their expression of these genes showed that samples with high expression of all or nearly all of these genes formed a cluster where all but one patient survived (Fig. 3A). Notably, this cluster of surviving patients consisted of the same five samples as those forming Cluster 1 in Fig. 2A, where 4/5 patients survived. This shows that the genes in the tumor inflammation signature drive clustering even when hundreds of genes are involved. To further illustrate the impact of this gene signature on the survival of patients with HPV-negative HPSCC, Fig. 3B shows a pathway enrichment plot for the 18-gene signature, which was significantly enriched in surviving patients.

Figure 3
figure 3

(A) Heatmap of all HPV-negative samples for the 18 genes in the tumor inflammation signature. The black box marks samples from cluster 1 in (B) Gene set enrichment plot for the 18-gene signature in relation to survival. Genes in the signature are marked on the x-axis and ranked from left to right based on their log2-fold change between tumors responding to treatment ant non-responding. Analyses of mRNA expression were performed in R 3.6.2. Heatmaps were plotted with the R package ComplexHeatmap40. (https://www.bioconductor.org/packages/release/bioc/html/ComplexHeatmap.html).

Full size image

Expression of S100A12 and HLA class I is linked to survival

Several S100 genes; S100A8, S100A9 and S100A12, showed an increased expression in HPV-negative tumors responding to treatment (Fig. 2B). Even though none of these genes was significantly overexpressed by itself, the fact that three related proteins all showed highly increased expression warranted further study. To validate this finding in a larger dataset, S100A12, showing the largest difference between responders and non-responders (3.5 times higher in responders) was chosen for further analysis by IHC (Supplementary Fig. S4). Seventy-four tumors from patients who had completed curative treatment (Table 1) were consequently analyzed for their expression of S100A12. A minority of tumors (10/74) showed complete absence or only weak expression in a minor part (< 33%) of the tumor. When progression-free, disease-specific and overall survival (PFS, DSS and OS) were analyzed in relation to the fraction of S100A12-positive cells, DSS was found to be significantly worse in patients with tumors having a low expression of S100A12 (p = 0.018) (Fig. 4A–C). For PFS and OS, low expression was also associated with low survival, although not significantly so (p = 0.098, p = 0.089, resp.).

Figure 4
figure 4

Kaplan–Meier curves presenting progression free (A,D,G), disease free (B,E,H) and overall survival (C,E,I) for patients with HPV-negative HPSCC in relation to S100A12 or HLA class I expression. (A–C) S100A12 expression dichotomized between fraction of positive cells more (red) or less (blue) than 33%. (D-F) HLA class I expression dichotomized between absent/low (blue) and medium/high (red). (G–I) Combined expression of S100A12 and HLA class I. Both high (green), one high, one low (red) or both low (blue). Notches denotes censored.

Full size image

Since immune signaling pathways were highly enriched, and several HLA antigens (HLA-A, HLA-DQA1 and HLA-DQB1) were overexpressed in patients with favorable prognosis, HLA class I A,B,C expression was also validated by IHC. PFS, DSS and OS were all significantly better for patients with medium or high HLA class I expression (p = 0.001, p < 0.001 and p = 0.001 for PFS, DSS and OS, resp.) (Fig. 4D–F).

Tissue expression levels of S100A12 and HLA class I were also analyzed together and found to be significantly correlated with PFS, DSS and OS (p = 0.005, p < 0.001 and p = 0.006, resp.) (Fig. 4G–I). Notably, none of the 7 patients with tumors lacking expression of both of these proteins survived for 3 years, and all had a PFS of at most 1.5 years.

FANCA, which was overexpressed in non-responders (Table 3) was also analyzed by IHC. However, no significant association between FANCA expression and survival was detected (data not shown).

Hazards analysis of S100A12 and HLA class I together with clinical parameters in relation to survival

S100A12 and HLA class I expression in relation to PFS, DSS and OS was also evaluated by univariate and multivariate Cox proportional hazards analyses (Table 4). Upon univariate analysis, hazards ratios for having a tumor with low S100A12 expression were 1.82, 2.51 and 1.85 for PFS, DSS and OS respectively, while the corresponding figures for low HLA class I expression were 2.31, 3.32 and 2.34. In the multivariate analysis, S100A12 and HLA class I were analyzed with sex, age, tumor size and nodal status as covariates. For PFS and OS, both HLA class I and age were significant in the multivariate analysis, while for DSS only HLA class I remained significant. When S100A12 and HLA class I expression were included as the only covariates in a multivariate analysis, HLA class I remained statistically significant in relation to PFS, DSS and OS, while S100A12 did not, demonstrating that the expression of these proteins is not quite independent (data not shown).

Table 4 Univariate and multivariate Cox regression analysis for 3-year survival in patients with HPV negative HPSCC.
Full size table

Combined analysis of HLA class I expression and CD8-positive tumor infiltrating lymphocytes predicts survival

Sixty-eight of the 73 tumors analyzed for HLA class I expression were in an earlier study evaluated for numbers of CD8-positive tumor infiltrating lymphocytes (CD8 + TILs)14. When data on CD8 + TILs was combined with HLA class I expression from the present study there was a very strong and significant difference in survival between patients having HPSCC with low HLA class I expression and low number of TILs and those with high HLA class I expression and high number of CD8 + TILs, PFS, DSS and OS (p = 0.004, p < 0.001 and p = 0.005, for PFS, DSS and OS resp., Supplementary Fig. S5).

Discussion

In the present study, expression of cancer and immune related genes in HPSCC biopsies was analyzed in relation to treatment outcome. A specific aim was to identify prognostic markers to help select HPSCC patients who, because of their poor prognosis, would benefit from surgery despite the side-effects. Gene expression analysis combined with enrichment analysis demonstrated that immune activation is linked to survival in HPV-negative HPSCC. The analysis indicated high expression of several HLA class I and II antigens as well as several calgranulins, S100A8, A9 and A12 as related to better treatment response. Expression of HLA class I (A,B,C) and S100A12 was further validated by IHC as markers of improved clinical outcome.

When performing gene set enrichment analysis (GSEA), immune signaling was significantly enriched in tumors responding to treatment while receptor tyrosine kinase (RTK) signaling was enriched in non-responding tumors. Tumor immunity and CD8 + lymphocyte infiltration are positive prognostic factors in many cancers including HPSCC, while RTK signaling is a negative prognostic factor14,15,16,17,18,19. The immune signatures here were characterized by the antigen presentation system through overexpression of interferons, HLA class I and CD8.

An 18-gene signature for inflammation response, suggested to be predictive for a favorable response to treatment with PD-L1 inhibitors9,10,11, and which was enriched in responding patients upon GSEA, was also studied separately in relation to survival. Clustering using these 18 genes identified 5/20 HPV-negative patients as potential responders to PD-L1 treatment. Notably, 4/5 of these patients survived for at least 3 years, indicating that these tumors respond better to RT/CRT treatment as well.

The correlation between high expression of HLA class I (A, B, C) and survival in HPSCC was not unexpected, as a number of studies have demonstrated similar associations for HPV-negative HNSCC from other subsites20,21,22,23,24,25. Despite the fact that HPSCC in general has a very poor prognosis, patients with tumors showing high HLA class I expression had a much improved treatment outcome, indicating the importance of tumor recognition by the immune system for these tumors. Further highlighting the importance of anti-tumor immunity, CD8 + TIL numbers have been linked to increased patient survival in several studies of HPSCC14,17,26. Here, the combined evaluation of HLA class I expression and CD8 + TILs demonstrated a major difference in survival between patients with tumors having high HLA class I expression and high numbers of CD8 + TILs vs tumors showing low values for both. This contradicts a previous HPSCC study, where no association was found between HLA class I expression and survival, although HLA class I was found to correlate with CD8. However, this study was limited to late-stage patients receiving neoadjuvant chemotherapy, which could explain the discrepancies in results17.

Notably, the HLA class II antigens HLA-DQA1 and HLA-DQB1 were also overexpressed in responders on the mRNA level. A positive correlation between high expression of HLA class II antigens and favorable clinical outcome has earlier been shown by IHC for HPV-negative oropharyngeal squamous cell carcinoma (OPSCC)22.

Low expression of S100A8, A9 and A12 (calgranulin A, B and C) mRNA was related to worse prognosis. Validation of S100A12, showing the largest expression difference between responders and nonresponders, showed tumors with absent or very low expression to have significantly worse clinical outcome. S100A8, A9 and A12 belong to a subgroup of S100 proteins, calgranulins, involved in senescence, aging and inflammation response. They are particularly expressed by cells in the innate immune system, e.g. neutrophils and macrophages27,28,29,30. S100 proteins have also been associated with epithelial differentiation in a single-cell study of head and neck cancer31. While S100A8 and A9 can act together by forming heterodimers, calprotectin, S100A12 functions as a homodimer. S100A8/S100A9 is often upregulated in cancer and is involved in inducing a pro-inflammatory response32. High S100A8 expression has been associated with increased disease specific and overall survival in both breast cancer and HNSCC33,34. In OPSCC and oral cancer, S100A8/S100A9 is often downregulated, associated with lower tumor grading, upregulation of apoptosis-related genes, decreased EGFR expression, epithelial differentiation and decreased migration and invasion35,36. While S100A12 is constitutively expressed by neutrophils and inducible in macrophages, expression in epithelial cells can lead to growth arrest27,37. S100A12 has received less attention than S100A8 and S100A9, but all three proteins were the focus of a study on OPSCC by Funk et al., where high expression of S100A8 and/or S100A12 were favorable prognostic factors38. A difference between the OPSCC study by Funk et al. and the present study was that much fewer HPSCC had low S100A12 expression. This may be similar to the lower levels of S100A8/S100A9 expression in OPSCC and OSCC as compared to other head and neck sites35.

There are several limitations caused by the small sample size. Since the number of samples was low and a large set of transcripts was analyzed, no transcript was significantly associated with survival after adjusting for multiple testing. In addition, the sample size led to patients receiving RT and CRT being analyzed together. It is thus possible that genes associated with response to only one of these treatments remained undetected. Another possible confounder is the higher age and lower frequency of CRT in the non-response group, while in the response group the tumors were on average of a higher stage. Thus, there is a need to further validate and extend the results presented in this study in a larger dataset.

In conclusion, we have identified an immune signature as a predictor of survival in HPV-negative HPSCC and validated S100A12 and HLA class I, particularly when combined with CD8, as markers for treatment response. This opens up prospects for applying immunotherapy and checkpoint inhibition in HPSCC, as well as using absent S100A12 and HLA class I expression as a means to indicate patients for whom the survival benefits of surgery probably outweigh the side effects.

Methods

Patients and tumor biopsies

Hypopharyngeal cancer biopsies with the ICD-10 codes C12.9 (pyriform sinus), C13.0 (postcricoid region), C13.1 (aryepiglottic fold, hypopharyngeal aspect), C13.2 (posterior wall of hypopharynx), C13.8 (overlapping sites of hypopharynx) and C13.9 (hypopharynx, unspecified location) were derived from two previously analyzed and published patient cohorts, covering in total 191 patients, treated 2000–2013 at the Karolinska University Hospital in Stockholm3,4. Pretreatment formalin-fixed paraffin-embedded (FFPE) biopsies from these cohorts were earlier examined for the presence of HPV DNA from 24 or 27 different HPV types by a Luminex bead-based assay and for the expression of p16INK4A (p16) in relation to clinical outcome3,4. Only samples both positive for HPV DNA and overexpressing p16 were considered as HPV positive. All three HPV positive samples were positive for HPV type 16. HPV DNA positive/p16 negative samples were not included in the study. Patient and tumor characteristics are presented in Table 1.

The study was performed according to permission 2009/1278-31/4 from the Regional Ethics Committee, Karolinska Institutet. Informed consent was obtained from all patients. All methods were carried out in accordance with relevant guidelines and regulations. All experimental protocols were approved by the Dept. Oncology-Pathology, Karolinska Institutet.

Selection of samples for gene expression analysis and immunohistochemistry

Inclusion criteria for selection of biopsies for gene expression analysis were complete patient records, small tumor size (T2 or T1), completed curative treatment and where patients either had a 3-year relapse free survival (group NED, no evidence of disease) or had died within 3 years due to the disease (group DOD, dead of disease). In addition, for HPV negative tumors, only patients who did not undergo local surgery were included. For immunohistochemistry, the criteria for inclusion were complete patient records, completed curative treatment and no local surgery.

Ten of the HPV-negative patients had 3 years disease-free survival, while the remaining 10 all died from the disease within 3 years. In the latter group, eight patients relapsed before death due to their cancer, while two had residual disease. All three patients with HPV positive tumors, defined as being both HPV DNA and p16 positive, survived 3 years and were relapse free.

Sample preparation and nanostring assay

Six 5 μm cuts/FFPE HPSCC biopsies were obtained using a microtome and mounted on PENMembrane 2.0 μm membrane slides (No.11505158) (Leica Microsystems, Wetzlar, Germany). Before microdissection cuts were deparaffinized for 2 × 15 min in fresh xylene, then rehydrated in decreasing percentages of ethanol (5 min each in 100%, 95%, 70%, 50%, 0% ethanol) and stained for 30 s with hematoxylin. Areas with tumors from the six cuts/tumor were then laser microdissected using a Leica LMD 7000 microscope (Leica Microsystems, Wetzlar, Germany) and the Laser Microdissection System (version 7.6.5684), pooled in microcentrifuge tubes containing PKD buffer from the RNeasy FFPE kit (Qiagen, Halden, Germany) and kept on ice.

RNA was extracted directly after microdissection using the RNeasy FFPE kit according to manufacturer's instructions. RNA concentrations were evaluated using the the Qubit RNA HS Assay Kit with the Qubit 4 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). All included samples had RNA concentrations > 7 ng/ul and all except one > 18 ng/ul. No amplification step was included before analysis on the nanostring panel.

Gene expression was analyzed at the KIgene facility, Center for Molecular Medicine, Karolinska Institute, using the nCounter IO360 panel, (NanoString Technologies, Seattle, WA, USA) according to manufacturer's instructions, using the nCounter Sample Prep Station with FLEX configuration (NanoString Technologies) and the nCounter Digital Analyzer 5 s (NanoString Technologies) for reading the samples.

Immunohistochemistry and evaluation of HLA class I, S100A12 and FANCA

Expression of HLA class I was analyzed using the mouse monoclonal HLA class I ABC antibody (ab70328, Abcam, Cambridge, UK), recognizing the heavy chains of human HLA class 1 A, B and C, diluted 1:1000. Staining was performed essentially as described previously with Vectastain Elite ABC HRP kit, PK-6100, and DAB Peroxidase (HRP) Substrate Kit SK-4100, both from Vector Laboratories, Burlingame, CA, USA24. As secondary antibody, biotinylated anti-mouse IgG (BA-2000, Vector Laboratories), diluted 1:200, was utilized.

Evaluation of staining was performed by three investigators (AN, TR and RU), blind to clinical outcome. In case of disagreement, a consensus decision was made. Fractions of HLA class I expressing cells were evaluated semi-quantitatively in five grades: 0 (0%), 1 (1–25%), 2 (26–50%), 3 (51–75%), and 4 (76–100%). Intensity of expression was scored on a four-tier scale as absent, weak, medium and strong staining. For analysis of survival, intensity of expression was dichotomized into absent/weak or medium/strong. Examples of staining are presented in Supplementary Fig. S4.

S100A12 expression was analyzed using the Calgranulin C (19F5): sc-101347 mouse monoclonal antibody from Santa Cruz Biotechnology (Dallas, TX, USA), diluted 1:60. Staining was performed as for HLA class I above and evaluated blind by three investigators (AN, TR and RU). Intensity of expression and fraction of positive cells was evaluated in a similar manner as Funk et al.38. The fraction of S100A12 positive cells was evaluated as 1 = 0%, 2 = 1–33%, 3 = 34–67% and 4 = 68–100%. Intensity was evaluated as 1 = negative, 2 = weak, 3 = medium and 4 = strong. However, fraction and intensity were not combined into a score as described by Funk et al.38. For evaluation of survival, the fraction of positive cells was dichotomized into 1 + 2 (0–33%) or 3 + 4 (34–100%). Examples of staining are presented in Supplementary Figur. S4.

Expression of FANCA was analyzed with the rabbit anti-FANCA polyclonal antibody ASJ-C2020F from Nordic Biosite (Stockholm, Sweden) at a dilution of 1:100. As secondary antibody biotinylated Goat Anti-Rabbit IgG Antibody (BA-1000), diluted 1:200, was utilized. The intensity and frequency of FANCA was evaluated as number of 10-percentiles of cells with absent, weak, intermediary or strong staining. For evaluation of FANCA expression in relation to survival, fractions of cells with intermediary and strong staining were combined.

Evaluation of patient survival in relation to IHC analysis

For analysis of survival in relation to expression of S100A12 or HLA class I expression by IHC, survival was measured in days from diagnosis until an event occurred, or until three years after diagnosis when patients were censored. Survival was evaluated as in the study by Landin et al.14. The Kaplan–Meier estimator was used for evaluation of three-year progression free survival (PFS), disease specific survival (DSS) and overall survival (OS). For PFS, disease recurrence, or death of any cause before three years were considered as events. For evaluation of disease specific survival (DSS), death with documented relapse was considered as an event while patients that had deceased before three years from other causes were censored at time of death. For OS, death of any cause before three years was considered as an event. Differences in the survival of patients in relation to S100A12 or HLA expression were evaluated using the log-rank test. Univariate and multivariate hazard ratios (HR) were calculated using Cox proportional hazards regression analysis. Calculations and analyses were performed using IBM SPSS Statistics, version 26.0.

Nanostring data analysis

Only endogenous mRNA transcripts were kept for downstream analysis, removing housekeeping genes. In order to increase power upon multiple testing, transcripts were filtered by mean and variance, keeping only those transcripts having above median variance and above median gene expression. This resulted in 332 out of 750 endogenous mRNA transcripts remaining (Supplementary Fig. S6). From the 18-gene tumor inflammation signature of response to PD-1 inhibition (9–13), seven genes that did not pass the mean and/or variance thresholds were nevertheless kept for downstream analysis. All raw nanostring counts are found in Supplementary Table S4. The expression data were log2-transformed and median-centered by sample before analysis, and re-centered whenever sub-group analyses not including all patients or transcripts were performed.

Student’s t-test was used for comparing continuous variables and Fisher’s exact test for comparing categorical variables. p values were adjusted for multiple testing using the Benjamini–Hochberg method. p values below 0.05 were considered significant.

For gene set enrichment analysis, enrichment scores and significance values were calculated as described in Subramanian et al.39. Highly similar pathways were collapsed into one in order to show only independently enriched pathways. All analyses of mRNA expression were performed in R 3.6.2. Heatmaps were plotted with the R package ComplexHeatmap40. For hierarchical clustering, Euclidean distance and complete linkage were used as parameters GSEA was performed using the fgsea R package41. Functional pathways were taken from Reactome42.

Data availability

The nanostring data generated and analysed during this study are included in Supplementary Table S4. Clinical data used for survival analysis based on evaluation by IHC are not publicly available, due to regulations on transfer of such data. However, the data utilized in the analysis in the present study is available from the corresponding author on reasonable request.