Investigation of an FGFR-Signaling-Related Prognostic Model and Immune Landscape in Head and Neck Squamous Cell Carcinoma

Background: There is accumulating evidence on the clinical importance of the fibroblast growth factor receptor (FGFR) signal, hypoxia, and glycolysis in the immune microenvironment of head and neck squamous cell carcinoma (HNSCC), yet reliable prognostic signatures based on the combination of the fibrosis signal, hypoxia, and glycolysis have not been systematically investigated. Herein, we are committed to establish a fibrosis–hypoxia–glycolysis–related prediction model for the prognosis and related immune infiltration of HNSCC.

Methods: Fibrotic signal status was estimated with microarray data of a discovery cohort from the TCGA database using the UMAP algorithm. Hypoxia, glycolysis, and immune-cell infiltration scores were imputed using the ssGSEA algorithm. Cox regression with the LASSO method was applied to define prognostic genes and develop a fibrosis–hypoxia–glycolysis–related gene signature. Immunohistochemistry (IHC) was conducted to identify the expression of specific genes in the prognostic model. Protein expression of several signature genes was evaluated in HPA. An independent cohort from the GEO database was used for external validation. Another scRNA-seq data set was used to clarify the related immune infiltration of HNSCC.

Results: Six genes, including AREG, THBS1, SEMA3C, ANO1, IGHG2, and EPHX3, were identified to construct a prognostic model for risk stratification, which was mostly validated in the independent cohort. Multivariate analysis revealed that risk score calculated by our prognostic model was identified as an independent adverse prognostic factor (p < .001). Activated B cells, immature B cells, activated CD4⁺ T cells, activated CD8⁺ T cells, effector memory CD8⁺ T cells, MDSCs, and mast cells were identified as key immune cells between high- and low-risk groups. IHC results showed that the expression of SEMA3C, IGHG2 were slightly higher in HNSCC tissue than normal head and neck squamous cell tissue. THBS1, ANO1, and EPHX3 were verified by IHC in HPA. By using single-cell analysis, FGFR-related genes and highly expressed DEGs in low-survival patients were more active in monocytes than in other immune cells.

Conclusion: A fibrosis–hypoxia–glycolysis–related prediction model provides risk estimation for better prognoses to patients diagnosed with HNSCC.

Introduction

Head and neck squamous cell carcinoma (HNSCC) has a worldwide incidence of more than 600,000 cases per year (Ferlay et al., 2015), including a heterogeneous group of tumors that arise from the oral cavity, oropharynx, larynx, hypopharynx, nasopharynx, and sinonasal cavity (Kim et al., 2021). Due to its special anatomical location, patients with HNSCC often experience vital dysfunction, especially in the aspects of swallowing, feeding, breathing, and psychological health (Zhu et al., 2017). Despite significant progress in available therapies, the 5-year overall survival (OS) rate of HNSCC patients has not obviously improved in recent decades (Siegel et al., 2019). Therefore, there is an urgent need for a way to predict the progression of HNSCC (Shield et al., 2017).

The fibroblast growth factor receptor (FGFR) is a receptor tyrosine kinase (RTK) signaling pathway involved in the regulation of angiogenesis, invasion, and metastasis of tumors (Chae et al., 2017). It is considered to be a contributory factor in the development of fibrosis, causing the exacerbation of liver fibrosis (Seitz and Hellerbrand, 2021), systemic sclerosis (SS) (Chakraborty et al., 2020), idiopathic pulmonary fibrosis (IPF) (Wollin et al., 2015), and so on. A recent study also revealed the clinically activity of an FGFR inhibitor against HNSCC (Schuler et al., 2019). However, few investigations focus on the potential mechanisms about this profibrotic mediator in HNSCC.

Hypoxia, as a hallmark of tumor, is a potent microenvironmental factor facilitating proliferation and progression in a variety of cancers (Rankin and Giaccia, 2016). Previous research suggests that hypoxic HNSCC cells trigger glycolysis to obtain energy and balance metabolic and bioenergetic (Zhu et al., 2017). Furthermore, activated fibroblasts synthesized excessive collagen, leading to a microenvironment of hypoxia, which might deteriorate the progression of disease, whereas there were few records on the prognosis of HNSCC by combining a profibrotic signal with hypoxia and glycolysis. The immune landscape associated with the factors mentioned above also demands exploration.

Therefore, in this study, we establish an FGFR-signaling–hypoxia–glycolysis–related prediction model for the prognosis of HNSCC with a series of bioinformatics analyses. Immune cell infiltration associated with prognosis and profibrotic signaling is also revealed.

Methods

Patient Cohort and Data Preparation

