Purpose:

To identify immune subtypes and investigate the immune landscape of squamous cell carcinomas (SCC), which share common etiology and histologic features.

Experimental Design:

Based on the immune gene expression profiles of 1,368 patients with SCC in the Cancer Genome Atlas (TCGA), we used consensus clustering to identify robust clusters of patients and assessed their reproducibility in an independent pan-SCC cohort of 938 patients. We further applied graph structure learning-based dimensionality reduction to the immune profiles to visualize the distribution of individual patients.

Results:

We identified and independently validated six reproducible immune subtypes associated with distinct molecular characteristics and clinical outcomes. An immune-cold subtype had the least amount of lymphocyte infiltration and a high level of aneuploidy, and these patients had the worst prognosis. By contrast, an immune-hot subtype demonstrated the highest infiltration of CD8+ T cells, activated NK cells, and elevated IFNγ response. Accordingly, these patients had the best prognosis. A third subtype was dominated by M2-polarized macrophages with potent immune-suppressive factors such as TGFβ signaling and reactive stroma, and these patients had relatively inferior prognosis. Other subtypes showed more diverse immunologic features with intermediate prognoses. Finally, our analysis revealed a complex immune landscape consisting of both discrete clusters and continuous spectrum.

Conclusions:

This study provides a conceptual framework to understand the tumor immune microenvironment of SCCs. Future work is needed to evaluate its relevance in the design of combination treatment strategies and guiding optimal selection of patients for immunotherapy.

This article is featured in Highlights of This Issue, p. 3471

Translational Relevance

Immunotherapy is being used to treat an increasing number of patients with cancer including several types of squamous cell carcinomas (SCC). However, the response rates are relatively low and survival benefit is limited to a subset of patients. A better understanding of the tumor immune microenvironment is needed to improve response and outcomes for immunotherapy. In this multicohort retrospective study, we present the identification, independent validation, and comprehensive molecular characterization of six reproducible immune subtypes of SCCs. We found that each of the immune subtypes was associated with distinct gene expression profiles, and accordingly demonstrated widely different patterns in tumor genetic aberrations, tumor-infiltrating immune cell composition and functional orientation (immune stimulating vs. suppressive), cytokine profiles, and, importantly, clinical outcomes. This study provides a conceptual framework to understand the tumor immune microenvironment of SCCs, and may have clinical implications for the design of novel immunotherapies and rational combination strategies.

Immunotherapy is becoming a pillar of modern cancer treatment. In particular, immune checkpoint blockade (ICB) such as anti-PD1 antibodies have demonstrated durable response and unprecedented clinical benefit in a subset of patients across multiple types of solid tumors (1–6). However, the response rates for single-agent ICB are relatively low, and not all patients benefit from immunotherapy (7). A critical unmet need is to identify mechanisms of response and resistance and design rational combination strategies with immunotherapy (8–10). However, because of its complex and dynamic nature, our understanding of the immune response in tumor microenvironment remains incomplete (11, 12).

Squamous cell carcinomas (SCC) arise from epithelial tissues of the aerodigestive or genitourinary tracts. They are frequently found in head and neck, esophagus, lung, and cervix. SCCs share common histologic features and certain risk factors such as smoking, alcohol consumption, and human papillomavirus (HPV) infection (13). Recent TCGA studies (14, 15) have revealed that SCCs also demonstrate similar molecular patterns that are distinct from other cancer types. These studies were primarily focused on tumor cell-intrinsic characteristics such as somatic mutations (16), copy number alternations, and dysregulated pathways. Although the immune microenvironment has been recently analyzed in a pan-caner or cancer-specific settings (17–20), there are no studies that provide a comprehensive immune characterization specifically for SCCs.

In this study, we identified six robust pan-SCC immune subtypes based on consensus clustering of immune-related gene expression profiles, and further validated their reproducibility in an independent meta-cohort. We showed that each of the six immune subtypes was associated with distinct gene expression patterns, molecular and cellular characteristics, as well as clinical outcomes. Finally, our analysis revealed a complex immune landscape consisting of both discrete clusters and continuous spectrum across patients.

Patients and datasets

This study was approved by the institutional review board (IRB) and conducted in accordance with ethical guidelines such as the Declaration of Helsinki. Patient informed consent was waived given the use of existing, de-identified public datasets. For the study design, please refer to Supplementary Materials and Methods and Supplementary Fig. S1. The discovery cohort for identifying the immune subtypes consists of 1,368 patients with SCC in TCGA (Supplementary Table S1). Four major cancer types were included: head and neck squamous cell carcinoma (HNSC), lung squamous cell carcinoma (LUSC), cervical squamous cell carcinoma (CESC), and esophageal squamous cell carcinoma (ESCA). Four independent cohorts (total n = 938), each representing the single largest public gene expression dataset outside TCGA for each of the four cancer types, were used to validate the immune subtypes (Supplementary Table S1). For details about data preprocessing, please refer to Supplementary Materials and Methods.

Discovery and validation of the immune subtypes

Based on the expression of 1,989 immune-related genes (Supplementary Materials and Methods; Supplementary Table S2), we used consensus clustering (Supplementary Materials and Methods; refs. 21, 22) to identify robust clusters of patients, that is, immune subtypes (IS) and immune gene modules (GM). Then, we validated the immune subtypes in a large independent meta-cohort collected from GEO. The in-group proportion (IGP; ref. 23) and Pearson correlation among centroids of gene module scores were used to quantitatively measure the consistency in subtype identification at both patient and subtype levels in the discovery and validation cohorts (Supplementary Materials and Methods).

Assessing clinical, molecular, cellular characteristics associated with the immune subtypes

