Immune cell clusters inside tumors linked to better survival and immunotherapy response in head and neck cancer

Derived from 12 chemokine gene signatures [13, 14] and 9 marker genes [6], we proposed a novel gene signature composed of 13 chemokine genes (LAMP3, CCL2, CCL3, CCL4, CCL5, CCL18, CCL19, CCL21, CXCL9, CXCL10, CXCL11, CXCK13, CXCR4) to specifically evaluate HNSCC TLS at the transcriptomic level from the dataset of TCGA (HNSCC) (Fig. 1A). We excluded BCL6 as its expression was not in accordance with that of the other genes in HNSCC samples (Supplementary Fig. 1A, B). Subsequently, the TLS score was calculated and divided into high and low for each HNSCC tumor specimen in TCGA. As shown in Fig. 1A, this stratification method was powerful enough to distinguish the signature gene expression and provide an easier way to determine the correlation between clinical stage and pathological grade among TLS subtypes. We then probed into the predicative role of TLS signature in HNSCC prognosis with the finding that TLS-hi tumors were associated with better overall survival (OS) rates (p < 0.05) than TLS-inter and TLS-low tumors (Fig. 1 B, C). However, there was no prognostic significance in disease-free survival (DFS), recurrence-free survival (RFS) or metastasis-free survival (MFS) rates among different TLS subtypes (Supplementary Fig. 1C–E). The positive correlation of TLS status with HNSCC prognosis was further verified in other two datasets and the merge data of TCGA-HNSCC with CPTAC-HNSCC (Supplementary Fig. 1F–H). To further confirm the predictive role of TLS status in overall survival for HNSCC, we integrated TLSs with other relevant risk factors, such as stage, oropharynx, HPV status, gender and age. In univariate analysis, TLSs (high), stage (I/II) and gender (female) were found to be positive associated with favorable survival for HNSCC. These discoveries were further confirmed by multivariate Cox regression model. HPV status failed to prove its predictive role in better OS for HNSCC in both univariate and multivariate analysis (Supplementary Fig. 2A, B).

Subsequently, we investigated the contribution ratio of every chemokine in TLSs to the prediction of OS rates. Interestingly, no single chemokine showed a significant correlation with favorable OS except for CCL5 (p = 0.048) (Fig. 1D), indicating a synergistic effect of the TLS signature instead of individual chemokines on prognostic prediction. To explore the expression of 13 chemokines in TME, HNSCC single-cell RNA-Seq data [15] were used to identify the major cell types contributing to the TLSs signature, which was composed of CD8/CD4 T cells, NK, B cells, monocyte, macrophage, dendritic, fibroblasts and endothelial (Supplementary Fig. 3A, B). Due to the heterogeneity of TLSs among different cancer types, we tried to identify whether this heterogeneity also exists between normal and tumor tissues, recurrent and primary tumors. Compared to paracarcinoma tissues, carcinoma tissues had higher expression levels of CXCL9, CXCL10, CXCL11, CXCL13, CCL3, CCL4, CCL5 and LAMP3, resulting in higher TLS scores (p < 0.001) (Fig. 1E, Supplementary Fig. 2D). A similar phenomenon was observed in the secondary versus primary groups, characterizing by higher level of chemokines except for CCL2, CCL19 and CCL21, even though with no statistical difference in TLSs score (Fig. 1F, Supplementary Fig. 3D), indicating that secondary tumors may evolve from primary tumors in a similar way to how carcinogenesis arises from normal tissues. However, there was no difference between recurrent and primary tumors in TLS signature expression and score (Supplementary Fig. 3C, D).

