Systematic Profiling of Immune Risk Model to Predict Survival and Immunotherapy Response in Head and Neck Squamous Cell Carcinoma

Background and Purpose Head and neck squamous carcinoma (HNSCC), characterized by immunosuppression, is a group of highly heterogeneous cancers. Although immunotherapy exerts a promising influence on HNSCC, the response rate remains low and varies in assorted primary sites. Immunological mechanisms underlying HNSCC pathogenesis and treatment response are not fully understood. This study aimed to develop a differentially expressed genes (DEGs)–based risk model to predict immunotherapy efficacy and stratify prognosis of HNSCC patients. Materials and Methods The expression profiles of HNSCC patients were downloaded from The Cancer Genome Atlas (TCGA) database. The tumor microenvironment and immune response were estimated by cell type identification via estimating relative subset of known RNA transcripts (CIBERSORT) and immunophenoscore (IPS). The differential expression pattern based on human papillomavirus status was identified. A DEGs-based prognostic risk model was developed and validated. All statistical analyses were performed with R software (version 3.6.3). Results By using the TCGA database, we identified DKK1, HBEGF, RNASE7, TNFRSF12A, INHBA, and IPIK3R3 as DEGs that were associated with patients’ overall survival (OS). Patients were stratified into the high- and low-risk subgroups according to a DEGs-based prognostic risk model. Significant difference in OS was found between the high- and low-risk patients (1.64 vs. 2.18 years, P = 0.0017). In multivariate Cox analysis, the risk model was an independent prognostic factor for OS (hazard radio = 1.06, 95% confidence interval [1.02–1.10], P = 0.004). More CD8+ T cells and regulatory T cells were observed in the low-risk group and associated with a favorable prognosis. The IPS analysis suggested that the low-risk patients possessed a higher IPS score and a higher immunoreactivity phenotype, which were correlated with better immunotherapy response. Conclusion Collectively, we established a reliable DEGs-based risk model with potential prognostic value and capacity to predict the immunophenotype of HNSCC patients.


INTRODUCTION
Squamous cell carcinomas, originating from the oral cavity, oropharynx, larynx, or hypopharynx, are collectively referred to as head and neck squamous cell carcinoma (HNSCC). HNSCC has become the sixth leading cancer by incidence worldwide, whereas only fewer than 50% of patients could survive for 5 years (Kamangar et al., 2006). About two-thirds of HNSCC patients are diagnosed at late stages with poor prognosis. Patients with locally advanced stage HNSCC are treated with surgery and postoperative radiotherapy. Targeted drugs, such as cetuximab, an epidermal growth factor receptor (EGFR)-specific antibody, have been utilized for HNSCC and reached a limited response rate, possibly due to its clinical heterogeneity (Licitra et al., 2011). Novel therapeutic strategies are in urgent need.
Recent studies documented that the programmed death 1 (PD-1) inhibitors such as pembrolizumab and nivolumab significantly improve the overall survival (OS) of recurrent/metastatic HNSCC patients. However, PD-1 inhibitors have merely low to moderate response rates (Chow et al., 2016;Seiwert et al., 2016). The factors that determine treatment response to PD-1 inhibitors in HNSCC are not clear. The tumor microenvironment (TME) of HNSCC is highly heterogeneous and predominantly immunosuppressive, characterized by macrophages and myeloid-derived suppressor cell recruitment (Canning et al., 2019), T cell and natural killer (NK) cell dysfunction, regulatory T cell (Treg) activation (Reichert et al., 2002), and alteration in cytokine release such as enhanced interleukin-10 (IL-10) and IL-6 production and reduced transforming growth factor β (TGF-β) and IL-12 secretion (Lathers and Young, 2004;Varilla et al., 2013). These changes in TME might lead to a limited PD-1 inhibitor treatment response. On the other hand, high-level expression of PD-1 and/or cytotoxic T-lymphocyte antigen 4 (CTLA4) on T cells, as well as programmed death ligand 1 (PD-L1) up-regulation on malignant cells and immune cells in a portion of HNSCC patients, has been reported, which might dictate the treatment efficacy of immune checkpoint inhibitors (ICIs) (Badoual et al., 2013;Concha-Benavente et al., 2016;Mandal et al., 2016). There is still a lack of a useful risk model that could predict the response to immunotherapy and provide prognostic information on HNSCC patients. Thus, deeper understanding of genomic alterations and immune-related markers might help develop such risk models.
Increased human papillomavirus (HPV) infection together with tobacco and alcohol abuse has been identified as the most important risk factors. HPV infection is mostly associated with oropharyngeal cancer (Tanaka and Alawi, 2018). A recent study in France showed 43.1% HPV-positive rate in HNSCC patients. Compared with HPV-negative HNSCC patients, HPV-positive HNSCC patients have a better OS and disease-free survival (Mirghani et al., 2019). HPV viral genome could integrate into the host genome and alter gene expression, affecting HNSCC TME. For example, HPV infection significantly impairs IL-6 and macrophage colony-stimulating factor release, creating an immunosuppressive microenvironment (Smola, 2017). Besides, HPV-positive HNSCC has high immune infiltrates such as Tregs and CD56 + NK cells (Mandal et al., 2016). Unfortunately, the mechanism underlying HPV-TME interaction is still unclear. Our study aims to construct a new prognostic risk model in HNSCC. The respective relations between risk model, clinical characteristics, and OS were investigated based on HPVrelated genomic expression alteration. Cell type identification by estimating relative subset of known RNA transcripts (CIBERSORT) was applied to quantify the immune cell infiltration status in HNSCC TME. Finally, we explored the hub immune biomarkers to have an in-depth understanding of HNSCC immunotherapy.