We first evaluated the distribution of the immune subtypes according to cancer type and HPV infection status. Next, we assessed the prognostic value of the immune subtypes using log-rank test and multivariable Cox regression with age, stage, cancer type, gender, and HPV infection status as covariates, using overall survival (OS) and progression-free survival (PFS) as the endpoint. Death and progression events after 60 months were censored due to a relatively short follow-up time and small number of late events in TCGA cohorts (24). The association of immune subtypes with a variety of immune-related molecular and cellular features (Supplementary Materials and Methods; ref. 17) was assessed with ANOVA (Supplementary Table S3).

Defining the immune landscape

Considering the dynamic nature of the immune system, we conducted dimensionality reduction analysis using a graph learning-based method to reveal the intrinsic structure and visualize the distribution of individual patients (Supplementary Materials and Methods). Briefly, this analysis projects the high-dimensional gene expression data to a tree structure in a low dimensional space, where the local geometric information is preserved (25). This approach was previously used to model cancer progression and define developmental trajectory using bulk and single-cell gene expression data (25–27). Here, we extend the analysis to the immune gene expression profiles. This immune landscape reflects the relationship among patients in a nonlinear manifold, which may complement the discrete immune subtypes defined in the linear Euclidean space.

After defining the immune landscape, the intra-cluster heterogeneity within immune subtype 1, 2, 4, and 6 was assessed in terms of gene module expression with ANOVA. The survival difference of three subgroups of immune subtype 4 was also compared using log-rank test.

Immune subtypes and functional gene modules

By performing consensus clustering of 1,368 tumors based on the immune-related gene expression profile, we identified six robust immune subtypes and seven gene modules in the TCGA discovery cohort (Supplementary Fig. S2). The patient-level annotation of the immune subtypes is shown in Supplementary Table S4. Gene modules appeared to be more closely clustered compared with immune subtypes. The functions of gene modules correspond to angiogenesis, inflammation, reactive stroma, T cell, IFNγ, TGFβ, and differentiation (Supplementary Table S5). Specifically, our gene module of reactive stroma was consistent with a previously proposed 25-gene stromal signature (28), in which 23 of 25 genes were included in our immune-related gene set and all 23 genes were assigned to this module. Gene module 6 was defined as TGFβ related due to its correlation with the TGFβ response score (Spearman ρ = 0.47; P < 2.2 × 10−16; ref. 29). The annotation of gene module 7 as a differentiation module was further supported by its significant correlation with histological grade (Supplementary Fig. S3).

Each of the six immune subtypes was associated with distinct immune gene expression patterns (Fig. 1A and B). Among all subtypes, immune subtype 3 (IS3) had the lowest expression in the gene modules of inflammation, T cell, and IFNγ, suggesting an immune-cold phenotype. This is closely followed by IS6, which also had a low inflammation signal but with an intermediate level of T cell and IFNγ module expression. By contrast, IS5 had the highest expression for most gene modules such as inflammation, reactive stroma, and TGFβ (except IFNγ), implying an immune-hot but suppressive microenvironment. However, IS4 had a more favorable antitumor immune response with the highest T cell and IFNγ gene expression and low reactive stroma and TGFβ, suggesting a favorable immune-activated phenotype. The remaining two subtypes IS1 and IS2 demonstrated an intermediate immune infiltration toward an immune-suppressive phenotype with elevated expression of reactive stroma and TGFβ modules.

Figure 1.

The immune subtypes and gene modules in TCGA pan-SCC cohort. A, Columns and rows represent patients and genes, respectively. Patients are arranged based on their immune subtypes, and genes are ordered based on the gene modules. Cancer type and HPV infection status are also annotated for each patient. B, Expression patterns of seven gene modules across six immune subtypes. The middle bar in each box represents the median expression level of corresponding gene module score in certain immune subtype. C, The distribution of immune subtypes among cancer types and HPV infection status in the TCGA cohort.

Figure 1.

The immune subtypes and gene modules in TCGA pan-SCC cohort. A, Columns and rows represent patients and genes, respectively. Patients are arranged based on their immune subtypes, and genes are ordered based on the gene modules. Cancer type and HPV infection status are also annotated for each patient. B, Expression patterns of seven gene modules across six immune subtypes. The middle bar in each box represents the median expression level of corresponding gene module score in certain immune subtype. C, The distribution of immune subtypes among cancer types and HPV infection status in the TCGA cohort.

Close modal

To validate our findings in TCGA cohort, we assessed reproducibility of the immune subtypes in an independent GEO meta-cohort. The gene module expression patterns were highly concordant between discovery and validation cohorts with an average linear correlation of 0.93 (Supplementary Table S6; Supplementary Fig. S4). At the individual patient level, there was moderate to good agreement between the two cohorts (IGP from IS1 to IS6: 0.61, 0.63, 0.61, 0.60, 0.72, and 0.55). All IGP values were significantly higher (P < 2 × 10−16) relative to a random partition into six groups (IGP = 0.17).

Clinical characteristics and prognoses of the immune subtypes

Each of the six immune subtypes included patients of all four types of SCC, and vice versa. The distribution of immune subtypes within each cancer type was tissue-specific (Chi-squared P < 2.2 × 10−16, Fig. 1C). For instance, a majority (∼75%) of lung SCCs were clustered into subtypes 1 and 5, whereas around 80% of cervical SCCs were clustered into subtypes 4 and 6. However, the distribution of head and neck SCCs was much more diverse. Of note, the majority of patients with HPV+ disease were clustered in two subtypes (IS4: 40.9% and IS6: 32.3%, Fig. 1C).

