Construction of the prognostic signature of alternative splicing revealed the prognostic predictor and immune microenvironment in head and neck squamous cell carcinoma

Background: Head and neck squamous cell carcinoma (HNSC) is a prevalent and heterogeneous malignancy with poor prognosis and high mortality rates. There is significant evidence of alternative splicing (AS) contributing to tumor development, suggesting its potential in predicting prognosis and therapeutic efficacy. This study aims to establish an AS-based prognostic signature in HNSC patients.

Methods: The expression profiles and clinical information of 486 HNSC patients were downloaded from the TCGA database, and the AS data were downloaded from the TCGA SpliceSeq database. The survival-associated AS events were identified by conducting a Cox regression analysis and utilized to develop a prognostic signature by fitting into a LASSO-regularized Cox regression model. Survival analysis, univariate and multivariate Cox regression analysis, and receiver operating characteristic (ROC) curve analysis were performed to evaluate the signature and an independent cohort was used for validation. The immune cell function and infiltration were analyzed by CIBERSORT and the ssGSEA algorithm.

Results: Univariate Cox regression analysis identified 2726 survival-associated AS events from 1714 genes. The correlation network reported DDX39B, PRPF39, and ARGLU1 as key splicing factors (SF) regulating these AS events. Eight survival-associated AS events were selected and validated by LASSO regression to develop a prognostic signature. It was confirmed that this signature could predict HNSC outcomes independent of other variables via multivariate Cox regression analysis. The risk score AUC was more than 0.75 for 3 years, highlighting the signature’s prediction capability. Immune infiltration analysis reported different immune cell distributions between the two risk groups. The immune cell content was higher in the high-risk group than in the low-risk group. The correlation analysis revealed a significant correlation between risk score, immune cell subsets, and immune checkpoint expression.

Conclusion: The prognostic signature developed from survival-associated AS events could predict the prognosis of HNSC patients and their clinical response to immunotherapy. However, this signature requires further research and validation in larger cohort studies.

Introduction

Head and neck squamous cell carcinoma (HNSC) is the sixth most common cancer globally, with an estimated 900,000 new cases and 450,000 deaths annually (Ferlay et al., 2019). HNSC patients are normally diagnosed at the later stages with a poor prognosis. Although chemotherapy, radiation, and targeted therapies have made great progress, the survival rate of advanced or recurrent HNSC remains low (Cohen et al., 2016). This is due to the heterogeneity of HNSC. Patients with advanced HNSC are treated with cetuximab, an anti-EGFR antibody, with a 13% success rate (Licitra et al., 2011). An immunotherapy agent, anti-PD1 antibody, successfully stimulated anti-tumor immunity and produced a significant clinical response in patients with aggressive HNSC. However, only a subset of patients (∼18%) benefitted from this strategy, while most HNSC patients displayed clinical resistance (Chow et al., 2016; Ferris et al., 2016; Seiwert et al., 2016). Therefore, it is crucial to identify new markers to provide an accurate HNSC prognosis and the appropriate treatment strategies.

Alternative splicing (AS) is a process in which mRNA precursors are alternatively spliced and ligated to produce mature mRNAs for protein diversity (Larochelle, 2016; Ule and Blencowe, 2019). Research has identified the contribution of AS to a variety of diseases, including cancer. The dysregulated expression of splice isoforms is considered a potential driver of tumor development and progression (Sciarrillo et al., 2020). Recent studies suggest that cancer-specific splice isoforms can be important signatures for predicting treatment efficacy. For instance, squamous cell carcinoma patients with EGFR isoform D splicing patterns showed better responses toward EGFR-TKIs (Tan et al., 2017). More importantly, AS was reported to contribute to the development of the immune microenvironment (Frankiw et al., 2019; Li et al., 2019; Yu et al., 2020). Changes in AS can alter both immune cell infiltration and tumor-associated immune cytolytic activity. Hence, AS signature may accurately predict the clinical response of HNSC patients to immunotherapy.

Herein, a comprehensive analysis of genome-wide AS events was performed with RNA sequencing (RNA-Seq) data obtained from TCGA-HNSC samples. We developed a novel prognostic signature from the data via LASSO regression analysis. HNSC patients could be categorized into either high or low-risk groups based on their risk scores calculated with the prognostic signature. Immune infiltration and functional analysis were then performed to identify the role of the signature in the tumor microenvironment. Our analysis suggested that this prognostic signature could be an effective tool for predicting the prognosis of HNSC patients and their clinical response to immunotherapy.

Methods

Data collection and pretreatment

The expression profiles, somatic mutation, and matching clinical follow-up information of HNSC patients were obtained from the TCGA database. AS data was obtained from the TCGA SpliceSeq database. A total of 486 patients were enrolled based on the following criteria (Ferlay et al., 2019): histologically confirmed HNSC (Cohen et al., 2016); patients with RNA-Seq data; and (Licitra et al., 2011) patients with basic clinical and follow-up information. Patients who lacked AS data in the TCGA SpliceSeq database had been excluded. HPV status information was collected from published data and the TCGA database (Cancer Genome Atlas Network, 2015). The clinical characteristics of patients are listed in Supplementary Table S1. For validation, RNA-Seq data of an independent cohort (86 oral cavity squamous cell carcinoma patients) was accessed from the European Nucleotide Archive (study accession: PRJEB24758) and processed by SpliceSeq software to obtain the AS profiles. The percent spliced-in index (PSI), a ratio of normalized read counts indicating the inclusion of a transcript element over the total normalized reads, was used to quantify AS events (Ryan et al., 2012). From here, seven different types of AS events were identified: exon skip (ES), alternate donor site (AD), alternate acceptor site (AA), retained intron (RI), exclusive exons (ME), alternate terminator (AT), and alternate promoter (AP). To obtain a set of more reliable AS events, several conditions were implemented (i.e., percentage of samples with PSI value ≥75, average PSI value ≥0.05). Figure 1 depicts the flow chart representative of this study.

FIGURE 1

FIGURE 1. The flow chart of AS-based prognostic signature construction and evaluation.

Identification and functional analysis of survival-associated AS events

The hazard ratios (HRs) and 95% confidence interval (95% CI) of the overall survival of AS events were calculated after performing a univariate Cox regression analysis. AS events with p < 0.05 were considered survival-associated AS events. The seven types of identified AS events were utilized to develop an UpSet plot with a “UpSetR” package (version 1.4.0). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed to determine the parental genes in survival associated AS events. To visualize the results, functionally organized networks were constructed through the Cytoscape plug-in, ClueGO (version 2.5.8), with network connectivity (κ-score) ≥ 0.5, as described by Bindea et al. (Bindea et al., 2009).

Splicing correlation network construction

The list of 71 experimentally validated splicing factor (SF) genes was obtained from the SpliceAid 2 database, and their expression profiles were obtained from the TCGA database. Pearson correlation analysis was performed between SF expression and PSI value of survival-associated AS events, with several conditions (i.e., |R| ≥ 0.7, and p < 0.001). The splicing correlation network was then established using the Cytoscape software (version 3.8.2).

Construction and evaluation of the prognostic signature

Our prognostic signature was developed based on a LASSO-regularized Cox regression model (Luo et al., 2020; Wang et al., 2021a; Guo et al., 2021). The survival-associated AS events obtained above were input into the LASSO regression analysis with a “glmnet” package (version 4.1.3) to remove highly correlated variables and prevent overfitting. The optimum shrinkage parameter (λ) was determined using 10-fold cross-validation. AS events with non-zero parameters estimated in the model were selected and used to fit a multivariate Cox regression model to obtain the regression coefficients. The risk score of each sample was calculated by the formula: $RiskScore={\displaystyle \sum }_{i=1}^n{Coef}_i\times {PSI}_i$, where ${Coef}_i$ is the coefficient and ${PSI}_i$ is the PSI value of selected AS events.