The presence of tumor-infiltrating TLSs may be associated with increased immunocyte invasion and immunological microenvironment. The majority of immune checkpoints and T cell receptors, including PD-L1, LAG3, TIGIT, HAVCR2, PD-1, CD80, CTLA4 and IDO1, showed higher expression in TLS-hi tumors (Fig. 2A). The expression relevance of these checkpoints was more significant in TLSs than in T cells, CD8 T cells, cytotoxic T, B lineages, NK, monocyte and mDC (Fig. 2B). TLSs were also associated with higher immune score and infiltration of cytotoxic T, monocytic lineage and T cells, but not stromal parameters, which directly indicates good responses to immunotherapy (Fig. 2C). Subsequently, we probed the tumor infiltrative immune landscape between the TLS-hi and TLS-low groups. A higher enrichment of memory, immature and activated B cells, but not naïve and plasma B cells, and increased proportions of memory and activated CD4 + T and CD8 + T cells (p < 0.001), and activated dendritic cells, type I T helper cells, MDSC, were found in TLS-hi tumors (Fig. 2D, E, Supplementary Fig. 4A, B). To explore whether tumor subsite poses a strong bias in the immune features between TLS-high and TLS-low phenotypes, we categorized HNSCC into three subtypes, including hypopharynx, oral cavity and oropharynx. Cancers from these tumor subsites showed similar immune cells infiltration and transcript level of immune checkpoints (Supplementary Fig. 5), and from which we can exclude the tumor subsite bias.

In addition, the TLS-hi subtype had a more powerful immunostimulatory profile, high expression of chemokines, a higher enrichment in human leukocyte antigens (MHC) and TCRs, and a higher score of tumor-infiltrating immune cells (TIICs) (CD8 T, dendritic cell, macrophage, NK, Th1) effectors (p < 0.0001), demonstrating a better function of these immune cells (Fig. 2F, G).

The relatively higher infiltration level of exhausted CD8 + T cells in TSL-hi tumors was also found in current analysis (p < 0.0001) (Fig. 3A). To further validate this discovery, CD8 + T cells are further classified into three groups: (1) Ly108⁺ progenitor T cells (Tpro), with high capacity to produce IFN-γ and TNF-α, (2) CX ₃ CR1⁻ Ly108⁻ exhausted T cells (Texh) with the lowest potential to produce inflammatory cytokines, represented as truly exhausted phenotype; (3) CX ₃ CR1⁺ cytolytic T cells (Teff) with capacity to expression Granzyme B and cytolytic function, during persistent infection or tumor progression [16]. We demonstrated that all these three subtypes showed obvious elevation in TLSs-hi groups (p < 0.0001) (Fig. 3B–D). BATF, a key transcription factor during progenitor to cytolytic effector CD8 + T cells, is specifically expressed in Tregs and exhausted CD8 + T cells among HNSCC scRNA-Seq data (Fig. 3E). As expected, BATF also highly enriched in TLSs-hi samples as compared to TLSs-low samples (p < 0.0001) (Fig. 3F). Exhausted CD8 + T cells hold a negative correlation with immune response and overall survival. We next probed into whether TLSs statue had a decisive role to the clinical significance of exhausted CD8 + T cells. In TLSs-hi group, patients with high infiltration of exhausted CD8 + T cells showed worse outcome (p < 0.01) (Fig. 3G). When it comes to TLSs-low group, there is no difference in OS rate between high and low infiltration groups (Fig. 3H). Immune evasion has an intimate correlation with T cell dysfunction and exclusion. An interesting finding was that TLSs-hi group hold more dysfunction T cells, TIDE score and higher concentration of IFN-γ, whereas TLSs-low group had higher exclusion score, indicating that TLSs-hi and TLSs-low patients experience immune scape through different mechanisms (Fig. 3I).

Subsequently, we investigated the interaction between TLS subtypes and the universal landscape of the TME by focusing on immune cells, immune scores, stromal cells and stromal scores. We discovered that TLSs mainly participated in immune reactions (monocytic lineage, immune score, CD8 T cells and B lineage) but showed little correlation with stromal cells and scores, even though stromal cells also expressed these chemokines of TLSs signature (Fig. 4A). We then integrated TLSs, immune-stromal parameters, tumor stage, gender and age as a whole panel for survival assessment. In the univariate analysis, TLSs, B lineage cells, T cells, tumor stage and gender were significantly different, and only T cells low (p < 0.01) and tumor stage (III/IV) (p < 0.01) were retained as prognostic predictive markers for HNSCC in the multivariate analysis (Fig. 4B, C). As expected, the TLS subtype played a significant role in the clinical outcome of advanced (stage III/IV) HNSCC (p < 0.05) but not in that of early-stage HNSCC (Fig. 4D, E). Interestingly, we also found that TLS-hi represented a favorable prognosis in male patients (p < 0.05). (Fig. 4F, G).