We observed significantly prognostic impact of the immune subtypes in TCGA cohort (P < 0.005, Fig. 2). Overall, IS4 was associated with the best prognosis for both OS and PFS. By contrast, the immune-cold subtype IS3 was associated with the worst outcomes among all subtypes. This survival difference was independent of cancer type, stage, age, gender, and HPV infection status (Table 1). The remaining subtypes had intermediate prognoses. Similar patterns were observed in TCGA CESC and HNSC cohorts (Supplementary Fig. S5A and S5B). In TCGA LUSC dataset, IS6 and IS4 had a marginally better OS compared with others (P = 0.071; Supplementary Fig. S5C). No clear pattern was observed in TCGA ESCA dataset, likely due to a limited sample size (Supplementary Fig. S5D). Of note, among HPV+ patients, there was significant difference in OS between two major subtypes (IS6 vs. IS4: HR = 2.04; 95% CI, 1.06–3.95; P = 0.034; Supplementary Fig. S6). The prognostic effect of gene modules was shown in Supplementary Fig. S7. Consistent with previous studies, we found that a higher expression score of T-cell module was significantly associated with improved survival, whereas higher scores of reactive stroma and TGFβ modules were related to inferior survival.

Figure 2.

Five-year Kaplan–Meier curves for (A) OS and (B) PFS of all patients stratified by the immune subtypes. P value was calculated by the log-rank test among subtypes.

Figure 2.

Five-year Kaplan–Meier curves for (A) OS and (B) PFS of all patients stratified by the immune subtypes. P value was calculated by the log-rank test among subtypes.

Close modal
Table 1.

Multivariable Cox regression analysis of OS and PFS, including immune subtype, cancer types, stage, age, gender, and HPV infection status

OSPFS
VariablesHR (95% CI)PHR (95% CI)P
Immune subtypes 
 1 1.40 (0.95–2.06) 0.087* 1.19 (0.81–1.75) 0.38 
 2 1.43 (0.95–2.16) 0.085* 1.15 (0.77–1.73) 0.50 
 3 1.63 (1.08–2.48) 0.021** 1.52 (1.01–2.29) 0.046** 
 4 1.00 – 1.00 – 
 5 1.35 (0.90–2.03) 0.14 1.13 (0.75–1.71) 0.56 
 6 1.55 (1.02–2.35) 0.040** 1.35 (0.89–2.04) 0.15 
Cancer type 
 CESC 1.00 – 1.00 – 
 ESCA 1.16 (0.65–2.07) 0.62 1.54 (0.88–2.67) 0.13 
 HNSC 0.57 (0.35–0.91) 0.020** 0.57 (0.35–0.92) 0.022** 
 LUSC 0.83 (0.51–1.36) 0.46 0.69 (0.42–1.14) 0.15 
Age 1.02 (1.01–1.03) 0.00096*** 1.01 (1.00–1.02) 0.089* 
Stage 
 I 1.00 – 1.00 – 
 II 1.23 (0.93–1.63) 0.14 1.34 (0.99–1.81) 0.059* 
 III 1.78 (1.33–2.38) 0.00012*** 1.75 (1.27–2.42) 0.00071*** 
 IV 2.70 (1.92–3.82) 1.58 × 10−8*** 2.88 (2.01–4.14) 1.00 × 10−8*** 
Gender (male) 1.00 (0.80–1.25) 0.97 1.09 (0.85–1.39) 0.51 
HPV+ 0.54 (0.37–0.79) 0.0017** 0.57 (0.39–0.83) 0.0030*** 
OSPFS
VariablesHR (95% CI)PHR (95% CI)P
Immune subtypes 
 1 1.40 (0.95–2.06) 0.087* 1.19 (0.81–1.75) 0.38 
 2 1.43 (0.95–2.16) 0.085* 1.15 (0.77–1.73) 0.50 
 3 1.63 (1.08–2.48) 0.021** 1.52 (1.01–2.29) 0.046** 
 4 1.00 – 1.00 – 
 5 1.35 (0.90–2.03) 0.14 1.13 (0.75–1.71) 0.56 
 6 1.55 (1.02–2.35) 0.040** 1.35 (0.89–2.04) 0.15 
Cancer type 
 CESC 1.00 – 1.00 – 
 ESCA 1.16 (0.65–2.07) 0.62 1.54 (0.88–2.67) 0.13 
 HNSC 0.57 (0.35–0.91) 0.020** 0.57 (0.35–0.92) 0.022** 
 LUSC 0.83 (0.51–1.36) 0.46 0.69 (0.42–1.14) 0.15 
Age 1.02 (1.01–1.03) 0.00096*** 1.01 (1.00–1.02) 0.089* 
Stage 
 I 1.00 – 1.00 – 
 II 1.23 (0.93–1.63) 0.14 1.34 (0.99–1.81) 0.059* 
 III 1.78 (1.33–2.38) 0.00012*** 1.75 (1.27–2.42) 0.00071*** 
 IV 2.70 (1.92–3.82) 1.58 × 10−8*** 2.88 (2.01–4.14) 1.00 × 10−8*** 
Gender (male) 1.00 (0.80–1.25) 0.97 1.09 (0.85–1.39) 0.51 
HPV+ 0.54 (0.37–0.79) 0.0017** 0.57 (0.39–0.83) 0.0030*** 

NOTE: Immune subtype 4 was used as the baseline for survival risk comparison for immune subtype variable. CESC was used as the baseline for survival risk comparison for cancer type variable. Stage I was used as the baseline for survival risk comparison for stage variable. *, P < 0.1; **, P < 0.05; ***, P < 0.01.

Molecular and cellular characteristics of the immune subtypes