From the median risk score, patients were classified into two risk groups. Kaplan–Meier survival analysis with a log-rank test was conducted to verify the statistical differences. Then, each selected AS event’s risk score distribution, survival time, and expression thermogram were visualized. Univariate and multivariate Cox regression models were fitted to analyze whether the risk scores were independent of clinical variables. A time-dependent receiver operating characteristic (ROC) curve analysis was performed to evaluate the prediction capability of the prognostic signature.

Validation of the prognostic signature was performed with an indirect method due to insufficient RNA-Seq data with survival information. Briefly, three published and well-validated prognostic signatures served as references and were used to calculate risk scores for the TCGA-HNSC and independent cohort samples (Wang et al., 2021c; Liu et al., 2021a; Chen et al., 2022). Our AS-based signature was also used to score the patient with the PSI profiles produced by SpliceSeq software. Correlation analysis and heatmaps were then used to demonstrate the consistency between the calculated risk scores of the prognostic signature and other signatures.

Immune cell infiltration and functional analysis

We used ESTIMATE, CIBERSORT, and ssGSEA to measure immune cell infiltration and functional differences between the risk subgroups. The presence of infiltrating stromal and immune cells in tumor samples was estimated with the “estimate” package (version 1.0.13) and classified into four scores: Stromal, Immune, ESTIMATE, and Tumor purity. The ESTIMATE score that infers tumor purity is the sum of immune and stromal enrichment scores and can be converted to tumor purity score as previously mentioned (Yoshihara et al., 2013). Immune cell infiltration in each sample was quantified by CIBERSORT, a validated deconvolution algorithm for characterizing cell composition based on the leukocyte signature matrix (LM22) (Newman et al., 2015). We used an R script to implement this algorithm with 1000 permutations, without quantile normalization. The quantified immune population was validated with ssGSEA using the “GSVA” package (version 1.40.1). The immune signatures that were identified using this approach included immune cell types (e.g., plasmacytoid DCs, inactivated DCs, activated DC, DCs, NK cells, B cells, mast cells, neutrophils, macrophages, CD8⁺ T cells, helper T-cells, Tfh, Th1, Th2, and Treg), immune-related functions (e.g., APC co-inhibition, APC co-stimulation, T cell co-inhibition, T cell co-stimulation, cytolytic activity, type I IFN response, type II IFN response, pro-inflammation, HLA, MHC class I, and CCR), and immune-related pathways (e.g. extracellular matrix, epithelial-mesenchymal transition, angiogenesis, and VEGF signaling pathway). These signatures were collected from publications and the MsigDB database (Barbie et al., 2009; Xue et al., 2019). A heatmap and a boxplot were plotted to visualize these results with the “heatmap” (version 1.0.12) and “ggplot2” (version 3.3.5) packages. The differential expression analysis for the identification of immune checkpoints between the two groups was performed with the “limma” package (version 3.48.3) and the Wilcoxon test. The correlation coefficient between risk score and immune checkpoint expression was calculated with Spearman correlation analysis. The tumor mutational burden (TMB) in each HNSC patient was calculated with the “maftools” package (version 1.0.12) from the somatic mutation data downloaded from the TCGA.

Functional analysis of the AS parent genes involved in the prognostic signature

We performed expression, survival, and functional analysis to determine the AS parent genes. For expression analysis, the log2-transformed Fragments Per Kilobase Million (FPKM) value of each gene was counted, and the differences between tumor and normal groups were analyzed by the GraphPad 8.0 software and the Mann-Whitney test. For survival analysis, HNSC patients were categorized into two groups based on their gene expression. The survival outcomes were determined through a Kaplan-Meier analysis. The pathways that were significantly linked to the expression of parent genes were identified with GSEA analysis (version 4.1.0) of the Hallmark gene sets. Using the gene expression lists ranked by Pearson, we calculated the weighted enrichment scores. Gene sets with a nominal p-value < 0.05 and FDR ≤0.1 were considered statistically significant. The first seven significantly correlated pathways were visualized using GraphPad 8.0 software. For immune function analysis, Spearman’s correlation method was used to calculate the correlation coefficient between gene expression and immune checkpoint or immune cell infiltration. The immune correlation network was then established using the Cytoscape software (version 3.8.2) with a p-value < 0.05 threshold.

Statistical analyses

R (version 4.1.0) and GraphPad 8.0 were utilized for all statistical analyses. p-value < 0.05 was statistically significant. The categorical clinical characteristics variables were analyzed by the Chi-square and Fisher’s exact tests, while the continuous variables were analyzed by the Mann-Whitney non-parametric and ANOVA tests. The non-normal data distribution was analyzed via Spearman correlation analysis, while the continuous variables with normal distribution were analyzed via the Pearson correlation analysis.

Result

Identification of survival-associated AS events in HNSC

Univariate Cox regression analysis were done on 486 HNSC patients to evaluate the correlation between AS events and OS. Out of 1714 genes, 2726 AS events were identified as survival-associated. ES was the predominant type (26.1%), followed by AP (23.5%) and AT (23.2%). Given the multiple splicing patterns for a single gene, the genes intersecting sets and survival-associated AS events were visualized using the UpSet plot (Figure 2A). Function enrichment analyses were performed to study the biological characteristics of AS events. Our findings revealed that the parent genes were enriched in KEGG and BP, which were related to important biological functions, including mRNA surveillance pathway, adherens junction, RNA transport, protein localization, intracellular transport, and mRNA metabolic process (Figure 2B).

FIGURE 2

FIGURE 2. Identification and functional analysis of survival-associated AS events in HNSC. (A) Interactive sets between the seven types of survival-associated AS events (n = 2726) shown in an UpSet plot. (B) Functionally grouped Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) categories for parent genes of survival-associated AS events. GO and KEGG categories were grouped based on their similarity (κ-score ≥ 0.5), and the most significant term in each group is shown in bold. (C) Correlation network of splicing factors and survival-associated AS events. This network was built based on significant correlations (|R| ≥ 0.7, p < 0.001) between the expression of 71 splicing factors and the PSI values of survival-associated AS events. Splicing factors are represented with orange triangles, and AS events are represented with circulars (red/blue represents favorable/inferior prognosis). The red/blue lines are represented with positive/negative correlation. Data were analyzed using Pearson’s correlation method.

Splicing factors (SFs) regulate AS events. To better understand the regulation of survival-associated AS events, we performed a Pearson correlation analysis between reported 71 expressions of SFs and PSI values. Following the significant correlation (|r|≥0.7, p < 0.001), a splicing regulatory network was built, which contained 206 SF-AS pairs, 21 SFs, and 88 survival-associated AS events (Figure 2C). Most SF-AS pairs were positively correlated in this network. Splicing factors such as DDX39B, PRPF39, and ARGLU1 might be key regulators in AS events.

The construction of the prognostic signature of AS events revealed the prognostic predictor in HNSC

Survival-associated AS events were fitted to LASSO Cox regression analysis to remove highly correlated variables and prevent overfitting when constructing the prognostic signature. The shrinkage parameter was determined using 10-fold cross-validation. This signature achieved minimal deviation with 11 OS-associated AS events, but only eight of them were retained after optimization by a stepwise multivariate Cox regression analysis (Figures 3A,B). The details of these AS events, including their corresponding coefficients and hazard ratios are listed in Supplementary Table S2. The risk scores of HNSC patients were calculated according to the PSI values and their coefficients as follows: $\mathit{riskscore}=-3.77\times AIG1\left|77971\right|AT-2.71\times PACS2\left|29633\right|AP-2.20\times PTGR1\left|87219\right|AA-2.46\times RHOT1\left|40176\right|ES-1.98\times AGTRAP\left|670\right|AA-0.98\times ABCC5\left|67820\right|RI-1.45\times SH3KBP1\left|88642\right|AP-1.02\times RBMX\left|90220\right|RI$.