Patients Data
The gene expression data and clinical information of 479 HNSCC patients were downloaded from The Cancer Genome Atlas (TCGA) data portal 1 . We further obtained 100 patients with known HPV status determined by p16 IHC staining and 83 patients with known HPV status determined by in situ hybridization (ISH). The differences between two methods were compared, and p16 testing was chosen for further analysis. The list of immune-related genes (IRGs) was achieved from the Immunology Database and Analysis Portal (ImmPort) database that provides more than 2,000 IRGs and annotations 2 .

Differential Gene Analysis and Functional Enrichment of DEGs
The differentially expressed genes (DEGs) in HPV-positive and HPV-negative HNSCC tissues were calculated using the limma package of R software 3 . mRNAs with an adjusted P < 0.05 and | log2 (fold change)| >1 were figured out as DEGs. Immune-related DEGs were identified by venn package. Heat maps of immune-related DEGs were drawn using the pheatmap package of the R software. ClusterProfiler package was applied to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of immune-related DEGs.

Identification of Prognosis-Related DEGs
Log-rank Kaplan-Meier survival analysis was performed to seek out prognosis-related genes from aforementioned immunerelated DEGs.

Construction and Validation of the Prognostic Risk Model
One hundred HNSCC patients with clear p16 status and complete clinical information were used for identifying hub prognostic immune-related signature and constructing prognostic risk model. Twelve prognostic-related DEGs were analyzed by a stepwise multivariate Cox proportional hazards regression analysis. Six genes were extracted, and the risk score was established with the following formula: risk score = expression of gene 1× coefficient 1+ expression of gene 2× coefficient 2+ . . . expression of gene n× coefficient n. We calculated the area under the curve (AUC) by using the survivalROC package of R software to validate the predictive ability of the prognostic risk model. Patients were stratified into the high-and low-risk groups with a threshold of the median risk score, and Kaplan-Meier survival analysis was performed to estimate the survival difference. A heat map was also created. To further validate and develop the prognostic risk model, we explored it in all 479 HNSCC patients. Relationship between risk score and clinical characteristics in HNSCC patients was also investigated and visualized by ggplot2 package.