We assessed the relation between the immune subtypes and 57 previously defined immune-related molecular features. Consistent with an immune-cold phenotype, tumors in IS3 had the least leukocyte and stromal fraction (Supplementary Fig. S8), the lowest lymphocyte infiltration signature score (30) and TCR diversity (Fig. 3A and B). Interestingly, IS3 was also associated with the highest degree of somatic copy number variation such as a high aneuploidy score (Fig. 3C; Supplementary Fig. S8). A closely related subtype is IS6, which also had low leucocyte and stromal fraction (Supplementary Fig. S8). However, compared with IS3, IS6 had an increased percentage of lymphocytes including CD8+ T cells and activated natural killer (NK) cells (Fig. 3D and E; Supplementary Fig. S9). The CD8+ T-cell repertoire still appeared to be highly restricted (Fig. 3B). In addition, IS6 expressed relatively high proliferation (17) and wound healing scores (31) as well as the lowest aneuploidy and TGFβ response score (Supplementary Fig. S8; Fig. 3C and F; ref. 29).

Figure 3.

A–I, Molecular and cellular characteristics associated with the immune subtypes. The middle bar in each box represents the median level of corresponding features in certain immune subtype. The FDR-adjusted P values for all features were less than 0.05.

Figure 3.

A–I, Molecular and cellular characteristics associated with the immune subtypes. The middle bar in each box represents the median level of corresponding features in certain immune subtype. The FDR-adjusted P values for all features were less than 0.05.

Close modal

Contrary to immune-cold IS3, IS5 had the highest leukocyte fraction, lymphocyte infiltration signature score, and TCR diversity, which is consistent with an immune-hot phenotype (Fig. 3A and B; Supplementary Fig. S8). Notably, IS5 was associated with the highest macrophage regulation score with macrophages consisting more than 40% of the leukocyte infiltration (Fig. 3G; Supplementary Fig. S8). However, the majority of tumor-associated macrophages tended to have a pro-tumor M2-polarized phenotype (Fig. 3H). Additionally, IS5 had one of the highest TGFβ response score (Fig. 3F). These data suggest that the tumor microenvironment of IS5 was immune-hot but highly immune-suppressive.

Different from the extreme cold/hot immune microenvironment in the above subtypes, IS4 demonstrated a more balanced and favorable immune profile. Similar to IS5, IS4 also had high leukocyte fraction, lymphocyte infiltration signature score, and a diverse TCR repertoire (Fig. 3A and B; Supplementary Fig. S8). One major difference from IS5 is that the immune composition in IS4 was enriched with lymphocytes, such as activated CD4+ memory T cells, Tfh cell, CD8+ T cell, and activated NK cells (Fig. 3D and E; Supplementary Fig. S9). Not surprisingly, IS4 had the highest local cytolytic activity (32), suggesting a preexisting antitumor immune microenvironment (Fig. 3I). Of note, IS4 was associated with the highest M1-to-M2 ratio (Fig. 3H). In addition, IS4 had an elevated level of IFNγ response (33) and suppressed TGFβ response (Supplementary Fig. S9; Fig. 3F). A low genome aberration was also observed in IS4 (Fig. 3C; Supplementary Fig. S8).

The remaining two immune subtypes (IS1 and IS2) were more diverse with intermediate levels of immune features. Both IS1 and IS2 had relatively high stromal fraction, with intermediate-to-low lymphocyte infiltration, TCR diversity, and cytolytic score (Fig. 3A, B, and I; Supplementary Fig. S8). IS1 had a high percentage of tumor-associated macrophages with an M2-polarized phenotype (Fig. 3H). IS1 also tended to be biased towards humoral immunity with high percentages of naïve B cells and plasma cells (Supplementary Fig. S9). Additionally, IS1 had high genome aberration, CTA score, and TGFβ response score (Fig. 3C and F; Supplementary Fig. S8). However, a low genome aberration and tumor mutation burden (16), and the highest IFNγ response score were observed in IS2 (Supplementary Fig. S8).

Immune landscape of SCC

To facilitate visualization and uncover the underlying structures of the distribution of individual patients, we applied the graph learning-based dimensionality reduction technique (25, 27) to the immune gene expression profiles. This analysis cast individual patients into a manifold with sparse tree structures, and defined the immune landscape of SCC (Fig. 4A). Consistent with our previously defined immune subtypes, many patients were segregated into distinct clusters, for example, subtypes 3, 5, and 6. The location of individual patients in the immune landscape represents the overall characteristics of the tumor immune microenvironment in the corresponding subtype. For instance, the immune-hot subtype IS5 and immune-cold subtype IS3 were distributed at the opposite end of the horizontal axis in the immune landscape. Therefore, we hypothesized that the horizontal axis in the immune landscape represents the overall immune infiltration. Indeed, the horizontal coordinate was highly correlated with the inflammation and T-cell modules (ρ = 0.91 and 0.78, respectively; both P < 2.2 × 10−16). However, the vertical coordinate of the immune landscape appeared to be more complex and may reflect multiple biological processes, mainly the reactive stroma and TGFβ gene modules (ρ = 0.61 and 0.62, respectively; both P < 2.2 × 10−16).

Figure 4.

The immune landscape of SCC and the intracluster heterogeneity within immune subtype 4. A, The immune landscape of SCC: each point represents a patient with colors corresponding to the immune subtype defined previously. B, Patients of immune subtype 4 could be further stratified into three subgroups based on their location in the immune landscape. C, Gene module expression patterns were shown to illustrate the intracluster heterogeneity of immune subtype 4. D, The three subgroups of patients in immune subtype 4 as stratified by the immune landscape were associated with distinct prognoses. Log-rank P value was calculated among subgroup stratification.

Figure 4.