We then tried to identify the determined role of TLSs on the clinical significance of immunocytes. Higher infiltration levels of T cells and B lineage cells in the TME signified better overall survival (p < 0.01) (Fig. 4H and J). However, cytotoxic T cells, neutrophils and mDCs were not predictive for OS outcomes (data not shown). Subsequently, we investigated whether the presence of TLSs has a measurable influence on the clinical significance of T cells and B lineage cells. As shown in Fig. 4I and K, TLS abundance was not associated with the level of B lineage cells both in the high and low groups. Regarding T cells, TLS-hi tumors showed relatively better OS outcomes than TLS-low tumors in the T cell high group but not in the T cell low group (Supplementary Fig. 6A, B). These results indicate that TLSs have a determined role in the clinical significance of T cell levels. To confirm the predictive role of T cell and B lineage in overall survival of HNSCC, we integrated T cell and B lineage with stage, HPV status, oropharynx, gender and age, with the finding that both HPV status and oropharynx have no determined role in the clinical significance of T cell and B lineage (Supplementary Fig. 7A–D).

To unravel the molecular mechanisms and functional changes between TLS-hi and TLS-low groups, differentially expressed genes (DEGs) were identified from the HNSCC TCGA gene set (Fig. 5A). Subsequently, intrinsic pathway alterations were analyzed based on these DEGs, characterized by significant enrichment of cytokine–cytokine receptor interaction, chemokine/TCR/BCR signaling, Th1/2/17 differentiation, NK cell-mediated cytotoxicity, and antigen processing and presentation (Fig. 5B). The bubble chart further confirmed these discoveries and emphasized enrichment in immune response, T cell activation, and lymphocyte differentiation and proliferation (Fig. 5C). We then tried to explore whether TLS subtype associated DEGs have correlation with stromal and immune scores. In HNSCC samples, stromal and immune scores have been categorized as stromal-hi/low and immune-hi/low according to the DEGs (Supplementary Fig. 8 A–C). Then the up- and down-expressed genes were presented in TLSs hi/low, stromal-hi/low and immune-hi/low, respectively (Supplementary Fig. 8 D–F). Finally, DEGs of different subtype of stromal and immune score were merged with that of TLS subtype, with the finding that TLS had more intimate association with immune score than stromal score (Supplementary Fig. 8G).

Subsequently, we investigated the level at which TLSs interact with other genes or signaling pathways in HNSCC tumors to affect the overall survival. We finally filtered 102 genes, among which 63 genes had positive coefficient and the remaining 39 genes hold a negative coefficient. Pathway enrichment showed that necroptosis, mitophagy, p53 and Wnt signaling, and the NOD-like receptor signaling pathway, complement activation, tumor necrosis factor-mediated signaling pathway, which served as anticancer effectors, enhanced the association between TLSs and better clinical outcomes, which showed the synergistic effects of TLSs and the enriched signaling pathways in contributing to a more favorable prognosis (Fig5D, E). When the selective threshold p value was adjusted to ≤ 0.001, the metabolic and AMPK signaling pathways were filtered and determined to be associated with TLSs (Supplementary Fig. 9A). Inspired by these enrichment results, we probed into the metabolic landscape in tumors with different TLS subtypes. As shown in Fig. 5F, TLS-hi tumors were featured by enhanced negative regulation of oxidative phosphorylation, glutathione metabolism, carbon metabolism, glycoprotein biosynthesis and positive regulation of lipid, phospholipid and phosphorus metabolic processes, indicating a metabolic shift from carbohydrate to lipid metabolism in TLS-hi tumors in comparison with TLS-low tumors.