The discovery cohort of the study contained 483 HNSCC patients from the Cancer Genome Atlas (TCGA, available at https://portal.gdc.cancer.gov/) data set. To obtain a validation cohort, RNA-seq data and related clinicopathological data were downloaded from the Gene Expression Omnibus (GEO, available at https://www.ncbi.nlm.nih.gov/geo/) database (GSE41613), including 97 patients with HPV-negative oral squamous cell carcinoma (OSCC). The microarray data of GSE41613 was built upon the GPL570 Platform (Affymetrix Human Genome U133 Plus 2.0 Array). A single-cell RNA sequencing (scRNA-seq) data set (GSE139324) was also used for analyzing tumor immune infiltration in HNSCC patients, including tumor-infiltrating immune cells from 18 HNSCC patients and tissue resident immune cells from five healthy donor tonsils (Cillo et al., 2020).

The mRNA expression profiles from the TCGA database were normalized using fragments per kilobase of exon per million reads mapped (FPKM). Background correction and normalization has been performed for each series before downloading GSE41613 from the GEO database. Additionally, the harmony algorithm was used in the integration of the single-cell RNA-seq data set GSE139324 considering biological and technical differences (Korsunsky et al., 2019). The general idea and methodologies of our study are shown in a flowchart (Figure 1).

FIGURE 1

FIGURE 1. Flowchart diagram of the study.

Distinction of Fibrotic Signal Status and FGFR-Related DEGs

An algorithm of uniform manifold approximation and projection (UMAP) was applied to deduce the fibrotic signal status of HNSCC patients. Based on the given hallmarks or signatures, UMAP, a nonlinear reductive dimension method, is able to assign a group of patients to diverse clusters. The gene set of the FGFR signal pathway was downloaded from the Molecular Signatures Database (MSigDB version 6.0), identifying the relative activation degree of fibrotic signal in patients. Based on the limma algorithm (Korsunsky et al., 2019) and functional enrichment analysis, two clusters including “fibrosis^low” and “fibrosis^high” were identified to estimate the fibrotic signal status. The limma algorithm was applied to identify DEGs between the two clusters based on the standards of false discovery rate (FDR) adjusted p-value <.0001 and | log2(Fold change) | > 1. To confirm biological functions and pathway enrichment of the fibrosis-related DEGs, we performed Gene Ontology (GO) functional analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using the “Metascape” website (https://metascape.org/) (Zhou et al., 2019).

Distinction of Hypoxia–Glycolysis–Related DEGs

The ssGSEA algorithm was applied to explore the hypoxia and glycolysis degree in the HNSCC expression profile of the TCGA database. According to hypoxia and glycolysis scores calculated previously, two groups of patients were stratified. An optimal cutting point for classifying was determined by maximally selected rank statistics using the “survival” and “survminer” R package. Subsequently, “hypoxia^high,” “hypoxia^low,” “glycolysis^high,” and “glycolysis^low” groups were identified, respectively. We further considered hypoxia and glycolysis together by combining them into a two-dimensional index; that is, patients were divided into three groups, i.e., hypoxia^low/glycolysis^low, hypoxia^high/glycolysis^low, and “mix” groups. Hypoxia–glycolysis–related DEGs were identified based on the standards of FDR adjusted p-value <.0001 and | log2(Fold change) | > 1 using the R package “limma”.

Prognosis Prediction Model of HNSCC Based on Fibrosis–Hypoxia–Glycolysis–Related DEGs

Hypoxia–glycolysis–related DEGs and fibrosis–related DEGs were intersected to obtain the 39 shared fibrosis–hypoxia–glycolysis–related DEGs by Venn analysis. To obtain prognostic shared DEGs, univariate Cox regression analyses were further performed among all 39 DEGs mentioned above to screen those risk or protective shared DEGs with p < .05. Thereafter, we applied the least absolute shrinkage and selection operator (LASSO) (Friedman et al., 2010; Liu et al., 2013) to preserve valuable variables in 21 prognostic shared DEGs, implementing a high-dimensional prediction and avoiding overfitting. The LASSO Cox regression model was scientifically built up depending on threefold cross-validation and 1000 iterations, which decreased the underlying instability of the results. The optimal tuning parameter λ was identified via 1-SE (standard error) criterion. Eventually, the selected prognostic gene signatures were used to establish the prognosis prediction model of HNSCC based on fibrosis–hypoxia–glycolysis–related DEGs. The risk score computing formula is:

$Risk score={\displaystyle \sum }_{i=1}^n\left(coefficien{t}_i\times \exp ression of signature gen{e}_i\right).$

Based on the risk scores, we computed the optimal cutting point to stratify HNSCC patients into high- and low-risk groups.

Identification of Immune Cell Infiltration Status

To predict the immune cell infiltration status, the ssGSEA algorithm was applied to identify the abundance of 28 immunocytes in each HNSCC patient, confirming the underlying correlation between immune infiltration status and prognostic model. The gene set from Charoentong et al. was obtained to calculate ssGSEA scores for immune cell populations. (Charoentong et al., 2017).

Evaluation of Immunohistochemical Staining

Formalin-fixed, paraffin-embedded tumor tissues were collected from eight patients with HNSCC diagnosed at the 3rd Xiangya Hospital from January 2020 to December 2020. This study was approved by the Ethics Committee of the 3rd Xiangya Hospital (No: 21158). The immunohistochemical process was performed as described previously (Guan et al., 2020). Our procedure used the following antibodies: Polyclonal rabbit anti-Semaphorin 3c (1:200 dilution; ab135842; Abcam Biochemicals, UK); monoclonal rabbit anti-human IgG2 (1:1000 dilution; ab134050; Abcam Biochemicals, UK). IHC results for IGHG2 and SEMA3C were evaluated by computerizing optical density (OD) measurements using ImageJ software, which depends on the degree and area of staining. Samples were scored by two trained pathologists according to the percentage contribution of high positive, positive, low positive, and negative. The immunoreactive score (IRS) was evaluated as follows: 4, high positive; 3, positive; 2, low positive; and 1 negative (Varghese et al., 2014).

The Human Protein Atlas (HPA, https://www.ptroteinatlas.org/) provides us the IHC staining data in HNSCC and normal tissue (Uhlen et al., 2010). The expression level of target protein was classified into high, medium, low, and not detected according to degree of staining (strong, moderate, weak, or negative) and the proportion of stained cells (>75%, 25%–75%, or <25%).

Single-Cell Analysis of Tumor Infiltrating Immune Cells From HNSCC Patients

After dimension reduction through principal component analysis (PCA), the t-distributed stochastic neighbor embedding (t-SNE) algorithm (Kobak and Berens, 2019), a technique that maps a set of high-dimensional points to two dimensions, was used to compute the degree of similarity between cells, which is visualized by the distances among the plotted points on the graph. It also potently governed how many of its nearest neighbors each point is attracted to. Here, each point represented a cell. The scores of individual cells for pathway activities were estimated by the R package “AUCell.” According to gene expression rankings in each cell for a certain gene set, area under the curve (AUC) values were calculated to represent the proportion of top-ranking genes in the gene set for each cell (Corridoni et al., 2020).

Statistical Analysis

Using R version 4.0.2 (www.r-project.org/) and the appropriate packages, all statistical analysis was carried out. We implemented the UMAP algorithm using R package “umapr” for nonlinear dimension reduction and the ssGSEA algorithm using R package “GSVA” for the hypoxia and glycolysis score. The Lasso Cox regression model was conducted, and standard statistical tests were guaranteed by using the R package “glmnet.” The FDR method was performed to adjust multiple tests. Risk factors were eventually identified through multivariate Cox regression analysis.

Results

Fibrosis Signal and Fibrosis-Related DEGs in HNSCC

The expression profiles and clinical information of 483 HNSCC patients were contained in the discovery cohort downloaded from the TCGA database. The clinical information of patients is shown in Table 1. Seventy genes positively correlated with the FGFR signaling pathway were used to evaluate the status of fibrosis signal activation in patients. Based on the algorithm UMAP, we divided the patients into two clusters using the fibrosis-related expression matrix, enabling us to assign each patient to the nearest cluster (Figure 2A). A Kaplan–Meier plot demonstrates that significant differences in survival were witnessed between the two clusters (Figure 2B). There were 309 and 176 patients included in clusters 1 and 2, respectively. To obtain fibrosis-related DEGs, we compared the expression profiles between the clusters. A total of 187 fibrosis-related DEGs overexpressed in cluster 1, where patients had worse survival, which were enriched in “response to wounding” (Figure 2C), “TGF-beta signaling pathway” (Figure 2D). This implied the patients in cluster 1 were in a higher state of fibrosis activation. Enrichment analysis showed 186 DEGs overexpressed in cluster 2 were enriched in “immunoglobulin complex” (Figure 2E), “metabolism of xenobiotics cytochrome P450” (Figure 2F). These findings are consistent with the previous research that patients with good immune status have a better prognosis.

TABLE 1

TABLE 1. Basic information of HNSCC patients in discovery cohort.

FIGURE 2

FIGURE 2. Grouping of patients according to FGFR signaling pathway. (A) The UMAP algorithm classifies HNSCC patients into two clusters, indicated by different colors. (B) The survivorship analysis plot of OS in two clusters. GO function enrichment (C,E) and KEGG pathway enrichment analysis (D,F) of differentially expressed genes in the fibrosis^low-survival group, colored by p-values. (C,D) Analysis of the high expression genes in the fibrosis^low-survival group. (E,F) Analysis of the low expression genes in the fibrosis^low-survival group.

Hypoxia Status, Glycolysis Status, and Hypoxia–Glycolysis–Related DEGs in HNSCC

Meanwhile, using the “GSVA” package, the ssGSEA algorithm was implemented to quantify the hypoxia or glycolysis enrichment score of each HNSCC patient in hypoxia or glycolysis hallmark genes from the MSigDB. To identify the effect of hypoxia and glycolysis on prognosis, univariate Cox regression analyses were further performed among patients’ hypoxia and glycolysis scores. Hypoxia and glycolysis, as illustrated in the forest diagram in Figure 3A, were considered risk factors for prognosis in HNSCC patients. Based on maximally selected rank statistics, we divided patients into two groups according to hypoxia (Figure 3B) and glycolysis (Figure 3C) scores. We further synthesized the hypoxia and glycolysis status into a two-dimensional index, dividing patients into three groups, i.e., hypoxia^high/glycolysis^high, hypoxia^low/glycolysis^low, and “mix” groups. Significant differences in survival were observed among three groups (Figure 3D, log rank test, p < .0001). The survivorship analysis (Kaplan–Meier) showed a better survival in the hypoxia^high/glycolysis^high group than the hypoxia^low/glycolysis^low group as we expected (Figure 3E). Besides this, the mix group was at an intermediate level. A total of 108 hypoxia–glycolysis–related DEGs were obtained after comparing expression profiles between hypoxia^high/glycolysis^high and hypoxia^low/glycolysis^low.

FIGURE 3

FIGURE 3. Grouping of patients according to their hypoxia and glycolysis status. (A) Forest plot of hypoxia and glycolysis scores by univariate Cox regression. (B,C) A vertical dashed line was used to mark the cutoff point with the maximum standard log-rank statistic based on hypoxia and glycolysis scores. (D) Kaplan–Meier plot of OS among hypoxia^high/glycolysis^high, hypoxia^low/glycolysis^low, and mix groups. (E) Kaplan–Meier plot of OS between hypoxia^high/glycolysis^high and hypoxia^low/glycolysis^low.

Construction of the Fibrosis–Hypoxia–Glycolysis–Related Prognostic Model in the TCGA Data set

The above 334 fibrosis-related and 108 hypoxia–glycolysis–related DEGs screened from HNSCC were intersected to obtain the 39 shared genes (Figure 4A). To further filtrate the prognostic DEGs, univariate Cox regression analysis was conducted on 39 shared DEGs, and 21 DEGs with p < .05 were identified (Figure 4B). Among them, most of them (18 out of 21, 85.7%) were risk DEGs. Six critical variables were selected from the above 21 prognostic DEGs using the LASSO regression method, among which four were risk DEGs and two were protective (Figures 4C,D). For each HNSCC patient, a risk score was calculated based on the expression levels of the six characteristic DEGs and corresponding coefficients from the LASSO Cox regression: risk score = 0.0369 × expression of AREG+ 0.03432 × expression of THBS1+ 0.02182 × expression of SEMA3C+ 0.07125 × expression of ANO1+(-0.07718) × expression of IGHG2+ (-0.09177) × expression of EPHX3. Using the maximum selective rank method as the basis of demarcation, patients were divided into high- and low-risk groups according to their risk scores (Figure 4E). The low-risk group showed a significantly better effect on prognosis compared with the high-risk group (Figure 4F, log rank test, p < .0001).

FIGURE 4

FIGURE 4. Constructing a prognostic survival model in patients with HNSCC. (A) Venn diagrams show the fibrosis–hypoxia–glycolysis related DEGs (39). (B) Forest plot of the prognostic effects of 21 DEGs with p < .05 on a univariate Cox regression analysis. (C) LASSO coefficient profiles of 21 screened DEGs. (D) Threefold cross-validation for LASSO analysis was performed to calculate the minimum criteria. (E) A vertical dashed line was used to mark the cutoff point with the maximum standard log-rank statistic. (F) Kaplan–Meier plot established the survival differences between high- and low-risk groups.

Supplementary Information on Prognostic Model and its Relationship With Immunocyte Infiltration

We performed survival analysis on 483 HNSCC patients, which reorganized according to the location of the primary tumor in an attempt to verify the reliability of the prognostic model. In several regions with high incidence of HNSCC, such as tongue and larynx (Figures 5B,D), survival comparison revealed that the high-risk group was associated with a worse prognosis of the patients. The same went for tonsil, hypopharynx, and so on (Figures 5A,C). Univariate Cox regression analyses indicate that the risk score, similar to other clinical characteristics, such as age, could be deemed as an independent risk factor to assess the prognosis of patients with HNSCC (Figure 5E). Additionally, ssGSEA was used to estimate the immune cell infiltration in the patients. To explore the correlation between prognostic model and immune cell infiltration, correlation analysis was performed between six optimal prognostic signatures and the immunocyte infiltration score (Figure 5F). Furthermore, significantly decreased infiltration of eight specific immune cells was observed in the high- compared with the low-risk group, that is, activated B cells, immature B cells, activated CD4⁺ T cells, activated CD8⁺ T cells, effector memory CD8⁺ T cells, MDSCs, and mast cells (Figure 5G).

FIGURE 5

FIGURE 5. (A–D) A Kaplan–Meier plot establishes the survival differences between high- and low-risk groups in HNSCC with different primary sites. (E) Forest plot of other clinical characteristic and risk scores on a univariate Cox regression analysis. (F) Heatmap plotted by the correlations between the expression of genes in prognostic model and immune infiltrate level in HNSCC. (G) Comparison of immunocyte infiltration between high- and low-risk groups.

External Independent Cohort Validation of the Fibrosis–Hypoxia–Glycolysis–Related Prognostic Model in the GEO Data set

The fibrosis–hypoxia–glycolysis–based prognosis model was further validated in an independent cohort “GSE41613.” Searching for six genes from the prognosis model in the expression matrix of 97 OSCC patients, five of them were found. The clinical information of patients is shown in Table 2. Univariate Cox regression analyses confirmed that AREG, THBS1, SEMA3C, and ANO1 were risk factors in the prognosis of patients. On the contrary, EPHX3 was a protective one, which is consistent with the model previously conducted (Figure 6A). Based on maximally selected rank statistics, patients were classified into two groups according to the level of each gene expression (Figures 6B,D,F,H,J), and survival analysis was performed (Figures 6C,E,G,I,K). Kaplan–Meier curves showing that patients with higher expression levels of four risk genes in the prognostic model had worse OS (Figures 6C,E,G,I), whereas the opposite was true for protective genes (Figure 6K).

TABLE 2

TABLE 2. Basic information of OSCC patients in validation cohort.

FIGURE 6

FIGURE 6. Validation in external data set (GSE416130). (A) Forest plot of the prognostic effects of five genes in prognostic model on a univariate Cox regression analysis. (B,D,F,H,J) A vertical dashed line was used to mark the cutoff point with the maximum standard log-rank statistic based on the expression of each gene. (C,E,G,I,K) Kaplan–Meier plot established the survival differences between high (high expression of target gene) and low (low expression of target gene) groups.

IHC Analysis of Prognostic Signatures

To further validate the expression of prognosis-related molecules in HNSCC specimens and normal squamous epithelium of head and neck, we performed IHC staining analysis of paraffin section of HNSCC. IHC staining analysis suggested that the expression of SEMA3C and IGHG2 were slightly higher in HNSCC tissue quantified by the antibodies ab135842 (Figures 7A–F) and ab134050 (Figures 8A–F).

FIGURE 7

FIGURE 7. IHC analysis of the protein expression of SEMA3C in HNSCC and normal tissues. (A–E) The expression of SEMA3C was detected by IHC in five patients with HNSCC (Magnification ×200). (F) The expression of SEMA3C was detected by IHC in normal head and neck squamous cell tissue (Magnification ×200).

FIGURE 8

FIGURE 8. IHC analysis of the protein expression of IGHG2 in HNSCC and normal tissues. (A–E) The expression of IGHG2 was detected by IHC in five patients with HNSCC (Magnification ×200). (F) The expression of IGHG2 was detected by IHC in normal head and neck squamous cell tissue (Magnification ×200).

According to the protein expression data from the HPA, we compared the protein expression of six-gene signatures in HNSCC tissue and squamous epithelium normally located in the head and neck, such as oral squamous epithelium. We preliminarily inferred that the protein expression of these genes differed between HNSCC and normal tissues. The detailed results are presented in Supplementary Figure S1.

Fibrosis-Related Immune Landscape of HNSCC Patients Based on scRNA-Seq

To further understand the correlation between the fibrosis signal and immunity, we analyzed scRNA-seq data downloaded from the GEO data set “GSE139324.” A total of 23 samples were used for analysis. Of these, 18 samples were tumor-infiltrating immune cells from HNSCC patients (HPV negative), and five were tissue-resident immune cells from healthy donor tonsils. After integrating data by the harmony algorithm and binning by the t-SNE algorithm, 20 cell clusters were successfully classified. Moreover, there was a significant difference in the number of each cell subset between HNSCC patients and healthy donors (HD) (Figure 9A). By the expression of several cell surface and intracellular markers, clusters 4 and 11 were defined as monocytes. Clusters 3, 9, 13, 14, and 20 expressed markers associated with B cells (e.g., CD79A), and the rest of the clusters expressed genes associated with T cells (e.g., CD3D) (Figure 9B). The cells with a high expression of FGF-receptor-signaling–related genes and highly expressed DEGs in low-survival patients were severally highlighted in Figures 9C,D, most of which were highly expressed in clusters 4 and 11. An obviously higher composition of monocytes is shown in patients compared with healthy controls (Figure 9E). Significant differences in the number of cells in clusters 4 and 11 were witnessed between HNSCC and HD (Figures 9F,G). The dot-plot heatmap implies that the genes mentioned are more enriched in monocytes than in other kinds of immune cells (Figure 9H).

FIGURE 9

FIGURE 9. (A) Twenty cell clusters were successfully classified by t-SNE algorithm. (B) Expression of several cell surface and intracellular markers in each cell subpopulation. (C,D) Individual cell AUC score overlay for FGF receptor signaling and highly expressed DEGs in low-survival patients (cells from n = 18 HNSCC, n = 5 HD). (E) Composition of clusters 4 and 11 in each sample. (F) Comparison of cluster 4 between HNSCC and HD. (G) Comparison of cluster 11 between HNSCC and HD. (H) Dot plot heatmap showing AUC1 and AUC2 enriched in monocyte (clusters 4 and 11).

Discussion

Considering the poor therapeutic effect and prognosis of patients with HNSCC, it is necessary to construct an accurate prognostic staging system, which might bring personalized treatment and timely follow-up to them. In this study, a total of 483 HNSCC patients from the TCGA data set were included in our discovery cohort to exploit the prognostic value of FGFR-signal, hypoxia, and glycolysis. The effect of immune status of the tumor microenvironment (TME) was also the focus of our research. We found that profibrotic signaling, hypoxia, and glycolysis were associated with the survival of patients with HNSCC. Furthermore, we constructed a new fibrosis–hypoxia–glycolysis–related prognostic classifier including a six-gene signature for HNSCC patients under the guarantee of an external independent validation cohort. IHC was used to determine the protein expression of these genes in HNSCC tissues. The single-cell analysis of monocyte infiltration also revealed marked differences between HNSCC patients and healthy controls. These findings represent a new insight into the prognosis and tumor immune microenvironment of patients with HNSCC.

Multiple studies suggest that the FGFR signal, hypoxia, and glycolysis play a critical role in the tumorigenesis and progression of HNSCC. On the one hand, it is reported that FGFR1 amplification is a frequent event (Göke et al., 2013) and might act as a candidate prognostic biomarker in primary and metastatic HNSCC (Koole et al., 2016). Rogaratinib, an inhibitor that effectively and selectively inhibits pan-FGFR, presents a broad antitumor activity in the FGFR-overexpressing preclinical HNSCC model, revealing the potentially important role of FGFR in disease development (Grünewald et al., 2019). In the mechanism of drug resistance in head and neck cancer stem cells, the FGFR signal also shows latent vitality (McDermott et al., 2018). In addition, the FGFR signal exerts high correlation with sustaining proliferative signaling, resisting cell death, inducing angiogenesis, and activating invasion by cell migration (Ipenburg et al., 2016). However, FGFR signal–related prognostic study was still lacking in HNSCC. On the other hand, to support rapid and unlimited proliferation, solid cancer cells adopt unique energy metabolism properties, such as anaerobic glycolysis (Pavlova and Thompson, 2016) with the faster rate of ATP production and reducing the generation of reactive oxygen species (ROS) mainly produced by the electron transport chain (ETC) in the mitochondria during respiration (Yamamoto et al., 2017). This rule was also applied to HNSCC; glycolysis occurs in HNSCC cells along with a hypoxia microenvironment (Zhu et al., 2017). In the current analysis, hypoxia and glycolysis were presented as risk factors in HNSCC, which was consistent with previous studies. The FGFR signal, hypoxia, and glycolysis play a synergistic role in HNSCC prognosis. Thus, the fibrosis signal, hypoxia, and glycolysis accompanied with their interaction and its relationship with the development of HNSCC could provide improved special insight about the prognosis.

The results of single-cell analysis show that the infiltration of monocytes in the HNSCC group was higher than the HD group. Fibroblasts could communicate with the tumor cells by secreting cytokines in TME. It is reported that monocyte chemotactic protein (MCP)-1, a kind of cytokine associated with poor long-term survival of HNSCC patients (Ji et al., 2014), could be produced by fibroblasts infiltrating in TME to facilitate the recruitment of monocytes into the local inflammatory tissues and regulate their functions (Kondoh et al., 2019). This provides an explanation for the increased infiltration of monocytes in the tumor group and its close interconnection with the FGFR signal.

Significant roles of the predictive signature genes identified above are reported previously in diversified types of cancers. AREG, a ligand of epidermal growth factor receptor (EGFR), is abnormally expressed in multiple types of cancers, such as pancreatic cancer, implicated in mediating the motility, metastasis, and proliferation of tumor cells (Liu et al., 2021). It is proved that the AREG mRNA levels in cancer cells was significantly correlated with the metastatic phenotype of HNSCC tissues (Zhang et al., 2015). THBS1, known as encoding thrombospondin 1, plays a vital role in angiogenesis and tumor progression, overexpression of which was significantly associated with tumor differentiation (Yang et al., 2020). TGFB1 was reported to induce the expression of THBS1, resulting in stimulating migration of cancer cells and driving the expression of MMP3 (matrix metalloprotease 3) via integrin signaling, conducive to OSCC intrusion (Pal et al., 2016). Though the activation of the p-ERK pathway, SEMA3C promotes cervical cancer growth, which is related to poor prognosis (Liu et al., 2019a). ANO1 encodes a calcium-dependent chloride channel protein and commonly amplifies to facilitate several cancers’ progression, including ovarian (Liu et al., 2019b), prostate (Liu et al., 2012), breast (Britschgi et al., 2013), and head and neck cancers (Filippou et al., 2021). After neoadjuvant chemoradiotherapy, the expression level of IGHG2 increased significantly in rectal cancer, indicating that IGHG2 was originally with a low expression in tumor cells and existed as a protective factor, which was consistent with our prognostic analysis. The protective gene EPHX3, known as epoxide hydrolase 3, whose hypermethylation is responsible for the development of OSCC (Morandi et al., 2017), contributes to predict the survival of HNSCC patients (Bai et al., 2019). However, six signature genes in this study were barely mentioned in the context of combination of the FGFR signal, hypoxia, and glycolysis. Thus, the abovementioned signature genes could provide therapeutic targets and directions for the elucidation of molecular mechanisms in HNSCC.

There were inevitable limitations in this study. The first limitation is that IGHG2 was not detected in the independent external cohort, so we could not validate the effect of it. More independent HNSCC cohorts should be used for the validation of the established prognostic model. Using expression profiles downloaded from publicly available databases, it is difficult for us to ensure that the validation samples including all primary tissue sites of tumors in TCGA. Thus, verification of findings above requires more well-designed, comprehensive, and thorough study.

Conclusion

In conclusion, the status of the FGFR signal, hypoxia, and glycolysis correlate with the prognosis of HNSCC patients. The prognostic model conducted above might provide potential application value for prognosis prediction and individualized treatment.

Data Availability Statement

The data sets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The studies involving human participants were reviewed and approved by the Ethics Committee of the 3rd Xiangya Hospital (No: 21158). The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author Contributions

The study was conceived and designed by QZ. Statistical analyses were performed by QC and XL. Software package was prepared by LC, YZ, and QC. Manuscript was written by QC and LC. QC, LC, and XL. contributed equally to this study. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by grants of the National Natural Science Foundation of China (81700658) and the Hunan Provincial Natural Science Foundation-Outstanding Youth Foundation (2020JJ3058).

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, orclaim 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/fcell.2021.801715/full#supplementary-material

References

Bai, G., Song, J., Yuan, Y., Chen, Z., Tian, Y., Yin, X., et al. (2019). Systematic Analysis of Differentially Methylated Expressed Genes and Site‐specific Methylation as Potential Prognostic Markers in Head and Neck Cancerfic Methylation as Potential Prognostic Markers in Head and Neck Cancer. J. Cel Physiol 234 (12), 22687–22702. doi:10.1002/jcp.28835

CrossRef Full Text | Google Scholar

Britschgi, A., Bill, A., Brinkhaus, H., Rothwell, C., Clay, I., Duss, S., et al. (2013). Calcium-activated Chloride Channel ANO1 Promotes Breast Cancer Progression by Activating EGFR and CAMK Signaling. Proc. Natl. Acad. Sci. USA 110 (11), E1026–E1034. doi:10.1073/pnas.1217072110

PubMed Abstract | CrossRef Full Text | Google Scholar

Chae, Y. K., Ranganath, K., Hammerman, P. S., Vaklavas, C., Mohindra, N., Kalyan, A., et al. (2017). Inhibition of the Fibroblast Growth Factor Receptor (FGFR) Pathway: the Current Landscape and Barriers to Clinical Application. Oncotarget 8 (9), 16052–16074. doi:10.18632/oncotarget.14109

PubMed Abstract | CrossRef Full Text | Google Scholar

Chakraborty, D., Zhu, H., Jüngel, A., Summa, L., Li, Y.-N., Matei, A.-E., et al. (2020). Fibroblast Growth Factor Receptor 3 Activates a Network of Profibrotic Signaling Pathways to Promote Fibrosis in Systemic Sclerosis. Sci. Transl. Med. 12 (563), eaaz5506. doi:10.1126/scitranslmed.aaz5506

PubMed Abstract | CrossRef Full Text | Google Scholar

Charoentong, P., Finotello, F., Angelova, M., Mayer, C., Efremova, M., Rieder, D., et al. (2017). Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cel Rep. 18 (1), 248–262. doi:10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Cillo, A. R., Kürten, C. H. L., Tabib, T., Qi, Z., Onkar, S., Wang, T., et al. (2020). Immune Landscape of Viral- and Carcinogen-Driven Head and Neck Cancer. Immunity 52 (1), 183–199. doi:10.1016/j.immuni.2019.11.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Corridoni, D., Antanaviciute, A., Gupta, T., Fawkner-Corbett, D., Aulicino, A., Jagielowicz, M., et al. (2020). Single-cell Atlas of Colonic CD8+ T Cells in Ulcerative Colitis. Nat. Med. 26 (9), 1480–1490. doi:10.1038/s41591-020-1003-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferlay, J., Soerjomataram, I., Dikshit, R., Eser, S., Mathers, C., Rebelo, M., et al. (2015). Cancer Incidence and Mortality Worldwide: Sources, Methods and Major Patterns in GLOBOCAN 2012. Int. J. Cancer 136 (5), E359–E386. doi:10.1002/ijc.29210

CrossRef Full Text | Google Scholar

Filippou, A., Pehkonen, H., Karhemo, P.-R., Väänänen, J., Nieminen, A. I., Klefström, J., et al. (2021). ANO1 Expression Orchestrates p27Kip1/MCL1-Mediated Signaling in Head and Neck Squamous Cell Carcinoma. Cancers 13 (5), 1170. doi:10.3390/cancers13051170

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedman, J., Hastie, T., and Tibshirani, R. (2010). Regularization Paths for Generalized Linear Models via Coordinate Descent. J. Stat. Softw. 33 (1), 1–22. doi:10.18637/jss.v033.i01

PubMed Abstract | CrossRef Full Text | Google Scholar

Göke, F., Bode, M., Franzen, A., Kirsten, R., Goltz, D., Göke, A., et al. (2013). Fibroblast Growth Factor Receptor 1 Amplification Is a Common Event in Squamous Cell Carcinoma of the Head and Neck. Mod. Pathol. 26 (10), 1298–1306. doi:10.1038/modpathol.2013.58

PubMed Abstract | CrossRef Full Text | Google Scholar

Grünewald, S., Politz, O., Bender, S., Héroult, M., Lustig, K., Thuss, U., et al. (2019). Rogaratinib: A Potent and Selective pan‐FGFR Inhibitor with Broad Antitumor Activity in FGFR‐overexpressing Preclinical Cancer Models. Int. J. Cancer 145 (5), 1346–1357. doi:10.1002/ijc.32224

CrossRef Full Text | Google Scholar

Guan, K., Liu, X., Li, J., Ding, Y., Li, J., Cui, G., et al. (2020). Expression Status and Prognostic Value of M6A-Associated Genes in Gastric Cancer. J. Cancer 11 (10), 3027–3040. doi:10.7150/jca.40866

PubMed Abstract | CrossRef Full Text | Google Scholar

Ipenburg, N. A., Koole, K., Liem, K. S., van Kempen, P. M. W., Koole, R., van Diest, P. J., et al. (2016). Fibroblast Growth Factor Receptor Family Members as Prognostic Biomarkers in Head and Neck Squamous Cell Carcinoma: A Systematic Review. Targ Oncol. 11 (1), 17–27. doi:10.1007/s11523-015-0374-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, W.-T., Chen, H.-R., Lin, C.-H., Lee, J.-W., and Lee, C.-C. (2014). Monocyte Chemotactic Protein 1 (MCP-1) Modulates Pro-survival Signaling to Promote Progression of Head and Neck Squamous Cell Carcinoma. PLoS One 9 (2), e88952. doi:10.1371/journal.pone.0088952

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, H. A. J., Zeng, P. Y. F., Shaikh, M. H., Mundi, N., Ghasemi, F., Di Gravio, E., et al. (2021). All HPV-Negative Head and Neck Cancers Are Not the Same: Analysis of the TCGA Dataset Reveals that Anatomical Sites Have Distinct Mutation, Transcriptome, Hypoxia, and Tumor Microenvironment Profiles. Oral Oncol. 116, 105260. doi:10.1016/j.oraloncology.2021.105260

PubMed Abstract | CrossRef Full Text | Google Scholar

Kobak, D., and Berens, P. (2019). The Art of Using T-SNE for Single-Cell Transcriptomics. Nat. Commun. 10 (1), 5416. doi:10.1038/s41467-019-13056-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kondoh, N., Mizuno-Kamiya, M., Umemura, N., Takayama, E., Kawaki, H., Mitsudo, K., et al. (2019). Immunomodulatory Aspects in the Progression and Treatment of Oral Malignancy. Jpn. Dental Sci. Rev. 55 (1), 113–120. doi:10.1016/j.jdsr.2019.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Koole, K., Brunen, D., van Kempen, P. M. W., Noorlag, R., de Bree, R., Lieftink, C., et al. (2016). FGFR1 Is a Potential Prognostic Biomarker and Therapeutic Target in Head and Neck Squamous Cell Carcinoma. Clin. Cancer Res. 22 (15), 3884–3893. doi:10.1158/1078-0432.Ccr-15-1874

PubMed Abstract | CrossRef Full Text | Google Scholar

Korsunsky, I., Millard, N., Fan, J., Slowikowski, K., Zhang, F., Wei, K., et al. (2019). Fast, Sensitive and Accurate Integration of Single-Cell Data with Harmony. Nat. Methods 16 (12), 1289–1296. doi:10.1038/s41592-019-0619-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, B., Fu, T., He, P., Du, C., and Xu, K. (2021). Construction of a Five-Gene Prognostic Model Based on Immune-Related Genes for the Prediction of Survival in Pancreatic Cancer. Biosci. Rep. 41 (7), BSR20204301. doi:10.1042/bsr20204301

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, R., Shuai, Y., Luo, J., and Zhang, Z. (2019a). SEMA3C Promotes Cervical Cancer Growth and Is Associated with Poor Prognosis. Front. Oncol. 9, 1035. doi:10.3389/fonc.2019.01035

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, W., Lu, M., Liu, B., Huang, Y., and Wang, K. (2012). Inhibition of Ca2+-Activated Cl− Channel ANO1/TMEM16A Expression Suppresses Tumor Growth and Invasiveness in Human Prostate Carcinoma. Cancer Lett. 326 (1), 41–51. doi:10.1016/j.canlet.2012.07.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X.-Y., Liang, Y., Xu, Z.-B., Zhang, H., and Leung, K.-S. (2013). AdaptiveL1/2Shooting Regularization Method for Survival Analysis Using Gene Expression Data. Scientific World J. 2013, 1–5. doi:10.1155/2013/475702

CrossRef Full Text | Google Scholar

Liu, Z., Zhang, S., Hou, F., Zhang, C., Gao, J., and Wang, K. (2019b). Inhibition of Ca2+-Activated Chloride Channel ANO1 Suppresses Ovarian Cancer through Inactivating PI3K/Akt Signaling. Int. J. Cancer 144 (9), 2215–2226. doi:10.1002/ijc.31887

CrossRef Full Text | Google Scholar

McDermott, S. C., Rodriguez-Ramirez, C., McDermott, S. P., Wicha, M. S., and Nör, J. E. (2018). FGFR Signaling Regulates Resistance of Head and Neck Cancer Stem Cells to Cisplatin. Oncotarget 9 (38), 25148–25165. doi:10.18632/oncotarget.25358

PubMed Abstract | CrossRef Full Text | Google Scholar

Morandi, L., Gissi, D., Tarsitano, A., Asioli, S., Gabusi, A., Marchetti, C., et al. (2017). CpG Location and Methylation Level Are Crucial Factors for the Early Detection of Oral Squamous Cell Carcinoma in Brushing Samples Using Bisulfite Sequencing of a 13-gene Panel. Clin. Epigenet 9, 85. doi:10.1186/s13148-017-0386-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Pal, S. K., Nguyen, C. T. K., Morita, K.-i., Miki, Y., Kayamori, K., Yamaguchi, A., et al. (2016). THBS1 Is Induced by TGFB1 in the Cancer Stroma and Promotes Invasion of Oral Squamous Cell Carcinoma. J. Oral Pathol. Med. 45 (10), 730–739. doi:10.1111/jop.12430

CrossRef Full Text | Google Scholar

Pavlova, N. N., and Thompson, C. B. (2016). The Emerging Hallmarks of Cancer Metabolism. Cel Metab. 23 (1), 27–47. doi:10.1016/j.cmet.2015.12.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Rankin, E. B., and Giaccia, A. J. (2016). Hypoxic Control of Metastasis. Science 352 (6282), 175–180. doi:10.1126/science.aaf4405

PubMed Abstract | CrossRef Full Text | Google Scholar

Schuler, M., Cho, B. C., Sayehli, C. M., Navarro, A., Soo, R. A., Richly, H., et al. (2019). Rogaratinib in Patients with Advanced Cancers Selected by FGFR mRNA Expression: a Phase 1 Dose-Escalation and Dose-Expansion Study. Lancet Oncol. 20 (10), 1454–1466. doi:10.1016/s1470-2045(19)30412-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Seitz, T., and Hellerbrand, C. (2021). Role of Fibroblast Growth Factor Signalling in Hepatic Fibrosis. Liver Int. 41 (6), 1201–1215. doi:10.1111/liv.14863

PubMed Abstract | CrossRef Full Text | Google Scholar

Shield, K. D., Ferlay, J., Jemal, A., Sankaranarayanan, R., Chaturvedi, A. K., Bray, F., et al. (2017). The Global Incidence of Lip, Oral Cavity, and Pharyngeal Cancers by Subsite in 2012. CA: A Cancer J. Clinicians 67 (1), 51–64. doi:10.3322/caac.21384

CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., and Jemal, A. (2019). Cancer Statistics, 2019. CA A. Cancer J. Clin. 69 (1), 7–34. doi:10.3322/caac.21551

CrossRef Full Text | Google Scholar

Uhlen, M., Oksvold, P., Fagerberg, L., Lundberg, E., Jonasson, K., Forsberg, M., et al. (2010). Towards a Knowledge-Based Human Protein Atlas. Nat. Biotechnol. 28 (12), 1248–1250. doi:10.1038/nbt1210-1248

PubMed Abstract | CrossRef Full Text | Google Scholar

Varghese, F., Bukhari, A. B., Malhotra, R., and De, A. (2014). IHC Profiler: an Open Source Plugin for the Quantitative Evaluation and Automated Scoring of Immunohistochemistry Images of Human Tissue Samples. PLoS One 9 (5), e96801. doi:10.1371/journal.pone.0096801

PubMed Abstract | CrossRef Full Text | Google Scholar

Wollin, L., Wex, E., Pautsch, A., Schnapp, G., Hostettler, K. E., Stowasser, S., et al. (2015). Mode of Action of Nintedanib in the Treatment of Idiopathic Pulmonary Fibrosis. Eur. Respir. J. 45 (5), 1434–1445. doi:10.1183/09031936.00174914

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamamoto, M., Inohara, H., and Nakagawa, T. (2017). Targeting Metabolic Pathways for Head and Neck Cancers Therapeutics. Cancer Metastasis Rev. 36 (3), 503–514. doi:10.1007/s10555-017-9691-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X., Chen, L., Mao, Y., Hu, Z., and He, M. (2020). Progressive and Prognostic Performance of an Extracellular Matrix-Receptor Interaction Signature in Gastric Cancer. Dis. Markers 2020, 1–23. doi:10.1155/2020/8816070

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Wang, Y., Chen, X., Zhou, Y., Jiang, F., Chen, J., et al. (2015). MiR-34a Suppresses Amphiregulin and Tumor Metastatic Potential of Head and Neck Squamous Cell Carcinoma (HNSCC). Oncotarget 6 (10), 7454–7469. doi:10.18632/oncotarget.3148

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Zhou, B., Pache, L., Chang, M., Khodabakhshi, A. H., Tanaseichuk, O., et al. (2019). Metascape Provides a Biologist-Oriented Resource for the Analysis of Systems-Level Datasets. Nat. Commun. 10 (1), 1523. doi:10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, G., Peng, F., Gong, W., She, L., Wei, M., Tan, H., et al. (2017). Hypoxia Promotes Migration/invasion and Glycolysis in Head and Neck Squamous Cell Carcinoma via an HIF-1α-MTDH Loop. Oncol. Rep. 38 (5), 2893–2900. doi:10.3892/or.2017.5949

PubMed Abstract | CrossRef Full Text | Google Scholar