The immune landscape of SCC and the intracluster heterogeneity within immune subtype 4. A, The immune landscape of SCC: each point represents a patient with colors corresponding to the immune subtype defined previously. B, Patients of immune subtype 4 could be further stratified into three subgroups based on their location in the immune landscape. C, Gene module expression patterns were shown to illustrate the intracluster heterogeneity of immune subtype 4. D, The three subgroups of patients in immune subtype 4 as stratified by the immune landscape were associated with distinct prognoses. Log-rank P value was calculated among subgroup stratification.

Close modal

The immune landscape further revealed significant intracluster heterogeneity within each subtype. We observed that certain subtypes appeared to be more diverse and heterogeneous than others. For instance, IS1 could be further divided into three subgroups based on their location in the immune landscape, which showed different immune gene expression profiles in specific modules (Supplementary Fig. S10A). Similar results were obtained for IS2 and IS6 (Supplementary Fig. S10B and S10C). Interestingly, the three subgroups of patients in IS4 as stratified by the immune landscape (Fig. 4B) were associated with distinct gene expression patterns and prognoses (Fig. 4C and D), and the same survival pattern was observed within cervical and head and neck cancer specifically (Supplementary Fig. S11). Within IS4, the subset of patients with the best prognosis (IS4A) was associated with the highest expression of T-cell module. Although these patients were located very close to IS5 in the immune landscape, their prognoses were rather different from IS5 (Fig. 4D). These results indicate that our immune landscape analysis provided complementary value to previously identified immune subtypes.

Immunotherapy is being used to treat an increasing number of cancers including SCCs in clinical practice or ongoing trials (4, 5, 34–36). However, response and survival benefit is typically limited to a subset of patients. A better understanding of the tumor immune microenvironment is needed for designing novel immunotherapies to improve response and outcomes. Here, we present the identification, independent validation, and comprehensive characterization of six reproducible immune subtypes of SCCs in a multicohort retrospective study. We found that each of the immune subtypes was associated with distinct gene expression profiles, and accordingly demonstrated widely different patterns in tumor genetic aberrations, cytokine profiles, tumor-infiltrating immune cell composition, functional orientation, and, importantly, clinical outcomes. This study provides a conceptual framework to understand the immune response of SCCs, and may have clinical implications for personalized immunotherapy.

Our work differs from recent immune landscape studies in several important aspects. First, we focused on squamous carcinomas that have the same histology and similar risk factors, for example, viral infection or exposures to exogenous carcinogens. In a recent pan-cancer analysis, Thorsson and colleagues discovered (17) six immune classes across 33 cancer types of various tissues and etiology, where the vast majority (90%–95%) of SCCs fell into two subtypes (C1: wound healing and C2: IFNγ dominant) with almost identical prognoses. Our analysis further stratified these patients within the C1/C2 classes (Supplementary Table S7). Second, instead of using established signatures, we carefully curated a comprehensive set of genes reflecting various immunological processes. Third, we assessed the reproducibility of our immune subtypes in an independent cohort. Finally, our study further extends beyond previous studies for patient subtyping based on simple clustering analyses. The discrete subtype information failed to capture inter- and intracluster relationships and did not provide the overall structure of the patient distribution. To remedy those shortcomings, we applied graph learning approaches to uncover the tree structures of immune profiles among patients, which provided complementary information to clustering analysis and offered new insight into the complex immune landscape of SCC.

In a recent study specifically focused on head and neck cancer, Chen and colleagues proposed three immune classes, that is, active, exhausted, and non-immune class (19). Consistent with these findings, all patients in our immune-cold subtype IS3 were classified as the non-immune class defined by Chen and colleagues (19) in TCGA cohort (Supplementary Table S8). Similarly, the vast majority (91%) of patients in our immune “favorable” IS4 belonged to the active immune class. However, the distribution for other immune subtypes was more heterogeneous with respect to the three immune classes. There was no dominant immune subtype within immune classes, and each class consisted of patients from at least four immune subtypes. Therefore, our study provides a different perspective of the complex immune landscape in head and neck cancer that complements previous analyses.

The impact of the tumor immune microenvironment on patient survival is well established in multiple cancer types (37). In our study, squamous carcinomas of IS4 demonstrated the highest levels of infiltration by immune effectors such as CD8+ T and activated NK cells, and elevated expression of IFNγ response. Accordingly, patients in subtype 4 had the best prognosis. In comparison, tumors of subtype 6 had reduced lymphocyte infiltration and thus these patients had relatively worse outcomes. Importantly, the same survival pattern was observed within the subgroup of patients with HPV-associated head and neck cancer and cervical cancer. These suggest that the immune profile may be a key determinant of outcomes across cancer types and could potentially be incorporated into future biomarker-based risk stratification strategy for personalized therapy of HPV-associated cancer. Finally, the immune-cold tumors of subtype 3 had minimal lymphocyte infiltration, which is likely due to its high level of aneuploidy (38). Consequently, these patients appeared to have the worst survival. These results are consistent with previous studies showing that preexisting antitumor immunity generally results in improved prognosis across cancers (32).

The relative dominance between immune stimulatory and suppressive factors is critical in determining prognosis. Although tumors of subtype 5 demonstrated a high level of immune infiltration similar to subtype 4, its immune composition was dominated by M2-polarized macrophages with highly immune-suppressive factors such as TGFβ signaling and reactive stroma. Accordingly, patients in subtype 5 had significantly worse prognosis compared with subtype 4. Similar patterns were also observed for tumors of subtypes 1 and 2, which had high macrophage/lymphocyte ratio, and low infiltration by CD8+ T cells and activated NK cells, leading to an inferior prognosis. These data add to the accumulating evidence that the immune composition, functional orientation (39), and immune-suppressive mechanisms such as TGFβ signaling (40) play critical roles in determining therapeutic response and outcomes.