High tumor mutational burden (TMB) was significantly associated with a higher survival rate among HNSCC patients (p < 0.05) (Supplementary Fig. 9B). We then probed into the influence of TLS subtype on the TMB in HNSCC. Unfortunately, the TLS subtype has no determined role in the clinical significance of TMB (data not shown). Subsequently, we explored the driver genes and their mutation diversity between the TLS-hi and TLS-low groups. The commonly mutated genes with high frequency across the two groups are listed in Fig. 5G, H, which showed a higher mutation frequency of MUC16 (22% vs. 17%), DNAH5 (21% vs. 14%), PKHD1L1 (19% vs. 12%), NSD1 (15% vs. 8%), KYR2 (15% vs. 9%), and COL11A1 (12% vs. 6%) in TLS-hi versus TLS-low tumors. However, CASP8 (22% vs. 10%) was found to have a strikingly higher mutation frequency in the TLS-low group. Regarding the mutation types, no obvious difference was found in the number of deletions, in frame and splice mutations (Supplementary Fig. 9C, D). A slight change in amplification mutation frequency and notable diversity in missense and truncating mutations were observed between TLS-hi and TLS-low tumors (Supplementary Fig. 9E). In addition, we also found that TLS-hi tumors hold a higher TCR and BCR richness (p < 0.0001) (Fig. 5I). These variations in the mutation landscape and TCR/BCR may contribute to the predictive effect of TLSs on survival outcome.

Innate lymphoid cells (ILCs) are tissue-resident lymphocytes composed of group 1 LCs (ILC1), including NK cells and non-NK cells, group 2 LCs (ILC2) and group 3 ILCs (ILC3), containing ILC3 cells and lymphoid tissue inducer cells (LTi), and participate in the TLS induction. In HNSCC samples, the signature of ILC1 and ILC3 was co-expressed at a high level with that of TLS-hi subtype, demonstrating that the presence of ILC1/3 cells may have a positive correlation with TLSs induction (Fig. 6A, B). More specifically, co-expression analysis of TLSs and ILC3 signatures indicated an enrichment of the LTα, LTβ and TNF families, including TNFRSF4, TNFSF13B and TNFSF14, which have an intimate association with TLS neogenesis (Fig. 6B). Further pathway enrichment confirmed the biological function of ILC3 cells featured by the activation of Th1/2/17 differentiation, cytokine–cytokine receptor interaction, positive regulation of immune effector process, programmed cell death and T cell activation (Fig. 6C, D). In regard to ILC1 cells, significant elevation of chemokine, Toll-like receptor signaling pathway, NK cell-mediated cytotoxicity and immune response were observed, implying that ILC1 cells and ILC3 cells exert distinct functions during TLS neogenesis and further immune response (Supplementary Fig. 10A, B).

Besides innate immunity, stromal cells, such as fibroblasts, have also been reported to participated in TLS neogenesis. To clarify which subtype of fibroblasts contributing to the TLS development, we collected the pivotal genes mediating the TLSs neogenesis, such as IL-6, TNF-α, LTA, LTB, CCL19, CCL21, CXCL12, CXCL13, TNFRSF1A/B and IL-17. Then we calculated the correlation of these marker genes with the fibroblasts from HNSCC single cell data and identified that cancer-associated fibroblasts (CAFs) and myofibroblast hold the positive association with TLSs neogenesis (Fig. 6E). Furthermore, inflammatory CAF (iCAF) and myofibroblast CAF (myCAF) showed the most relevant with this gene set as compared with antigen-presenting CAF (apCAF), demonstrating that these two subtypes of CAFs may have an important role in TLS neogenesis (Fig. 6F).