FIGURE 3

FIGURE 3. Construction and evaluation of prognostic signature. (A) Distribution of LASSO regression coefficients for survival-associated AS events. (B) Shrinkage parameter selection in the LASSO model used ten fold cross-validation via minimum criteria. Vertical dotted lines were drawn at the optimum λ values. (C) The risk score distribution, overall survival status, and expression profile of selected AS events of each HNSC sample. (D) Kaplan-Meier survival curve of overall survival. Log-rank test was used for data analysis. (E,F) Forest plots of hazard ratios (HRs) calculated by univariate (E) and multivariate (F) Cox regression for risk score and clinical features associated with overall survival. (G) Receiver operating characteristic (ROC) curve and the area under the curve (AUC) of risk score and clinical features in predictive performance for HNSC patients. (H) ROC plots and AUC of the risk score at 1, 2, and 3 years.

To study the correlation between risk scores and clinical features, HNSC patients were categorized into two risk groups. The clinical characteristics of these groups are shown in Table 1. Notably, there were more HPV-positive patients in the low-risk group. The score distribution, overall survival status, and expression profile of AS events are as plotted in Figure 3C. The high-risk group had a lower OS rate as the Kaplan-Meier survival analysis indicated (p < 0.001, Figure 3D), which validated the prediction capability of our prognostic signature. To establish that risk score is independent of other variables in predicting HNSC outcome, univariate and multivariate Cox regression analyses were employed using clinical features like age, gender, pathologic stage, and HPV status. The analysis reported that the risk score variable was statistically significant (Figures 3E,F). From the ROC analysis, the AUC of risk score was higher than the other variables and remained above 0.75 for 3 years, which suggested a powerful predicting capability of our signature (Figures 3G,H). Together, these results revealed the satisfactory efficiency of our signature in predicting prognosis for HNSC patients.

TABLE 1

TABLE 1. Correlation between risk score based subgroups and clinical characteristics.

To validate our signature in an independent cohort, we utilized three published and well-validated prognostic signatures as references in our study. Correlation analysis reported that risk scores correlated significantly with the reference signatures in the TCGA-HNSC and independent cohorts (Figures 4A–F). The risk status of HNSC patients was visualized with a heatmap, revealing the consistency of our signature with the reference signatures in predicting the prognosis of HNSC (Figures 4G,H). These results validated the credibility of our prognostic signature.

FIGURE 4

FIGURE 4. Validation of our signature in a independent HNSC cohort. (A–F) Scatter plots of correlations between the risk scores calculated by our signature and that calculated by reference signatures in TCGA-HNSC and independent cohort. Data were analyzed using Spearman’s correlation method. (G,H) Heatmaps showing the risk status of each sample predicted by our signature and reference signatures.

Correlation of risk score with the immune microenvironment signatures

Given the pivotal role of the immune microenvironment in tumor development and progression, we investigated the correlation between risk score and immune features. As shown in Figures 5A–D, the immune and stromal cells in tumors were quantified based on the ESTIMATE algorithm. This finding indicated that the high-risk group had lower immune and stromal cell concentration and higher tumor purity (Figures 5A–D). Immune subpopulation analysis demonstrated that tumors in the high-risk group were highly infiltrated with immune cells such as M2 and M0 macrophages, and lowly infiltrated with CD8⁺ T cells, Treg, memory CD4 T cells, naïve B cells, and M1 macrophages (Figure 5E). Further correlation analysis indicated that the infiltration of activated mast cells, eosinophils, M2 macrophages, and M0 macrophages had a positive correlation with the risk score. On the other hand, naïve B cells, plasma cells, T follicular helper cells, Tregs, resting mast cells, CB8+ T cells, memory CD4 T cells, and γδ T cells had a negative correlation with the risk score (Supplementary Figure S1A, ranked by the Spearman correlation coefficients). Related immune function signatures were quantified using ssGSEA analysis and visualized with a heatmap (Figures 5F,G). The result was consistent with the distribution of immune cell subpopulations, where the immune functions such as cytolytic activity, T cell co-stimulation, and APC co-stimulation were lower in the high-risk group. We also analyzed the expression of immune checkpoints and demonstrated that most immune checkpoints were lowly expressed in the high-risk group (Supplementary Figure S1B). The correlation analysis indicated that the risk scores were negatively correlated with PDCD1, CD274, CTLA4, and CD80, and positively correlated with CD276 (Figures 5H–L). TMB was linked to immune infiltration of tumors (Goodman et al., 2017; Jardim et al., 2021), but no significant difference was found between the two groups (Supplementary Figure S1C).

FIGURE 5

FIGURE 5. Immune cells infiltration and functional analysis. (A–D) The analysis of the stromal score, immune score, ESTIMATE score, and tumor purity between the high- and low-risk groups. (E) Immune cell subpopulations analysis between the high- and low-risk groups based on CIBERSORT algorithm. (F) Heatmap of the enrichment score of immune cells subpopulations and immune-related pathways signatures based on ssGSEA algorithm. (G) Differential enrichment analysis of immune cells subpopulations and immune-related pathways between the high- and low-risk groups based on ssGSEA algorithm. (H–L) Scatter plot of correlations between representative immune checkpoint expression and risk score. Data were analyzed using Spearman’s correlation method. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.

AS events involved in our signature might play a role in the HNSC tumorigenesis, prognosis, and immune regulation.

We studied the expression, function, and prognosis of parent genes in AS events. Expression analysis indicated that five genes had significant differential expressions between tumor and normal samples (Figure 6A). AGTRAP1, ABCC5, SH3KBP1, and RBMX were upregulated, while PTGR1 was downregulated in tumor samples as compared to normal samples (Figure 6A). Survival analysis showed that AIG1, PACS2, PTGR1, AGTRAP, SH3KBP1, and RBMX expressions significantly affected the overall survival (Figures 6B–I). High expressions of AGTRAP, SH3KBP1, and RBMX were found in tumor tissues and correlated with a poor prognosis. This finding suggested that these genes could have a role in tumor initiation and progression. The lack of correlation between the PTGR1 expression and prognosis suggested its complex regulatory mechanism. We then analyzed the pathways and functions of these genes using GSEA with Hallmark gene sets. The seven significantly enriched gene sets are shown in Figures 6J–M. These results revealed that these parent genes could be involved in HNSC development. Notably, SH3KBP1 had a significantly positive correlation with immune-related pathways, such as complement system, allograft rejection, and IL-2-STAT5 signaling (Figure 6K). In addition, PTGR1 had a positive correlation with the metabolic pathways and a negative correlation with the immune pathways (Figure 6M).

FIGURE 6

FIGURE 6. Expression and functional analysis of the parent genes of AS involved in prognostic signature. (A) Heatmap and differential expression analysis of the parent genes between HNSC and normal samples. The filled color represents the log2 (FPKM value) for each gene. Data were analyzed using Mann-Whitney nonparametric test. (B–I) Kaplan-Meier survival curve of overall survival for each parent gene. Log-rank test was used for data analysis. (J–M) GSEA plots for the top seven Hallmark gene sets significantly correlated with the parent gene expression. The enrichment scores were calculated with gene expression lists ranked by pearson. NES, normalized enrichment score. (N) The network diagram shows the correlation between parent genes, immune checkpoints, and immune cell subsets (CIBERSORT). Colored triangles, squares, and circles represent respectively parent genes, immune checkpoints, and immune cells, respectively. The thickness of the line represents the strength of correlation. Red and blue represent positive and negative correlation, respectively. Data were analyzed using Spearman’s correlation method.