Traditionally, an individual-based model is often used to develop predictive and prognostic biomarkers, which requires the response to therapy and clinical outcomes to be known for each individual patient. By contrast, our approach is “unsupervised,” which relies on the immune-related gene expression profiles to reveal the underlying structures of the immune landscape within tumors. In future, the intrinsic properties of immune landscape may be incorporated when developing biomarkers. It is conceivable that a hierarchical model, first stratifying patients into subgroups and then applying individual-based risk stratification, might be used to predict clinical outcomes with biological relevance. The idea of “subtype-specific” biomarkers has been successfully applied to improve outcome prediction in breast and colorectal cancer (41, 42). Therefore, integrating subtype analyses and individual-based model could be a promising approach to developing clinically relevant biomarkers.

Although the immune landscape by and large recapitulated the immune subtypes based on clustering analyses, it also uncovered previously unappreciated intracluster heterogeneity with potential clinical relevance. For instance, a fraction of patients in IS4 were shown to have a superior prognosis relative to others in the same subtype. This insight could be only gained by integrating information from both immune subtyping and landscape analyses. Similarly, subtypes 1 and 2 showed divergent behaviors with intermediate transition states in the immune landscape. This raises the intriguing question of how to optimally modulate the host immune response so that patients are mobilized toward more favorable states, providing a roadmap to more successful immunotherapy.

Our study has potential therapeutic implications for the rational design of combination immunotherapy strategies. For patients with a favorable immune microenvironment (e.g., subtype 4), ICB may be used to enhance the preexisting antitumor immunity of these patients and further improve their survival. However, for patients in other subtypes, ICB alone may be insufficient due to the suboptimal immune activation or presence of immune-suppressive mechanisms. Combination of ICB with immune costimulatory antibodies such as OX-40 and 4-1BB may be used to amplify or boost the dampened immune response for patients in subtype 6. However, for those with an immune-cold tumor (i.e., subtype 3), ICB should be optimally combined with cancer vaccines, oncolytic viruses, radiotherapy, chemotherapy, which may convert non-inflamed tumors into inflamed tumors by triggering an inflammatory response (43). For the remaining patients in subtypes 1, 2, and 5, depending on their specific immune and stromal microenvironment, macrophage (44) or NK-cell targeted (45, 46) therapies, CAF-targeted therapy (47), anti-TGFβ (48), or anti-angiogenic (49) therapies might be used together with ICB to revert the ineffective antitumor immune response.

In summary, we identified 6 reproducible immune subtypes of SCCs with distinct molecular characteristics and clinical outcomes. Our study provides a conceptual framework for the future design of rational combination treatment strategies and optimal selection of patients to improve immunotherapy outcomes.

J.B. Sunwoo holds ownership interest (including patents) in and is a consultant/advisory board member for Indapta Therapeutics. No potential conflicts of interest were disclosed by the other authors.

Conception and design: B. Li, Y. Cui, J.B. Sunwoo, R. Li

Development of methodology: B. Li, Y. Cui, R. Li

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): B. Li

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): B. Li, Y. Cui, D.K. Nambiar, J.B. Sunwoo, R. Li

Writing, review, and/or revision of the manuscript: B. Li, Y. Cui, D.K. Nambiar, J.B. Sunwoo, R. Li

Study supervision: Y. Cui, R. Li

