Abstract
To identify immune subtypes and investigate the immune landscape of squamous cell carcinomas (SCC), which share common etiology and histologic features.
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.
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.
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
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.
Introduction
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.
Materials and Methods
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.
Results
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.
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.
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.
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.
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.
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.
Multivariable Cox regression analysis of OS and PFS, including immune subtype, cancer types, stage, age, gender, and HPV infection status
. | OS . | PFS . | ||
---|---|---|---|---|
Variables . | HR (95% CI) . | P . | HR (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*** |
. | OS . | PFS . | ||
---|---|---|---|---|
Variables . | HR (95% CI) . | P . | HR (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).
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.
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.
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).
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.
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.
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.
Discussion
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.
Disclosure of Potential Conflicts of Interest
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.
Authors' Contributions
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
Acknowledgments
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.