We also analyzed the correlation of these genes with immune cells and the expression of immune checkpoints. The analysis is illustrated in Figure 6N. It was worth noting that SH3KBP1 was positively correlated with immune checkpoints such as CD80, CD86, HAVCR2, PDCD1LG2, CTLA4, PDCD1 and CD274 (r = 0.55, 0.53, 0.52, 0.47, 0.47, 0.44, and 0.38, respectively). These results highlighted its role in the progression and prognosis of HNSC. Furthermore, there was a negative correlation between ABCC5 and PDCD1LG2 (r = -0.36), between AGTRAP and B and plasma cells (r = −0.28, −0.22), and between PACS2 and memory CD4 T cells and neutrophils (r = −0.23, 0.21). These genes could play a role in carcinogenesis and the tumor immune microenvironment.

Discussion

HNSC is a diverse group of cancers occurring in the oral cavity, oropharynx, and larynx. The limitations of the anatomic region, TNM stage, and HPV status as predictors for prognosis and treatment outcomes have limited clinical applications. Genome-wide research revealed that HNSC patients can be categorized into separate AS subgroups according to their AS patterns, which exhibited uneven distribution of survival status, EGFR mutation/amplification, TP53 mutation, and immune characteristics (Li et al., 2019). These data indicated that AS events may help stratify high-risk patients and predict the treatment response. In this study, we used the LASSO-regularized Cox regression model to screen eight prognostic AS events and established a prognostic signature for HNSC. Patients can be easily classified as high- or low-risk, according to our signature. Our results reported that the high-risk group had a significantly worse survival outcome. The multivariate Cox regression analysis revealed that risk score was independent of clinical variables such as age, pathological stage, and HPV status. The ROC curve displayed the AUC of risk score above 0.75 in both the short- and long-term. These results highlighted the predictive capability of our signature.

Previously, a genome-wide AS profiling analysis revealed that clusters of AS with distinct patterns correlated with different immune statuses (Li et al., 2019). This was consistent with our finding that the AS risk score was connected to distinct immune cell populations and immune activation. The findings indicated that the predicted low-risk group had a higher level of immune cell infiltration (e.g., CD8⁺ T cells, Tregs, memory CD4 T cells, naïve B cells, and M1 macrophages) and immune activation signatures (e.g., cytolytic activity, CD8+/Treg ratio, and IFN-γ signaling). Conversely, most immune cells and activation signatures were lowly expressed in high-risk patients. However, we observed an increased infiltration of M2 and M0 macrophages in these patients. These results characterized low-risk patients as “Immune Class”, with high immune cell infiltration, enhanced cytolytic activity, and active IFN signaling (Mandal et al., 2016; Chen et al., 2019). In contrast, high-risk patients were labeled as “non-Immune Class”, with “Exhausted Subtype” characteristics such as increased M2 macrophage infiltration (Chen et al., 2019).

Anti-tumor immunity is dependent on several key aspects:1) recognition of tumor-specific antigens; 2) immune cell infiltration; 3) freedom from immunosuppressive effects of immunoregulatory pathways and cells (Kim and Chen, 2016). Non-synonymous mutation load, which can result in the expression of neoantigens, may contribute to immunogenicity and inflamed cancer phenotype (Goodman et al., 2017; Jardim et al., 2021). However, we did not find a difference in the mutation load between high- and low-risk groups. In addition to the neoantigens generated by mutations, viral antigens can compromise the immune system and are targeted by T-cells (Tashiro and Brenner, 2017). HPV + HNSC has a better prognosis than HPV-cancer, reflecting the high level of viral activity (Mandal et al., 2016; Gameiro et al., 2018; Johnson et al., 2020). Our research found that the proportion of HPV + HNSC was lower in the high-risk group, indicating that the difference in the HPV + cancer distribution may be part of the low immunity in the high-risk group. Antigens released by tumor cells are taken up by antigen-presenting cells (APC), such as dendritic cells (DCs), which become activated and migrate to the tumor-draining lymph node. The activation of DCs is known to initiate the cancer-immunity cycle (Steinman, 2012). In comparison with the low-risk group, DC signatures (e.g., activated and plasmacytoid DCs) were lowly expressed in the high-risk group, contributing to the low immune infiltration in these patients. For successful trafficking or infiltration, the immune cells need to overcome the pressure from the tumor itself and the surrounding stroma cell (Ariffin et al., 2014; Joyce and Fearon, 2015; Turley et al., 2015). Stromal cells can reduce lymphocyte adhesion by regulating the expression of adhesion molecules (e.g., ICAM-1 and VCAM-1) (Griffioen et al., 1996; Bouzin et al., 2007). Recent studies have shown that activated stromal cells combined with increased TGF-β and M2 infiltration can lead to an immunosuppressive phenotype. This correlated to poor prognosis and good response to immunotherapy (Chen et al., 2019). Nevertheless, no difference was found in stroma cell scores and extracellular matrix signatures between the two risk groups. This indicated that the exclusion exerted by stromal cells may not be the reason for the difference in immune infiltration between the two groups.

Anti-tumor immunity also relieves the immunosuppressive effects from immune pathways and regulatory cells. Our research demonstrated that in the low-risk group, most immune checkpoints were highly expressed, and some were negatively correlated with the risk score. These results are consistent with a previous observation that immune-high HNSC is often accompanied by a high immunoregulatory response (Mandal et al., 2016). Interestingly, we found a high expression of CD276 in the high-risk group. CD276, also known as B7 homolog three protein (B7-H3), is recently considered a T cell inhibitor that promotes tumor proliferation and invasion (Zhou and Jin, 2021). Compared to other immune checkpoints, CD276 affects immunity, as well as regulates the cancer cell aggressiveness via multiple non-immune pathways (Liu et al., 2011; Li et al., 2017; Shi et al., 2019). Overexpression of CD276 was found in various tumors, including HNSC, and conferred a poor prognosis (Liu et al., 2021b). Recent research reported a negative correlation between CD276 expression and CD8⁺ T cell infiltration in human HNSC (Mao et al., 2017). Inhibition of CD276 increased the infiltration of CD8⁺ T cells and T cell activation in mice models (Mao et al., 2017; Wang et al., 2021b). Given these results, the abnormal expression of CD276 might be responsible for excluding the immune cells in the high-risk group. Immune cells are an important source of immunosuppression. Our results found that M2 macrophage infiltration was increased in high-risk patients. M2 macrophages have been demonstrated to reduce tumor-infiltrating lymphocytes, especially CD8⁺ T cells, by reducing the expression of chemokines and promoting the production of extracellular matrix (Peranzoni et al., 2018; Pan et al., 2020). The depletion of M2 macrophages had been reported to restore CD8⁺ T cell migration and infiltration, and improve the efficacy of immunotherapies (Peranzoni et al., 2018). Therefore, high expression of M2 macrophages may led to low immune infiltration in the high-risk group.

Immune hot tumors are better recognized by the immune system and can trigger a better response to checkpoint therapies (Kim and Chen, 2016). Therefore, we hypothesized that checkpoint blockade may benefit low-risk HNSC patients, with high immune infiltration and immune checkpoint expression (e.g., PD-1 and CTLA-4). Tregs are also highly infiltrated in these patients, indicating that a combination of molecular antagonists (e.g., CTLA-4, CCR4, and STAT3 antagonist) can attain a better response rate (Selby et al., 2013; Sugiyama et al., 2013; Woods et al., 2018). In contrast, the immune landscape of the high-risk group is characterized by lower immune cell infiltration and higher expressions of CD276 and M2 macrophages. This kind of immunologically cold tumor can be responsive to the checkpoint blockade by a combination of therapies that promote tumor immune cell infiltration and convert tumors into an immunologically hot phenotype (Ribas et al., 2018). Therefore, we hypothesized that the strategy of combining the following therapies can be more effective in the high-risk group:

1) Cytokine and tumor vaccine therapies. A randomized phase III trial of IL-2 given via the perilymph after oral cavity surgery reported that cancer patients exhibited >25% improvement in OS (De Stefani et al., 2002). Furthermore, perilymphatic delivery of cytokines increased lymphocyte infiltration and improved prognosis in HNSC patients (Berinstein et al., 2012). Cancer vaccines can bypass the immune cold tumor microenvironment and deliver antigens directly to the APC (Tan et al., 2018; Shibata et al., 2021). This approach displayed significant potential to upregulate CD8⁺ T cells and sustain their function (Ott et al., 2017). HPV E6 and E7 are ideal vaccine targets for HNSC because of their highly immunogenic epitopes. Clinical trials assessing E6/E7 protein-based vaccines, DNA-based vaccines encoding E6/E7, or pathogen vector-based vaccines containing E6/E7-encoding DNA are currently in progress (Tan et al., 2018; Shibata et al., 2021).

2) Immunotherapy based on targeting CD276. CD276 expression is associated with resistance to anti-PD-1 immunotherapy in non-small cell lung and ovarian cancer (Yonesaka et al., 2018; Cai et al., 2020). In solid tumor mice models, blocking of CD276 with mAbs increased CD8⁺ T and NK cell infiltration, decreased tumor development, and prolonged survivability (Wang et al., 2021b). In addition to mAbs, chimeric antigen receptor (CAR) T cell technology is another option for targeted CD276 immunotherapy. Autologous T cells are designed with CARs to target tumor antigens and destroy cancer cells. Majzner et al. described a CAR-T cell system targeting CD276 which showed significant effectiveness against various xenograft cancer types, including osteosarcoma, medulloblastoma, and Ewing sarcoma (Majzner et al., 2019). They also demonstrated that the efficacy of CAR-T cells depended on the density of CD276 on the tumor surface. Recently, Tang et al. demonstrated that the local administration of CAR-T cells limited tumor development without off-target toxicity or major adverse effects in recurrent anaplastic meningioma (Tang et al., 2020). This evidence supports the application of CD276 CAR-T cell therapy in NHSC patients. Hence, future research should focus on this aspect.

3) Therapy to inhibit M2 macrophages. M2 macrophages, derived from peripheral blood monocytes, are recruited at the tumor site via the CCL2-CCR2 axis. Blockade of the CCL2-CCR2 axis could reduce the incidence of tumors by preventing M2 macrophage recruitment and enhancing the efficacy of CD8⁺ T cells in the tumor microenvironment (Yang et al., 2020). Targeting immunosuppressive molecules on M2 macrophages is also an effective method. CSF-1R, a tyrosine kinase transmembrane receptor on M2 macrophages, is currently the most studied tumor-associated macrophages (TAMs) (Peyraud et al., 2017). Several clinical trials have reported a decrease in M2 macrophages and an increase in CD8/CD4+T cell ratio in advanced solid tumors treated with single or combined anti-CSF-1R agents (Strachan et al., 2013; Ries et al., 2014; Peyraud et al., 2017). Other promising targets include MerTK, Axl, and Tyro3, which have yielded encouraging results in preclinical studies (Myers et al., 2019).

Finally, we performed expression and functional analysis on the parent genes in AS events. Our results highlighted the role of SH3KBP1 in tumor immune regulation in HNSC. SH3KBP1, also called Cbl-interacting 85 kD (CIN85), encodes an adaptor protein that is involved in many signaling pathways, connecting multiple cellular compartments and processes (Havrylov et al., 2010). We observed that the expression of SH3KBP1 was higher in tumors than in normal tissues, indicating a poor prognosis in HNSC. Importantly, our correlation analysis found a correlation between its expression and the immune checkpoint expression. Although SH3KBP1 mediates the inhibition of T cell activation (Kong et al., 2019), the regulatory mechanism remains unclear. Yakymovych et al. proved that SH3KBP1 promoted the expression of TGFβ receptor on the cell surface and positively regulated TGFβ signaling (Yakymovych et al., 2015). TGFβ signaling induced the expression of both PD-1 and PD-L1 in infiltrating T cells in multiple preclinical models (Baas et al., 2016; Strait and Wang, 2020). Furthermore, a recent study demonstrated that TGFβ induced PD-L1 in vitro on human non-small cell lung cancer cell lines by Smad2-mediated canonical TGFβ signaling (David et al., 2017). Therefore, TGFβ signaling may regulate the expression of immune checkpoints via SH3KBP1, which should be verified in future studies.

Several AS-based prognostic signatures have been identified in HNSC (Xing et al., 2019; Li et al., 2020; Zhao et al., 2020). However, these studies have not focused on the prognostic signatures and the tumor microenvironment. In this research, the AS events were analyzed in a large-scale HNSC cohort using the TCGA database. We established a predictive AS event signature for the prognosis and treatment of HNSC. This research had several limitations. Due to the lack of RNA-Seq data with survival information in the public database, we were unable to perform survival analysis in the independent cohort. Thus, these results need further validation by other available datasets and future research. Secondly, the two risk groups exhibited distinct immune landscapes, but their differences were insignificant. Hence, greater subgroup differentiation is required for precision therapy. Despite these limitations, our signature can successfully differentiate between high- and low-risk patients and suggest accurate treatment strategies.

Conclusion

This research established a prognostic factor for HNSC patients according to the eight survival-associated AS events. The risk score was confirmed to be connected to HNSC prognosis and immune infiltration, suggesting that the prognostic signature could be an effective tool for the prognosis of HNSC patients and their clinical response to immunotherapy. However, this signature requires further research and validation in larger cohort studies.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

FY, PW, and YG designed the study. FY and YZ collected and assembled the data. FY, PW, YZ, GH, YT, ZL, and YG analyzed and interpreted data. YG and PW provided financial support. All authors read and approve the final manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This work was supported by the High-Level Hospital Program, Health Commission of Guangdong Province, China (No. HKUSZH201901033).

Acknowledgments

We acknowledge the cancer genome atlas and SpliceAid 2 database for providing their platforms and contributors for uploading their meaningful datasets.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.989081/full#supplementary-material

Abbreviations

HNSC, head and neck squamous cell carcinoma; AS, alternative splicing; TCGA, the cancer genome atlas; SF, splicing factors; LASSO, least absolute shrinkage and selection operator; HPV, human papilloma virus; OS, overall survival; ROC, receiver operating characteristic; GSEA, gene set enrichment analysis; TMB, tumor mutational burden.

References

Ariffin, A. B., Forde, P. F., Jahangeer, S., Soden, D. M., and Hinchion, J. (2014). Releasing pressure in tumors: What do we know so far and where do we go from here? A review. Cancer Res. 74 (10), 2655–2662. doi:10.1158/0008-5472.can-13-3696

PubMed Abstract | CrossRef Full Text | Google Scholar

Baas, M., Besançon, A., Goncalves, T., Valette, F., Yagita, H., Sawitzki, B., et al. (2016). Tgfβ-dependent expression of Pd-1 and Pd-L1 controls Cd8(+) T cell anergy in transplant tolerance. eLife 5, e08133. doi:10.7554/eLife.08133

PubMed Abstract | CrossRef Full Text | Google Scholar

Barbie, D. A., Tamayo, P., Boehm, J. S., Kim, S. Y., Moody, S. E., Dunn, I. F., et al. (2009). Systematic rna interference reveals that oncogenic kras-driven cancers require Tbk1. Nature 462 (7269), 108–112. doi:10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar

Berinstein, N. L., Wolf, G. T., Naylor, P. H., Baltzer, L., Egan, J. E., Brandwein, H. J., et al. (2012). Increased lymphocyte infiltration in patients with head and neck cancer treated with the irx-2 immunotherapy regimen. Cancer Immunol. Immunother. 61 (6), 771–782. doi:10.1007/s00262-011-1134-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Bindea, G., Mlecnik, B., Hackl, H., Charoentong, P., Tosolini, M., Kirilovsky, A., et al. (2009). Cluego: A Cytoscape plug-in to decipher functionally grouped gene Ontology and pathway annotation networks. Bioinformatics 25 (8), 1091–1093. doi:10.1093/bioinformatics/btp101

PubMed Abstract | CrossRef Full Text | Google Scholar

Bouzin, C., Brouet, A., De Vriese, J., Dewever, J., and Feron, O. (2007). Effects of vascular endothelial growth factor on the lymphocyte-endothelium interactions: Identification of caveolin-1 and nitric oxide as control points of endothelial cell anergy. J. Immunol. 178 (3), 1505–1511. doi:10.4049/jimmunol.178.3.1505

PubMed Abstract | CrossRef Full Text | Google Scholar

Cancer Genome Atlas Network (2015). Comprehensive genomic characterization of head and neck squamous cell carcinomas. Nature 517 (7536), 576–582. doi:10.1038/nature14129

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, D., Li, J., Liu, D., Hong, S., Qiao, Q., Sun, Q., et al. (2020). Tumor-expressed B7-H3 mediates the inhibition of antitumor T-cell functions in ovarian cancer insensitive to Pd-1 blockade therapy. Cell. Mol. Immunol. 17 (3), 227–236. doi:10.1038/s41423-019-0305-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., Lu, T., Zhong, F., Lv, Q., Fang, M., Tu, Z., et al. (2022). A signature of N(6)-methyladenosine regulator-related genes predicts prognoses and immune responses for head and neck squamous cell carcinoma. Front Immunol 13 1664–3224. doi:10.3389/fimmu.2022.809872

CrossRef Full Text | Google Scholar

Chen, Y. P., Wang, Y. Q., Lv, J. W., Li, Y. Q., Chua, M. L. K., Le, Q. T., et al. (2019). Identification and validation of novel microenvironment-based immune molecular subgroups of head and neck squamous cell carcinoma: Implications for immunotherapy. Ann. Oncol. 30 (1), 68–75. doi:10.1093/annonc/mdy470

PubMed Abstract | CrossRef Full Text | Google Scholar

Chow, L. Q. M., Haddad, R., Gupta, S., Mahipal, A., Mehra, R., Tahara, M., et al. (2016). Antitumor activity of pembrolizumab in biomarker-unselected patients with recurrent and/or metastatic head and neck squamous cell carcinoma: Results from the phase ib keynote-012 expansion cohort. J. Clin. Oncol. 34 (32), 3838–3845. doi:10.1200/jco.2016.68.1478

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, E. E., LaMonte, S. J., Erb, N. L., Beckman, K. L., Sadeghi, N., Hutcheson, K. A., et al. (2016). American cancer society head and neck cancer survivorship care guideline. Ca. Cancer J. Clin. 66 (3), 203–239. doi:10.3322/caac.21343

PubMed Abstract | CrossRef Full Text | Google Scholar

David, J. M., Dominguez, C., McCampbell, K. K., Gulley, J. L., Schlom, J., and Palena, C. (2017). A novel bifunctional anti-Pd-L1/Tgf-Β trap fusion protein (M7824) efficiently reverts mesenchymalization of human lung cancer cells. Oncoimmunology 6 (10), e1349589. doi:10.1080/2162402X.2017.1349589

PubMed Abstract | CrossRef Full Text | Google Scholar

De Stefani, A., Forni, G., Ragona, R., Cavallo, G., Bussi, M., Usai, A., et al. (2002). Improved survival with perilymphatic interleukin 2 in patients with resectable squamous cell carcinoma of the oral cavity and oropharynx. Cancer 95 (1), 90–97. doi:10.1002/cncr.10654

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferlay, J., Colombet, M., Soerjomataram, I., Mathers, C., Parkin, D. M., Pineros, M., et al. (2019). Estimating the global cancer incidence and mortality in 2018: GLOBOCAN sources and methods. Int. J. Cancer 144 (8), 1941–1953. doi:10.1002/ijc.31937

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferris, R. L., Blumenschein, G., Fayette, J., Guigay, J., Colevas, A. D., and Licitra, L., (2016). Nivolumab for recurrent squamous-cell carcinoma of the head and neck. N. Engl. J. Med. 375 (19), 1856–1867. doi:10.1056/NEJMoa1602252

PubMed Abstract | CrossRef Full Text | Google Scholar

Frankiw, L., Baltimore, D., and Li, G. (2019). Alternative mrna splicing in cancer immunotherapy. Nat. Rev. Immunol. 19 (11), 675–687. doi:10.1038/s41577-019-0195-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Gameiro, S. F., Ghasemi, F., Barrett, J. W., Koropatnick, J., Nichols, A. C., Mymryk, J. S., et al. (2018). Treatment-naïve Hpv+ head and neck cancers display a T-cell-inflamed phenotype distinct from their hpv- counterparts that has implications for immunotherapy. Oncoimmunology 7 (10), e1498439. doi:10.1080/2162402x.2018.1498439

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodman, A. M., Kato, S., Bazhenova, L., Patel, S. P., Frampton, G. M., Miller, V., et al. (2017). Tumor mutational burden as an independent predictor of response to immunotherapy in diverse cancers. Mol. Cancer Ther. 16 (11), 2598–2608. doi:10.1158/1535-7163.mct-17-0386

PubMed Abstract | CrossRef Full Text | Google Scholar

Griffioen, A. W., Damen, C. A., Blijham, G. H., and Groenewegen, G. (1996). Tumor angiogenesis is accompanied by a decreased inflammatory response of tumor-associated endothelium. Blood 88 (2), 667–673. doi:10.1182/blood.v88.2.667.bloodjournal882667

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, T., He, K., Wang, Y., Sun, J., Chen, Y., and Yang, Z. (2021). Prognostic signature of hepatocellular carcinoma and analysis of immune infiltration based on m6a-related lncrnas. Front. Oncol. 11, 691372. doi:10.3389/fonc.2021.691372

PubMed Abstract | CrossRef Full Text | Google Scholar