This project is partially supported by the NIH grant 1R01 CA222512-01.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Pardoll
DM
. 
The blockade of immune checkpoints in cancer immunotherapy
.
Nat Rev Cancer
2012
;
12
:
252
64
.
2.
Topalian
SL
,
Hodi
FS
,
Brahmer
JR
,
Gettinger
SN
,
Smith
DC
,
McDermott
DF
, et al
Safety, activity, and immune correlates of anti–PD-1 antibody in cancer
.
New Engl J Med
2012
;
366
:
2443
54
.
3.
Brahmer
JR
,
Tykodi
SS
,
Chow
LQ
,
Hwu
WJ
,
Topalian
SL
,
Hwu
P
, et al
Safety and activity of anti-PD-L1 antibody in patients with advanced cancer
.
N Engl J Med
2012
;
366
:
2455
65
.
4.
Brahmer
J
,
Reckamp
KL
,
Baas
P
,
Crinò
L
,
Eberhardt
WE
,
Poddubskaya
E
, et al
Nivolumab versus docetaxel in advanced squamous-cell non–small-cell lung cancer
.
New Engl J Med
2015
;
373
:
123
35
.
5.
Garon
EB
,
Rizvi
NA
,
Hui
R
,
Leighl
N
,
Balmanoukian
AS
,
Eder
JP
, et al
Pembrolizumab for the treatment of non–small-cell lung cancer
.
New Engl J Med
2015
;
372
:
2018
28
.
6.
Wolchok
JD
,
Chiarion-Sileni
V
,
Gonzalez
R
,
Rutkowski
P
,
Grob
JJ
,
Cowey
CL
, et al
Overall survival with combined nivolumab and ipilimumab in advanced melanoma
.
New Engl J Med
2017
;
377
:
1345
56
.
7.
Ribas
A
,
Wolchok
JD
. 
Cancer immunotherapy using checkpoint blockade
.
Science
2018
;
359
:
1350
5
.
8.
Sharma
P
,
Hu-Lieskovan
S
,
Wargo
JA
,
Ribas
A
. 
Primary, adaptive, and acquired resistance to cancer immunotherapy
.
Cell
2017
;
168
:
707
23
.
9.
Gibney
GT
,
Weiner
LM
,
Atkins
MB
. 
Predictive biomarkers for checkpoint inhibitor-based immunotherapy
.
Lancet Oncol
2016
;
17
:
e542
51
.
10.
Ayers
M
,
Lunceford
J
,
Nebozhyn
M
,
Murphy
E
,
Loboda
A
,
Kaufman
DR
, et al
IFN-γ–related mRNA profile predicts clinical response to PD-1 blockade
.
J Clin Invest
2017
;
127
:
2930
40
.
11.
Chen
DS
,
Mellman
I.
. 
Elements of cancer immunity and the cancer–immune set point
.
Nature
2017
;
541
:
321
30
.
12.
Binnewies
M
,
Roberts
EW
,
Kersten
K
,
Chan
V
,
Fearon
DF
,
Merad
M
, et al
Understanding the tumor immune microenvironment (TIME) for effective therapy
.
Nat Med
2018
;
24
:
541
50
.
13.
Dotto
GP
,
Rustgi
AK
. 
Squamous cell cancers: a unified perspective on biology and genetics
.
Cancer Cell
2016
;
29
:
622
37
.
14.
Campbell
JD
,
Yau
C
,
Bowlby
R
,
Liu
Y
,
Brennan
K
,
Fan
H
, et al
Genomic, pathway network, and immunologic features distinguishing squamous carcinomas
.
Cell Rep
2018
;
23
:
194
212
.
e6
.
15.
Hoadley
KA
,
Yau
C
,
Hinoue
T
,
Wolf
DM
,
Lazar
AJ
,
Drill
E
, et al
Cell-of-origin patterns dominate the molecular classification of 10,000 tumors from 33 types of cancer
.
Cell
2018
;
173
:
291
304
.
e6
.
16.
Riaz
N
,
Havel
JJ
,
Makarov
V
,
Desrichard
A
,
Urba
WJ
,
Sims
JS
, et al
Tumor and microenvironment evolution during immunotherapy with nivolumab
.
Cell
2017
;
171
:
934
49
.
e16
.
17.
Thorsson
V
,
Gibbs
DL
,
Brown
SD
,
Wolf
D
,
Bortone
DS
,
Ou Yang
TH
, et al
The immune landscape of cancer
.
Immunity
2018
;
48
:
812
30
.
e14
.
18.
Tamborero
D
,
Rubio-Perez
C
,
Muiños
F
,
Sabarinathan
R
,
Piulats
JM
,
Muntasell
A
, et al
A pan-cancer landscape of interactions between solid tumors and infiltrating immune cell populations
.
Clin Cancer Res
2018
;
24
:
3717
28
.
19.
Chen
Y
,
Wang
Y
,
Lv
J
,
Li
Y
,
Chua
M
,
Le
Q
, et al
Identification and validation of novel microenvironment-based immune molecular subgroups of head and neck squamous cell carcinoma: implications for immunotherapy
.
Ann Oncol
2019
;
30
:
68
75
.
20.
Keck
MK
,
Zuo
Z
,
Khattri
A
,
Stricker
TP
,
Brown
CD
,
Imanguli
M
, et al
Integrative analysis of head and neck cancer identifies two biologically distinct HPV and three non-HPV subtypes
.
Clin Cancer Res
2015
;
21
:
870
81
.
21.
Monti
S
,
Tamayo
P
,
Mesirov
J
,
Golub
T
. 
Consensus clustering: a resampling-based method for class discovery and visualization of gene expression microarray data
.
Machine learning
2003
;
52
:
91
118
.
22.
Wilkerson
MD
,
Hayes
DN
. 
ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking
.
Bioinformatics
2010
;
26
:
1572
3
.
23.
Kapp
AV
,
Tibshirani
R
. 
Are clusters found in one dataset present in another dataset?
Biostatistics
2007
;
8
:
9
31
.
24.
Cancer Genome Atlas Research Network
,
Albert Einstein College of Medicine
,
Analytical Biological Services
,
Barretos Cancer Hospital
,
Baylor College of Medicine
,
Beckman Research Institute of City of Hope
et al 
Integrated genomic and molecular characterization of cervical cancer
.
Nature
2017
;
543
:
378
84
.
25.
Mao
Q
,
Wang
L
,
Goodison
S
,
Sun
Y
.
Dimensionality reduction via graph structure learning
. In:
KDD '15 Proceedings of the 21st ACM SIGKDD International Conference on Knowledge Discovery and Data Mining; 2015 Aug 10–13
;
Sydney, Australia. New York
:
ACM
; 
2015
. p.
765
74
.
26.
Trapnell
C
,
Cacchiarelli
D
,
Grimsby
J
,
Pokharel
P
,
Li
S
,
Morse
M
, et al
The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells
.
Nat Biotechnol
2014
;
32
:
381
6
.
27.
Sun
YJ
,
Yao
J
,
Nowak
NJ
,
Goodison
S
. 
Cancer progression modeling using static sample data
.
Genome Biol
2014
;
15
:
440
57
.
28.
Moffitt
RA
,
Marayati
R
,
Flate
EL
,
Volmar
KE
,
Loeza
SG
,
Hoadley
KA
, et al
Virtual microdissection identifies distinct tumor- and stroma-specific subtypes of pancreatic ductal adenocarcinoma
.
Nat Genet
2015
;
47
:
1168
78
.
29.
Teschendorff
AE
,
Gomez
S
,
Arenas
A
,
El-Ashry
D
,
Schmidt
M
,
Gehrmann
M
, et al
Improved prognostic classification of breast cancer defined by antagonistic activation patterns of immune response pathway modules
.
BMC Cancer
2010
;
10
:
604
24
.
30.
Calabro
A
,
Beissbarth
T
,
Kuner
R
,
Stojanov
M
,
Benner
A
,
Asslaber
M
, et al
Effects of infiltrating lymphocytes and estrogen receptor on gene expression and prognosis in breast cancer
.
Breast Cancer Res Treat
2009
;
116
:
69
77
.
31.
Chang
HY
,
Sneddon
JB
,
Alizadeh
AA
,
Sood
R
,
West
RB
,
Montgomery
K
, et al
Gene expression signature of fibroblast serum response predicts human cancer progression: similarities between tumors and wounds
.
PLoS Biol
2004
;
2
:
e7
.
32.
Rooney
MS
,
Shukla
SA
,
Wu
CJ
,
Getz
G
,
Hacohen
N
. 
Molecular and genetic properties of tumors associated with local immune cytolytic activity
.
Cell
2015
;
160
:
48
61
.
33.
Wolf
DM
,
Lenburg
ME
,
Yau
C
,
Boudreau
A
,
van‘t Veer
LJ
. 
Gene co-expression modules as clinically relevant hallmarks of breast cancer diversity
.
PLoS One
2014
;
9
:
e88309
.
34.
Ferris
RL
,
Blumenschein
G
 Jr