Tumor-Infiltrating Immune Cells Fraction Calculation
We calculated tumor-infiltrating immune-cell fraction of each HNSCC patient by CIBERSORT algorithm (Newman et al., 2015). The reference 22 leukocyte expression data were downloaded from the website. Differences of immune-cell fraction between the high-and lowrisk groups were visualized by pheatmap package and vioplot package. Single sample Gene Set Enrichment Analysis (GSEA) was applied to validate the differences of immune-cell fraction by GSVA and GSEAbase packages.

Estimation of the Immunoreactivity
As hub immune response biomarkers, the expressions of PD-1, PD-L1, PD-L2, CTLA4, CD4, CD8A, and CD8B between the high-and low-risk groups were analyzed using Wilcoxon test. Based on machine learning and tumor genotype, immunophenoscore (IPS), is calculated and normalized with a range of 0 to 10. Higher scores are correlated with higher immunoreactivity (Charoentong

Statistical Analysis
A χ 2 test or Fisher exact test was used for calculating the differences of clinical characteristics. Wilcoxon rank sum 4 https://tcia.at/home test was used to estimate the statistical significance of continuous variables. For constructing the prognostic risk model, univariate Cox regression analysis and multivariate Cox regression analysis were used. Survival analysis was performed by survival package (Holleczek and Brenner, 2013). All the aforementioned statistics were analyzed by R software (version 3.6.3). A two-sided P < 0.05 was considered to be statistically significant.  Frontiers in Genetics | www.frontiersin.org

Clinical Characteristics of HNSCC Patients
We downloaded the clinical information of 479 HNSCC patients from TCGA data portal. In this study, we compared the results of HPV status identified by ISH or p16 testing, and there were few discrepancies (Supplementary Figure 1 and Supplementary  Table 1). Thus, we chose p16 status, a widely used marker in clinic, to identify HPV-driven HNSCC tumor. We identified 30 p16-positive and 70 p16-negative HNSCC patients (Figure 1 and Table 1). Compared with p16-negative patients, p16-positive patients had more early T-stage diseases and fewer death events, albeit the latter did not meet a statistical difference.

Identification of DEGs and Immune-Related DEGs
Gene expression profiles of p16-positive and p16-negative patients were obtained from the TCGA database. With a filter criterion of adjusted P < 0.05 and |log2 (fold change)| >1, 658 DEGs were identified by limma package (Figure 2A), of which 305 genes were up-regulated and 353 genes were downregulated. Seventy-three immune-related DEGs were extracted from the intersection of DEGs and IRGs downloaded from ImmPort database (Figure 2B), containing 28 up-regulated genes and 45 down-regulated genes (Figure 3). GO and KEGG enrichment analyses were conducted to uncover the allocated biological process, cellular component, molecular function, and pathway of immune-related DEGs. The identified immunerelated DEGs were found to be chiefly enriched in cytokine and receptor related functions and pathways. As shown in Figure 4A, the most significantly enriched terms were "cell chemotaxis, " "clathrin-coated vesicle membrane, " and "receptor ligand activity" under GO analysis. Top six enriched pathways were "cytokine-cytokine receptor interaction, " "Viral protein interaction with cytokine and cytokine receptor, " "EGFR tyrosine kinase inhibitor resistance, " "ErbB signaling pathway, " "MAPK signaling pathway, " and "chemokine signaling pathway" under KEGG analysis ( Figure 4B).

Identification of Prognosis-Related DEGs
Kaplan-Meier survival analysis was conducted to identify genes highly associated with OS of HNSCC patients. As shown in Figure 5, low expressions of DKK1, HBEGF, AREG, TNFRSF12, TGFA, RNASE7, PLAU, F2RL1, and IHNBA were significantly related with a favorable prognosis, whereas low expressions of CD79A, FAM3B, and PIK3R3 were linked to an ominous clinical outcome.

Construction of Immune-Related Risk Model
To explore the predictive value of aforementioned prognostic genes, we conducted a stepwise multivariate Cox regression analysis. Six of the 12 genes were extracted and considered significantly linked to the OS of 100 HNSCC patients. Forest plot showed DKK1, HBEGF, RNASE7, and TNFRSF12A were associated with poor outcomes, and INHBA and PIK3R3 were associated with favorable outcomes (Figure 6 and Table 2 According to the formula, we calculated the risk score of every HNSCC patient. We then divided the patients into lowrisk group (n = 50) and high-risk group (n = 50) by the median risk score. A receiver operating characteristic curve of 3-year OS (AUC = 0.815) was drawn to estimate the predictive ability of the risk model (Figure 7A). High-risk score group revealed a worse clinical outcome compared with low-risk group by Kaplan-Meier analysis (Figure 7B). The distributions of risk score and patients' survival status are displayed in Figures 7C,D. Heat map revealed gene expression differences between the high-and low-risk groups (Figure 7E).

Validation of the Risk Model in Total HNSCC Patients
To further validate the reliability of the constructed risk model, the risk score of whole 479 HNSCC patients were calculated and divided into either low-risk group (n = 187) or high-risk group (n = 292). The AUC of 3-year OS was 0.624 ( Figure 8A). Kaplan-Meier survival curves revealed  statistically significant difference between the high-and lowrisk groups (P < 0.001), and the high-risk group was associated with a worse clinical outcome ( Figure 8B). Based on the risk score, the distribution, survival status, and gene  expression differences were also evaluated and visualized (Figures 8C-E). Similar results were observed in whole 479 TCGA HNSCC patients, illustrating the reliable predictive ability of risk model.

Association Between the Risk Score and Clinical Characteristics and Their Prognostic Roles
The relationship between the risk model and clinical characteristics were analyzed by univariate and multivariate Cox regression analyses. Age older than 65 years, advanced N stage, and metastatic disease were significantly associated with unfavorable prognosis (Table 3). Risk score could serve as an independent predictive factor for HNSCC patients. Relationship between risk score and clinical characteristics were also analyzed (Figure 9), in which risk score and T stage seemed closely linked to each other.

Differential Immune Landscape in the High-and Low-Risk Groups
CIBERSORT algorithm was applied to demonstrate the relationship between risk score and tumor-infiltrating immunecell fractions. Samples with P < 0.05 were considered statistically different. Immune-cell fraction of 448 of 479 HNSCC patients   was calculated for further analyses ( Figure 10A). The abundance of immune cells between 268 high-risk and 180 low-risk patients were normalized and compared by Wilcoxon rank sum test ( Figure 10B). CD8 + T cells, Tregs, naive B cells, naive CD4 + T cells, and resting mast cells were evidently abundant in low-risk TME, whereas fewer resting NK cells, activated dendritic cells, activated mast cells, eosinophils, and neutrophils were present. A univariate Cox regression analysis-based forest plot displayed the association between immune cells and OS of HNSCC patients ( Figure 10C). Using ssGSEA analysis, we validated the differences of immune cells between the low-and high-risk groups, and similar results were observed (Supplementary Figure 2). We further investigated the expression of hub biomarkers of ICI response ( Figure 11A). The expressions of PD-1, CTLA4, PD-L2, CD4, CD8A, and CD8B were significantly higher in low-risk samples. Differences between two groups demonstrated that the specific changes of tumor-infiltrating immune cells might be associated with the OS of HNSCC patients. The scores of IPS and IPS with CTLA4 blockers, IPS with PD1/PDL1/PDL2 blockers, IPS with CTLA4, and PD1/PDL1/PDL2 blockers were calculated (Figures 11B-E), and IPS score in low-risk group was significantly higher compared with that in the high-risk group.

DISCUSSION
The incidence of HNSCC is continuously increasing around the world and may surpass that of cervical cancer (Marur and Forastiere, 2016). However, only a small proportion of HNSCC patients respond to and gain benefits from targeted drugs and ICI therapies (Concha-Benavente et al., 2016;Ferris et al., 2018). It is urgent and critical to find inventive ways to predict response to ICI treatment. Accumulating evidences have documented the potential effects of HPV on HNSCC patients in several aspects including genomic integration, carcinogenesis, tumor angiogenesis, and TME (Troy et al., 2013;Parfenov et al., 2014;Cao et al., 2019;Wookey et al., 2019). Therefore, HPV status of HNSCC patients provides a promising entry point for predicting response to ICI therapy and unearthing the elusive molecular mechanism. Within this study, we aim to further inspect how HPV infection impacts malignant gene expression profile and response of HNSCC to ICI therapy. In addition, it is of vital importance to develop a prognostic immune signature based on immune-related biomarkers, which might enhance the therapeutic efficacy of various HNSCC patients. It has been demonstrated that the p16 protein is overexpressed in the majority (82.2%) of HPV-associated (defined as HPV DNA positive) oropharyngeal carcinoma, and the detection of p16 by immunohistochemical staining is routinely used as a reliable surrogate marker in clinical practice and research (Li et al., 2004;El-Naggar and Westra, 2012;Ndiaye et al., 2014). Compared with other detection methods such as HPV DNA and E6/E7 mRNA detection by ISH that are usually expensive and time-consuming, p16 immunohistochemical staining is relatively inexpensive, timesaving, and convenient, so it is widely used for clinical detection and recommended by the modified eighth AJCC/UICC cancer staging system regarding HPV positivity (Lydiatt et al., 2017).
Previous studies pointed out that HPV infection and TME seem impossibly remote. HPV-positive status is a proficient effector in lymphocyte recruitment in the TME of cervical cancer (Wood et al., 2016). Along similar lines, in head and neck cancer, HPV-positive cell lines are competent in lymphocyte recruitment and cytokine secretion compared to HPV-negative cell lines (Stone et al., 2014). Moreover, HPV-positive tumor could polarize macrophages toward classically activated phenotype (M1) and augment antitumoral IL-6 secretion (Chen et al., 2019). In our study, it is highly likely that immune-related DEGs were allocated to cytokines and receptors-related functions and pathways, indicating that HPV infection might serve a dominant role in shaping the TME by altering the expression and cell recognition of cytokines, with the molecular mechanisms remaining uncertain.
Because of the complexity and heterogeneity of tumor biological behavior and immune response, it is unreliable to apply one single biomarker to fully illustrate and predict the prognosis and therapeutic response. Therefore, prediction model based on multibiomarkers may be a more effective and accurate tool, stratifying HNSCC patients into high-and low-risk subgroups and identifying those that might achieve clinical benefits from ICIs or other immunotherapies.
We developed and validated an immune-related prognostic risk model that gave an accurate prediction of survival outcome for HNSCC patients based on immune-related DEGs. Most of these genes hold a decisive place in expression of cytokine and its receptors. Overexpression of DKK1, an identified proinflammatory cytokine in quite a number of cancers (Chae and Bothwell, 2019), could boost malignant cell proliferation, migration, and invasiveness and give rise to unfavorable clinical outcomes in HNSCC (Shi et al., 2014;Gao et al., 2018). HBEGF serves as a crucial component of EGFR, mediating tumor cell proliferation and bringing about cetuximab resistance (Hatakeyama et al., 2010;Huang et al., 2014). Therefore, HBEGF expression contributes to locating the potential beneficiaries of cetuximab treatment. INHBA encodes proteins of TGFβ superfamily members, which is considered as hub gene in lymphatic metastasis and predictor of unsatisfying OS in HNSCC patients (Kelner et al., 2015;Chang et al., 2016). Moreover, in terms of PIK3R3, RNASE7, and TNFRSF12A, similar results have been obtained in other researches, indicating that these genes might be engaged in immune-related pathways and promote tumor progression (Scola et al., 2012;Wang et al., 2015;Eichler et al., 2016;Yang et al., 2018;Sun and Feng, 2020). In this study, expression patterns of these genes in HNSCC patients were suggestive to be capable of modifying the TME, adjusting immune response, initiating tumorigenesis, and predicting prognosis. Given the high AUC and statistically significant differences of OS between two groups, immunerelated risk model functioned as a strong indicator of OS in HNSCC patients. Furthermore, the risk score, correlated with T stage, might act as an independent prognostic factor.
Composition of tumor-infiltrating immune cells was estimated by CIBERSORT and ssGSEA for each sample. Abundance of CD8 + T cells and Tregs, was richer in the low-risk group and was correlated with favorable prognosis. In this study, patients with higher CD8 + T cell infiltration had a better clinical outcome, which was consistent with previous investigations (Balermpas et al., 2014;de Ruiter et al., 2017). It is a well-acknowledged fact that Tregs work as an adverse factor in antitumor immunity. Quite the reverse, recent studies demonstrated that compared with HPV-negative HNSCC patients, HPV-positive patients have a higher level of Treg infiltration (Mandal et al., 2016) and usually a longer survival time (Mirghani et al., 2019). Likewise, higher infiltration of Tregs and expression of immune checkpoints are revealed in HNSCC microenvironment, indicating a complicated relationship between Tregs and the antitumor immune response (Mandal et al., 2016;Saloura et al., 2019). High-degree infiltration of Tregs and its correlation with better prognosis observed in this study were possibly due to the interaction between Tregs and other immune components in the TME.
The hub biomarkers of ICI were explored. Low-risk patients displayed up-regulation of PD-1 and CTLA4. Meanwhile, in the risk model, patients in the low-risk group had a remarkably increased IPS score. Aleix and colleagues described that PD-1 gene expression is associated with longer progression-free survival (Prat et al., 2017). Furthermore, anti-CTLA4 therapy could enhance antitumor effect via Treg exhaustion (Selby et al., 2013;Simpson et al., 2013). Hence, high infiltration of Tregs in the low-risk group may be conducive to combination of ICIs and depletion strategies targeting Tregs. Collectively, these results demonstrated that the TME of low-risk group had a higher immunoreactivity. The risk model might hold the capacity to determine tumor response to immunotherapy in HNSCC patients. IPS with CTLA4 and/or PD1/PDL1/PDL2 blocker scores were of no statistical significance, revealing that the prediction value of the risk model in response to immunotherapy was not fully understood and required further research.
In conclusion, we constructed a gene model based on HPV status, providing a possible method to estimate the TME and predict OS and response to immunotherapy in HNSCC patients. Despite the aforementioned promising results, there still exist some limitations. First, the datasets are based on a public database. Second, the number of HPV-positive HNSCC patients is relatively small. Whether a larger amount of cases enrolled would affect the study results remains unknown. In addition, we focused on the DEGs within p16-positive and -negative groups, which is only one of the potentially useful stratifying strategies in head and neck cancers. Further, we focused only on the biological function and pathway enrichment of immune-related DEGs due to limited research scope, and many potential DEGs failed to be investigated. However, the promising results carry more weight over these limitations.

CONCLUSION
This study has examined the differential expression pattern between the HPV-positive group and HPV-negative group defined by p16 status, and 12 DEGs were perceived as prognosisrelated genes. A six-immune-related DEGs-based risk model was constructed, and the predictive ability, TME, and ICI response were consequently estimated. The risk model is reliable and provides a possible method for clinical outcome prediction and immunotherapy improvement.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: https://portal.gdc.cancer.gov/.

ETHICS STATEMENT
Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

AUTHOR CONTRIBUTIONS
XL: substantial contributions to conception, design, and data analysis. JC: bioinformatics algorithms and analysis instruction. WL: datasets acquisition and data preprocessing. ZZ and JL: construction of stepwise multivariate Cox analysis and validation. YPG, XJ, and QW: drafting the article. YG and QW: supervision of the study and critical revision of manuscript. CX: conception of the study, and design the workflow. All authors contributed to manuscript revision, read and approved the submitted version.