We next proposed a novel HNSCC TME classifications (HNSCC-TCs) based on the comprehensive estimation of the whole tumor landscape, including stromal cells and stromal score, TLSs, T cell dysfunction and immune evasion features, antitumor effective cells and malignant tumor cells (Fig. 7A). Five HNSCC-TCs presented the high heterogeneous TME profiles. HNSCC-TC 1, the “immune desert” subtype, scarce infiltration of immune-promoting and suppressive cells, low level of TMB and TLSs, and high proportion of malignant cells. HNSCC-TC 2, a few infiltrations of both immune cell subtypes, high level of TMB and low level of TLSs, and high proportion of malignant cells. HNSCC-TC 3 features high infiltration level of both types of immune cells, low level of TMB and high level of TLSs, and a high density of stromal cells and high stromal score. HNSCC-TC 4 features scarce infiltration of immune cells, a high density of stromal cells and high stromal score, low level of TMB and TLSs. HNSCC-TC 5, high infiltration of both types of immune cells, a moderate density of stromal cells and stromal score, and high level of TMB and TLSs (Fig. 7B). Subsequently, we verified the feasibility and rationality of the current classification by association analysis with patient survival outcomes. Patients with HNSCC-TC 4 exhibited the worst OS outcomes in comparison with those with HNSCC-TC 1, indicating that a high density of stromal cells aggravated the immune desert environment and exerted synergistic effects on immune evasion and tumor progression. This finding was also proven by the OS rate difference between HNSCC-TC 3 and HNSCC-TC 5 patients, which was also characterized by increased infiltration of stromal cells and a higher stromal score. The prognosis of HNSCC-TC 2 was comparable with that of HNSCC-TC 3, signifying that the higher infiltration of immunologically promoted cells could be counteracted by the abundance of stromal cells and other immunosuppressive factors. Patients with HNSCC-TC 5 showed the best OS prognosis as compared to other groups, featuring with high infiltration of immune cells and high score of TLSs, low score of stromal (Fig. 7C). In addition, a positive correlation of enriched stromal elements with low TMB, such as in HNSCC-TCs 3 and 4, indicated that stromal factor is also an influencing parameter for determining immunotherapy responses and patient prognoses (p < 0.05) (Fig. 7D).

We then explored the TLS distribution in HNSCC-TC 1–5 and found the highest level of TLSs in HNSCC-TC 5, followed by HNSCC-TC 3, 2, 4, and 1, showing identical trends with PD-L1 expression but not TMB distribution, which further confirmed the co-expression of TLSs and immune checkpoints (Fig. 7F). As for immune evasion, HNSCC-TC 3 presented higher T cell dysfunction score than HNSCC-TC 5, and HNSCC-TC 4 was estimated with higher exclusion score, signifying that stromal score play a determined role in T cell dysfunction and exclusion, thus promoting the immune scape (Supplementary Fig. 10C). This TME classification in HNSCC TCGA was re-verified in another database and showed good repeatability in cluster distinguish and survival prediction (Supplementary Fig. 11A). The difference of overall survival among different clusters were similar to that of TCGA-HNSCC (Supplementary Fig. 11B). Furthermore, we discussed the dynamic change of TLS level during ICB treatment with the findings that TLS level showed no significant difference between CR/PR and SD groups, but revealed higher TLS enrichment in CR/PR patients than PD patients (Supplementary Fig. 11C).

We then probed into the underlying mechanism contributing to the difference in TME composition between HNSCC-TC 5 (the subtype with the most abundant immune-promoting factor and the best prognosis) and HNSCC-TC 1 and 4 (the subtype with a scarcity of immune-promoting factors and the worst prognosis). A significant reduction in cytokine-cytokine interaction, chemokine signaling, antigen processing and presentation, and T cell differentiation, was observed in the HNSCC-TC 1 and 4 groups, which indicated a lower immune response activity and stronger stromal biological process (Fig. 7E). We also performed further in-depth analyses to identify markers that were specifically expressed in tumor cells and found that 9 genes (COL4A6, LANCL2, PLA2R1, PLAC1, IL1A, MAGEA12, RNF128, GLI2 and DMP1) participated in the inhibition of TLS neogenesis, which deserved further molecular investigation.