,
Fayette
J
,
Guigay
J
,
Colevas
AD
,
Licitra
L
, et al
Nivolumab for recurrent squamous-cell carcinoma of the head and neck
.
N Engl J Med
2016
;
375
:
1856
67
.
35.
Doi
T
,
Piha-Paul
SA
,
Jalal
SI
,
Saraf
S
,
Lunceford
J
,
Koshiji
M
, et al
Safety and antitumor activity of the anti-programmed death-1 antibody pembrolizumab in patients with advanced esophageal carcinoma
.
J Clin Oncol
2018
;
36
:
61
7
.
36.
Chung
HC
,
Schellens
JHM
,
Delord
JP
,
Perets
R
,
Italiano
A
,
Shapira-Frommer
R
, et al
Pembrolizumab treatment of advanced cervical cancer: updated results from the phase 2 KEYNOTE-158 study
.
J Clin Oncol
2018
;
36
(
15 suppl
):
5522
.
37.
Li
B
,
Cui
Y
,
Diehn
M
,
Li
R
. 
Development and validation of an individualized immune prognostic signature in early-stage nonsquamous non–small cell lung cancer
.
JAMA Oncol
2017
;
3
:
1529
37
.
38.
Davoli
T
,
Uno
H
,
Wooten
EC
,
Elledge
SJ
. 
Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy
.
Science
2017
;
355
:
eaaf8399
413
.
39.
Fridman
WH
,
Zitvogel
L
,
Sautès–Fridman
C
,
Kroemer
G
. 
The immune contexture in cancer prognosis and treatment
.
Nat Rev Clin Oncol
2017
;
14
:
717
34
.
40.
Mariathasan
S
,
Turley
SJ
,
Nickles
D
,
Castiglioni
A
,
Yuen
K
,
Wang
Y
, et al
TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells
.
Nature
2018
;
554
:
544
8
.
41.
Callari
M
,
Cappelletti
V
,
D'Aiuto
F
,
Musella
V
,
Lembo
A
,
Petel
F
, et al
Subtype-specific metagene-based prediction of outcome after neoadjuvant and adjuvant treatment in breast cancer
.
Clin Cancer Res
2016
;
22
:
337
45
.
42.
Bramsen
JB
,
Rasmussen
MH
,
Ongen
H
,
Mattesen
TB
,
Orntoft
MW
,
Arnadottir
SS
, et al
Molecular-subtype-specific biomarkers improve prediction of prognosis in colorectal cancer
.
Cell Rep
2017
;
19
:
1268
80
.
43.
Ott
PA
,
Hodi
FS
,
Kaufman
HL
,
Wigginton
JM
,
Wolchok
JD
. 
Combination immunotherapy: a road map
.
J Immunother Cancer
2017
;
5
:
16
31
.
44.
Mantovani
A
,
Marchesi
F
,
Malesci
A
,
Laghi
L
,
Allavena
P
. 
Tumour-associated macrophages as treatment targets in oncology
.
Nat Rev Clin Oncol
2017
;
14
:
399
416
.
45.
Guillerey
C
,
Huntington
ND
,
Smyth
MJ
. 
Targeting natural killer cells in cancer immunotherapy
.
Nat Immunol
2016
;
17
:
1025
36
.
46.
Barry
KC
,
Hsu
J
,
Broz
ML
,
Cueto
FJ
,
Binnewies
M
,
Combes
AJ
, et al
A natural killer–dendritic cell axis defines checkpoint therapy–responsive tumor microenvironments
.
Nat Med
2018
;
24
:
1178
91
.
47.
Hanley
CJ
,
Mellone
M
,
Ford
K
,
Thirdborough
SM
,
Mellows
T
,
Frampton
SJ
, et al
Targeting the myofibroblastic cancer-associated fibroblast phenotype through inhibition of NOX4
.
J Natl Cancer Inst
2018
;
110
20
.
48.
Holmgaard
RB
,
Schaer
DA
,
Li
Y
,
Castaneda
SP
,
Murphy
MY
,
Xu
X
, et al
Targeting the TGFβ pathway with galunisertib, a TGFβRI small molecule inhibitor, promotes anti-tumor immunity leading to durable, complete responses, as monotherapy and in combination with checkpoint blockade
.
J Immunother Cancer
2018
;
6
:
47
62
.
49.
Khan
KA
,
Kerbel
RS
. 
Improving immunotherapy outcomes with anti-angiogenic treatments and vice versa
.
Nat Rev Clin Oncol
2018
;
15
:
310
24
.

Supplementary data