Havrylov, S., Redowicz, M. J., and Buchman, V. L. (2010). Emerging roles of ruk/cin85 in vesicle-mediated transport, adhesion, migration and malignancy. Traffic (Copenhagen, Den. 11 (6), 721–731. doi:10.1111/j.1600-0854.2010.01061.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jardim, D. L., Goodman, A., de Melo Gagliato, D., and Kurzrock, R. (2021). The challenges of tumor mutational burden as an immunotherapy biomarker. Cancer Cell 39 (2), 154–173. doi:10.1016/j.ccell.2020.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, D. E., Burtness, B., Leemans, C. R., Lui, V. W. Y., Bauman, J. E., and Grandis, J. R. (2020). Head and neck squamous cell carcinoma. Nat. Rev. Dis. Prim. 6 (1), 92. doi:10.1038/s41572-020-00224-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Joyce, J. A., and Fearon, D. T. (2015). T cell exclusion, immune privilege, and the tumor microenvironment. Science 348 (6230), 74–80. doi:10.1126/science.aaa6204

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, J. M., and Chen, D. S. (2016). Immune escape to Pd-L1/Pd-1 blockade: Seven steps to success (or failure). Ann. Oncol. 27 (8), 1492–1504. doi:10.1093/annonc/mdw217

PubMed Abstract | CrossRef Full Text | Google Scholar

Kong, M. S., Hashimoto-Tane, A., Kawashima, Y., Sakuma, M., Yokosuka, T., Kometani, K., et al. (2019). Inhibition of T Cell activation and function by the adaptor protein Cin85. Sci. Signal. 567, eaav4373. doi:10.1126/scisignal.aav4373

PubMed Abstract | CrossRef Full Text | Google Scholar

Larochelle, S. (2016). Systems biology. Protein isoforms: More than meets the eye. Nat. Methods 13 (4), 291. doi:10.1038/nmeth.3828

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Guo, G., Song, J., Cai, Z., Yang, J., Chen, Z., et al. (2017). B7-H3 promotes the migration and invasion of human bladder cancer cells via the pi3k/akt/stat3 signaling pathway. J. Cancer 8 (5), 816–824. doi:10.7150/jca.17759

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Chen, X., Wei, M., Liu, G., Tian, Y., Zhang, X., et al. (2020). Systemic analysis of rna alternative splicing signals related to the prognosis for head and neck squamous cell carcinoma. Front. Oncol. 10, 87. doi:10.3389/fonc.2020.00087

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z. X., Zheng, Z. Q., Wei, Z. H., Zhang, L. L., Li, F., Lin, L., et al. (2019). Comprehensive characterization of the alternative splicing landscape in head and neck squamous cell carcinoma reveals novel events associated with tumorigenesis and the immune microenvironment. Theranostics 9 (25), 7648–7665. doi:10.7150/thno.36585

PubMed Abstract | CrossRef Full Text | Google Scholar

Licitra, L., Mesia, R., Rivera, F., Remenár, É., Hitt, R., Erfán, J., et al. (2011). Evaluation of egfr gene copy number as a predictive biomarker for the efficacy of cetuximab in combination with chemotherapy in the first-line treatment of recurrent and/or metastatic squamous cell carcinoma of the head and neck: Extreme study. Ann. Oncol. 22 (5), 1078–1087. doi:10.1093/annonc/mdq588

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, B., Su, Q., Ma, J., Chen, C., Wang, L., Che, F., et al. (2021b). Prognostic value of eight-gene signature in head and neck squamous carcinoma. Front Oncol 11, 2234–943X. doi:10.3389/fonc.2021.657002

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H., Tekle, C., Chen, Y. W., Kristian, A., Zhao, Y., Zhou, M., et al. (2011). B7-H3 silencing increases paclitaxel sensitivity by abrogating jak2/stat3 phosphorylation. Mol. Cancer Ther. 10 (6), 960–971. doi:10.1158/1535-7163.mct-11-0072

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, S., Liang, J., Liu, Z., Zhang, C., Wang, Y., Watson, A. H., et al. (2021a). The role of Cd276 in cancers. Front. Oncol. 11, 654684. doi:10.3389/fonc.2021.654684

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, H., Ma, C., Shao, J., and Cao, J. (2020). Prognostic implications of novel ten-gene signature in uveal melanoma. Front. Oncol. 10, 567512. doi:10.3389/fonc.2020.567512

PubMed Abstract | CrossRef Full Text | Google Scholar

Majzner, R. G., Theruvath, J. L., Nellan, A., Heitzeneder, S., Cui, Y., Mount, C. W., et al. (2019). Car T cells targeting B7-H3, a pan-cancer antigen, demonstrate potent preclinical activity against pediatric solid tumors and brain tumors. Clin. Cancer Res. 25 (8), 2560–2574. doi:10.1158/1078-0432.ccr-18-0432

PubMed Abstract | CrossRef Full Text | Google Scholar

Mandal, R., Şenbabaoğlu, Y., Desrichard, A., Havel, J. J., Dalin, M. G., Riaz, N., et al. (2016). The head and neck cancer immune landscape and its immunotherapeutic implications. JCI Insight 1 (17), e89829. doi:10.1172/jci.insight.89829

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, L., Fan, T. F., Wu, L., Yu, G. T., Deng, W. W., Chen, L., et al. (2017). Selective blockade of B7-H3 enhances antitumour immune activity by reducing immature myeloid cells in head and neck squamous cell carcinoma. J. Cell. Mol. Med. 21 (9), 2199–2210. doi:10.1111/jcmm.13143

PubMed Abstract | CrossRef Full Text | Google Scholar

Myers, K. V., Amend, S. R., and Pienta, K. J. (2019). Targeting Tyro3, Axl and mertk (tam receptors): Implications for macrophages in the tumor microenvironment. Mol. Cancer 18 (1), 94. doi:10.1186/s12943-019-1022-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12 (5), 453–457. doi:10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

Ott, P. A., Hu, Z., Keskin, D. B., Shukla, S. A., Sun, J., Bozym, D. J., et al. (2017). An immunogenic personal neoantigen vaccine for patients with melanoma. Nature 547 (7662), 217–221. doi:10.1038/nature22991

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, Y., Yu, Y., Wang, X., and Zhang, T. (2020). Tumor-associated macrophages in tumor immunity. Front. Immunol. 11, 583084. doi:10.3389/fimmu.2020.583084

PubMed Abstract | CrossRef Full Text | Google Scholar

Peranzoni, E., Lemoine, J., Vimeux, L., Feuillet, V., Barrin, S., Kantari-Mimoun, C., et al. (2018). Macrophages impede Cd8 T cells from reaching tumor cells and limit the efficacy of anti-Pd-1 treatment. Proc. Natl. Acad. Sci. U. S. A. 115 (17), E4041–e50. doi:10.1073/pnas.1720948115

PubMed Abstract | CrossRef Full Text | Google Scholar

Peyraud, F., Cousin, S., and Italiano, A. (2017). Csf-1r inhibitor development: Current clinical status. Curr. Oncol. Rep. 19 (11), 70. doi:10.1007/s11912-017-0634-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ribas, A., Dummer, R., Puzanov, I., VanderWalde, A., Andtbacka, R. H. I., Michielin, O., et al. (2018). Oncolytic virotherapy promotes intratumoral T cell infiltration and improves anti-Pd-1 immunotherapy. Cell 174 (4), 1031–1032. doi:10.1016/j.cell.2018.07.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Ries, C. H., Cannarile, M. A., Hoves, S., Benz, J., Wartha, K., Runza, V., et al. (2014). Targeting tumor-associated macrophages with anti-csf-1r antibody reveals a strategy for cancer therapy. Cancer Cell 25 (6), 846–859. doi:10.1016/j.ccr.2014.05.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryan, M. C., Cleland, J., Kim, R., Wong, W. C., and Weinstein, J. N. (2012). Spliceseq: A resource for analysis and visualization of rna-seq data on alternative splicing and its functional impacts. Bioinformatics 28 (18), 2385–2387. doi:10.1093/bioinformatics/bts452

PubMed Abstract | CrossRef Full Text | Google Scholar

Sciarrillo, R., Wojtuszkiewicz, A., Assaraf, Y. G., Jansen, G., Kaspers, G. J. L., Giovannetti, E., et al. (2020). The role of alternative splicing in cancer: From oncogenesis to drug resistance. Drug resist. updat. 53, 100728. doi:10.1016/j.drup.2020.100728

PubMed Abstract | CrossRef Full Text | Google Scholar

Seiwert, T. Y., Burtness, B., Mehra, R., Weiss, J., Berger, R., Eder, J. P., et al. (2016). Safety and clinical activity of pembrolizumab for treatment of recurrent or metastatic squamous cell carcinoma of the head and neck (Keynote-012): An open-label, multicentre, phase 1b trial. Lancet. Oncol. 17 (7), 956–965. doi:10.1016/s1470-2045(16)30066-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Selby, M. J., Engelhardt, J. J., Quigley, M., Henning, K. A., Chen, T., Srinivasan, M., et al. (2013). Anti-Ctla-4 antibodies of Igg2a isotype enhance antitumor activity through reduction of intratumoral regulatory T cells. Cancer Immunol. Res. 1 (1), 32–42. doi:10.1158/2326-6066.cir-13-0013

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, T., Ma, Y., Cao, L., Zhan, S., Xu, Y., Fu, F., et al. (2019). B7-H3 promotes aerobic glycolysis and chemoresistance in colorectal cancer cells by regulating Hk2. Cell Death Dis. 10 (4), 308. doi:10.1038/s41419-019-1549-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Shibata, H., Zhou, L., Xu, N., Egloff, A. M., and Uppaluri, R. (2021). Personalized cancer vaccination in head and neck cancer. Cancer Sci. 112 (3), 978–988. doi:10.1111/cas.14784

PubMed Abstract | CrossRef Full Text | Google Scholar

Steinman, R. M. (2012). Decisions about dendritic cells: Past, present, and future. Annu. Rev. Immunol. 30, 1–22. doi:10.1146/annurev-immunol-100311-102839

PubMed Abstract | CrossRef Full Text | Google Scholar

Strachan, D. C., Ruffell, B., Oei, Y., Bissell, M. J., Coussens, L. M., Pryer, N., et al. (2013). CSF1R inhibition delays cervical and mammary tumor growth in murine models by attenuating the turnover of tumor-associated macrophages and enhancing infiltration by CD8+ T cells. Oncoimmunology 2 (12), e26968. doi:10.4161/onci.26968

PubMed Abstract | CrossRef Full Text | Google Scholar

Strait, A. A., and Wang, X-J. (2020). The role of transforming growth factor-beta in immune suppression and chronic inflammation of squamous cell carcinomas. Mol. Carcinog. 59 (7), 745–753. doi:10.1002/mc.23196

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugiyama, D., Nishikawa, H., Maeda, Y., Nishioka, M., Tanemura, A., Katayama, I., et al. (2013). Anti-Ccr4 mab selectively depletes effector-type Foxp3+Cd4+ regulatory T cells, evoking antitumor immune responses in humans. Proc. Natl. Acad. Sci. U. S. A. 110 (44), 17945–17950. doi:10.1073/pnas.1316796110

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, D. S. W., Chong, F. T., Leong, H. S., Toh, S. Y., Lau, D. P., Kwang, X. L., et al. (2017). Long noncoding rna egfr-as1 mediates epidermal growth factor receptor addiction and modulates treatment response in squamous cell carcinoma. Nat. Med. 23 (10), 1167–1175. doi:10.1038/nm.4401

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, Y. S., Sansanaphongpricha, K., Prince, M. E. P., Sun, D., Wolf, G. T., and Lei, Y. L. (2018). Engineering vaccines to reprogram immunity against head and neck cancer. J. Dent. Res. 97 (6), 627–634. doi:10.1177/0022034518764416

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, X., Liu, F., Liu, Z., Cao, Y., Zhang, Z., Wang, Y., et al. (2020). Bioactivity and safety of B7-H3-targeted chimeric antigen receptor T cells against anaplastic meningioma. Clin. Transl. Immunol. 9 (6), e1137. doi:10.1002/cti2.1137

PubMed Abstract | CrossRef Full Text | Google Scholar

Tashiro, H., and Brenner, M. K. (2017). Immunotherapy against cancer-related viruses. Cell Res. 27 (1), 59–73. doi:10.1038/cr.2016.153

PubMed Abstract | CrossRef Full Text | Google Scholar

Turley, S. J., Cremasco, V., and Astarita, J. L. (2015). Immunological hallmarks of stromal cells in the tumour microenvironment. Nat. Rev. Immunol. 15 (11), 669–682. doi:10.1038/nri3902

PubMed Abstract | CrossRef Full Text | Google Scholar

Ule, J., and Blencowe, B. J. (2019). Alternative splicing regulatory networks: Functions, mechanisms, and evolution. Mol. Cell 76 (2), 329–345. doi:10.1016/j.molcel.2019.09.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, C., Li, Y., Jia, L., Kim, J. K., Li, J., Deng, P., et al. (2021). Cd276 expression enables squamous cell carcinoma stem cells to evade immune surveillance. Cell stem Cell 28 (9), 1597–1613. doi:10.1016/j.stem.2021.04.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, M., Zhou, Z., Zheng, J., Xiao, W., Zhu, J., Zhang, C., et al. (2021). Identification and validation of a prognostic immune-related alternative splicing events signature for glioma. Front. Oncol. 11, 650153. doi:10.3389/fonc.2021.650153

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z. A-O., Yuan, H., Huang, J., Hu, D., Qin, X., Sun, C., et al. (2021c). Prognostic value of immune-related genes and immune cell infiltration analysis in the tumor microenvironment of head and neck squamous cell carcinoma. Head Neck 43 182–197. doi:10.1002/hed.26474

PubMed Abstract | CrossRef Full Text | Google Scholar

Woods, D. M., Ramakrishnan, R., Laino, A. S., Berglund, A., Walton, K., Betts, B. C., et al. (2018). Decreased suppression and increased phosphorylated Stat3 in regulatory T cells are associated with benefit from adjuvant Pd-1 blockade in resected metastatic melanoma. Clin. Cancer Res. 24 (24), 6236–6247. doi:10.1158/1078-0432.ccr-18-1100

PubMed Abstract | CrossRef Full Text | Google Scholar

Xing, L., Zhang, X., and Tong, D. (2019). Systematic profile Analysis of prognostic alternative messenger rna splicing signatures and splicing factors in head and neck squamous cell carcinoma. DNA Cell Biol. 38 (7), 627–638. doi:10.1089/dna.2019.4644

PubMed Abstract | CrossRef Full Text | Google Scholar

Xue, Y., Tong, L., LiuAnwei Liu, F., Liu, A., Zeng, S., Xiong, Q., et al. (2019). Tumor-infiltrating M2 macrophages driven by specific genomic alterations are associated with prognosis in bladder cancer. Oncol. Rep. 42 (2), 581–594. doi:10.3892/or.2019.7196

PubMed Abstract | CrossRef Full Text | Google Scholar

Yakymovych, I., Yakymovych, M., Zang, G., Mu, Y., Bergh, A., Landström, M., et al. (2015). Cin85 modulates tgfβ signaling by promoting the presentation of tgfβ receptors on the cell surface. J. Cell Biol. 210 (2), 319–332. doi:10.1083/jcb.201411025

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, H., Zhang, Q., Xu, M., Wang, L., Chen, X., Feng, Y., et al. (2020). Ccl2-Ccr2 Axis recruits tumor associated macrophages to induce immune evasion through Pd-1 signaling in esophageal carcinogenesis. Mol. Cancer 19 (1), 41. doi:10.1186/s12943-020-01165-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yonesaka, K., Haratani, K., Takamura, S., Sakai, H., Kato, R., Takegawa, N., et al. (2018). B7-H3 negatively modulates ctl-mediated cancer immunity. Clin. Cancer Res. 24 (11), 2653–2664. doi:10.1158/1078-0432.ccr-17-2852

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, S., Hu, C., Liu, L., Cai, L., Du, X., Yu, Q., et al. (2020). Comprehensive analysis and establishment of a prediction model of alternative splicing events reveal the prognostic predictor and immune microenvironment signatures in triple negative breast cancer. J. Transl. Med. 18 (1), 286. doi:10.1186/s12967-020-02454-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, X., Si, S., Li, X., Sun, W., and Cui, L. (2020). Identification and validation of an alternative splicing-based prognostic signature for head and neck squamous cell carcinoma. J. Cancer 11 (15), 4571–4580. doi:10.7150/jca.44746

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, W. T., and Jin, W. L. (2021). B7-H3/Cd276: An emerging cancer immunotherapy. Front. Immunol. 12, 701006. doi:10.3389/fimmu.2021.701006

PubMed Abstract | CrossRef Full Text | Google Scholar