Finally, we probed whether TLS subtypes are intimately associated with different therapeutic regimens. Unsurprisingly, TLS abundance was not associated with the regulation of radiotherapy, chemotherapy or EGFR targeted therapy processes (Supplementary Fig. 12). We then investigated whether the established TLS signature plays a role in OS outcome and immunotherapy response prediction. In HNSCC- immunotherapy cohort (GSE159067↗), a total of 102 patients were treated with PD-1/PD-L1inhibitors, and TLS status also showed its predicative role in overall survival of HNSCC with a p value of 0.028, and the TLS score was highest in the CR patients, followed by the PR, SD/PD patients (Fig. 8A–C). As expected, TLS level presented a gradually increased trend among desert, excluded and inflamed immune phenotypes (Fig. 8D).

Subsequently, we verified the signature using urothelial carcinoma data from GSEA, which containing 348 patients who received Atezolizuma, a PD-L1 antibody. In this ICB-treated cohort, patients were categorized into two groups according to the PD-L1 expression in tumor biopsies, IC1(1 < = PDL1 < 5%) and IC2/3 (PDL1 > 5%), with the finding that IC2/3 patients reached objective response rate (ORR) at 26% and CR rate at 11%, and a higher OS rate in CR/PR patients as compared with that of SD/PD patients. Data analysis presented that patients with high TLS levels hold significantly better OS rates than those with low TLS levels (p < 0.01) (Fig. 8E). The CR group had the highest abundance of TLSs than that in the PR (p < 0.05) and PD groups (p < 0.001), whereas there was no significant difference in TLSs level between the CR and SD groups (Fig. 8F). What’s more, we verified the expression of different gene sets represented various treatments (not performed in HNSCC-immunotherapy cohort due to the lack of transcript data of relative genes), such as immunotherapy, anti-EGFR therapy, anti-FGFR3 therapy and anti-PPARγ therapy, and found that TLSs was intimately correlated only with immunotherapy (Fig. 8G).

Publicly available HNSCC datasets were obtained from the TCGA (HNSCC), GEO (GSE41613↗, GSE42743↗ and E-MTAB-1328, GSE137564↗, GSE68858↗, GSE159067↗), CPTAC-HNSCC and IMvigor210 databases, which contain clinical characteristic, therapeutic regimen, corresponding response, follow-up, RNA-Seq and somatic mutation data. The HNSCC RNA sequencing data (FPKM values) were downloaded via the R package TCGA biolinks and were subsequently transformed into transcripts per kilobase million (TPM), which was identical to the microarray data downloaded from the Array Express database. The somatic mutation data were retrieved, and TMB was calculated by using VarScan2. The raw data from GEO databases were normalized based on the chip platform. A total 482 HNSCC samples in TCGA data were enrolled for overall survival analysis. Due to the data missing of follow-up in 121 patients, 361 patients were finally filtered for DFS, RFS and MFS estimation. TMB data was available in 477 patients. This project was approved by the ethics committee of Shanghai Ninth People’s Hospital Affiliated to Shanghai Jiao Tong University.

We retrieved the standard data after log2 scale and RMA normalization from the TCGA database. Then, the geometric mean of the genes (candidates for the TLS signature) in every sample was assessed by an MCP counter method to filter specific markers (CCL2, CCL3, CCL4, CCL5, CCL18, CCL19, CCL21, CXCL9, CXCL10, CXCL11, CXCL13, CXCR4, LAMP3) for HNSCC and to form the TLS score. Samples with scores greater than the third quartile were classified as TLS-hi, and those with scores lower than the third quartile were defined as TLS-low [9].

Differentially expressed genes between TLS-hi and TLS-low cohorts were determined by employing the R package limma and setting the cutoff of statistical significance at a p value < 0.05 and an absolute fold-change > 1.5. These up- or downregulated genes were then subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses using the R package clusterProfiler (v3.14.3) based on the criteria for a cutoff p value of < 0.05 and an adjusted p value of < 0.2. Those GO categories FDR < 0.05 were considered as significant enriched. While pathways with a p value < 0.05 were regarded as enriched. Only those GO categories or pathways contains≧5 DEGs were kept for further analysis.

We applied Cox-PH survival regression to analyze how the CTLs interact with the other genes in the tumor to affect survival outcomes. We used a linear hazard model (= a × CTL + b × V + d × CTL × V + c) with the Cox-PH regression model [32]. In the Cox-PH model, the death hazard was estimated through patient survival information. The variable V represents the expression level of a candidate gene in the test. Since we selected datasets where CTL levels correlated with favorable survival outcomes, the coefficient a was always negative. The association slope between CTLs and hazard was a + d × V. If the coefficient d was positive, a higher V level would flatten the slope between CTLs and hazard, indicating a reduced association between the cytotoxic T cell level and better survival outcomes. If d was negative, a higher V level would sharpen the slope between CTLs and hazard, indicating an increased association between the cytotoxic T cell level and better survival outcomes. Collectively, a positive coefficient “d” indicates a reduced association between the TLS level and better survival outcomes, whereas a negative coefficient “d” signifies an enhanced correlation between the TLS level and better OS outcomes.

The immunological features of the TME were composed of immunomodulators, chemokine-chemokine interactions, tumor-infiltrating lymphocytes (TILs), immune checkpoint expression and ICB resistance. First, a panel of 122 immunomodulators containing immune stimulators, MHCs, receptors and chemokines was collected and analyzed as previously reported [33]. Second, the infiltration status of distinct immune cells in the TME of HNSCC between the TLS-hi and TLS-low groups was quantified by using four independent algorithms, including the “CIBERSORT” R package, MCP counter, TIMER and XCELL, which were all analyzed based on the LM22 signature. The calculation results were then remodeled by scRNA sequencing data through LM22. C2 analysis[15]. Third, a signature of 20 immune checkpoints was retrieved to investigate the association between TLSs and the expression of immune repressive molecules [34]. Finally, “ESTIMATE” was further applied to calculate the immune and stromal scores for HNSCC individuals by calculating the single-sample gene-set enrichment analysis (ssGSEA).

Unsupervised clustering was applied to clarify different associated patterns of immune and stromal factors. Consensus parameters, such as cancer-associated fibroblasts (CAFs), other stromal cells, stromal score, tumor-associated macrophages (TAMs-M2/M1), ICB resistance, MDCSs, Tregs, exhausted CD8 + T cells, T cells/CD8 + T cells/cytotoxic T cells, TLS, immune score, IFNγ, total TMB and malignant cells, which represent the whole TME landscape, were extracted from the TCGA database and then integrated and clustered according to their intrinsic interactions. Gene signatures related to T cell dysfunction and tumor immune evasion used for clustering were identified from tumor cohorts through association analysis of each marker with TIL infiltration in influencing patient survival after ICB therapy [32].

The somatic mutation datasets of patients in the TLS-hi and TLS-low groups were extracted from the TCGA data portal (https://www.cancer.gov/tcga/↗). The total number of nonsynonymous mutations was selected as the representative to calculate the total mutational burden of the tumors. The mutation frequency and somatic alterations of common driver genes were listed in both groups. The driver genes of HNSC were identified using the “maftool” R package.

Hierarchical clustering analysis was performed on the R “hclust” function using the complete method to identify the number of clusters in TCGA-HNSC based on the expression pattern of signatures. Univariate and multivariate Cox proportional hazards regression models were used to assess the association between the signatures and overall survival with/without clinical variables. The hazard ratio (HR) and 95% confidence interval (CI) were calculated. One-tailed or two-tailed Wilcoxon rank-sum or student tests were used to compare two groups. The Kaplan–Meier method and log-rank test were conducted to compare survival differences between two tumor groups. All statistical analysis was performed using R/ Bioconductor (version 3.6.1).

None of the authors has a commercial or financial conflict of interest involved in current research.

The datasets generated during the current study are available from the corresponding author on reasonable request.