A Neural-Related Gene Risk Score for Head and Neck Squamous Cell Carcinoma

doi:10.21203/rs.3.rs-1735098/v1

Download PDF

Research Article

A Neural-Related Gene Risk Score for Head and Neck Squamous Cell Carcinoma

https://doi.org/10.21203/rs.3.rs-1735098/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Objectives

This study was aimed to establish a neural-related gene risk score (NRGRS) for the prediction of head and neck squamous cell carcinoma (HNSCC) prognosis, and explore its prognostic value on the benefit of immune checkpoint inhibitor (ICI) therapy.

Materials and Methods

Based on the transcriptome data of HNSCC patients (n = 546) from The Cancer Genome Atlas (TCGA), 37 neural-related hub genes were identified by weighted gene co-expression network analysis. Four genes (ITGA5, PYGM, GNG7 and ATP2A3) were identified to construct NRGRS using Lasso-Cox regression method based on the TCGA cohort and validated in the Gene Expression Omnibus (GEO) cohort (n = 109). The survival analysis and Cox regression were performed to validate the prognostic value of NRGRS. Afterward, the molecular and immune characteristics in NRGRS-defined subgroups were analyzed and the prediction value of NRGRS on the benefit of ICI therapy was compared with other biomarkers.

Results

NRGRS-high patients had a worse overall survival than NRGRS-low patients. A low NRGRS was correlated with DNA repair-related and immune response-related pathways, while a high NRGRS was associated with Wnt-beta-catenin signaling, neurotrophic signaling and metabolism-related pathways. Tumors with high NRGRS were more likely to have high TP53, CDKN2A, FAT1 and NOTCH1 mutation rates and have high infiltration of naive CD4 + T cells, resting NK cells, M0, M2 macrophages and resting mast cells, which illustrated suppressive immunity and less benefit from ICI therapy.

Conclusion

In conclusion, NRGRS is a promising prognostic biomarker to predict the survival, the molecular and immune characteristics, and also immunotherapy benefit for HNSCC patients.

Clinical Relevance

This study provides evidence for the potential correlation between neural-related transcriptome alteration and immune activity in HNSCC.

head and neck cancer

prognostic signature

nerve-cancer crosstalk

overall survival

immunotherapy

Head and neck squamous cell carcinoma (HNSCC), the most common type of head and neck cancer, has been widely recognized as a neurotrophic tumor with detection of nerve involvement in around 80% of cases [42]. Perineural invasion (PNI) and high intra-tumoral nerve density are associated with poor clinical outcomes in HNSCC [40]. Within the tumor microenvironment (TME), the reciprocal interaction between tumor cells and the peripheral nervous system leads to not only the cancer cells disseminating along nerves but also the neurogenesis and reprogramming of nerve fibers within the tumor. Neural signaling from intra-tumoral nerve terminals facilitates the aggressiveness of tumor, including tumorigenesis, angiogenesis, invasion and metastasis [23, 40]. Tumor cells secrete some neurotrophic factors [18] and axon guidance regulators [15], which stimulate neurogenesis. A recent study has identified that HNSCC with loss of TP53 contributes to the trans-differentiation of tumor-infiltrated sensory neurons to an adrenergic phenotype, which in turn promotes tumor growth [1]. Growing evidence indicates that the presence of nerve fibers within the tumor is an independent prognostic factor for survival, but whether the expression of neural-related genes is correlated with HNSCC patient prognosis remains largely unknown. There is a need to bridge this gap and develop a novel prognostic model based on neural-related genes.

The specific mechanisms underlying the negative role of nerve-cancer crosstalk in tumor initiation and progression are still poorly understood. A possible explanation for neural infiltration leading to poor prognosis in HNSCC might be the immunosuppressive microenvironment induced by neural signaling. However, there is limited evidence that nerve infiltration or PNI alters immune activity within TME. Thus, it is worthwhile to explore whether neural-related gene transcriptome alteration will affect immune characteristics in HNSCC. Immune checkpoint inhibitor (ICI) therapies targeting programmed death-1 (PD-1)/programmed death-ligand 1 (PD-L1) and cytotoxic T-lymphocyte -associated protein 4 (CTLA4), have been shown to improve the survival of cancer patients with a variety of malignant tumors [37], but only around one-fifth of HNSCC patients respond to ICI due to inadequate immune activation [43]. Thus, the role of neural-related genes in estimating immunotherapy benefits in HNSCC should also be investigated, as they may provide insight into the immunosuppression condition of neural-infiltrated tumors, as well as a promising biomarker for predicting clinical outcomes in HNSCC.

Therefore, this study aimed to construct a neural-related gene risk score (NRGRS) based on neural-related hub genes, and explore the prognostic ability of NRGRS in survival as well as the efficacy of ICI therapy in HNSCC patients, providing evidence for the relationship between neural-related signaling and immune activity in HNSCC.

2.1 Patients and datasets

As the flow chart shown in Fig. 1, RNA sequencing (RNA-seq) data of 546 HNSCC samples, including 502 cancer samples and 44 para-cancer samples, and their clinicopathologic information (one tumor sample was excluded since the follow-up time was 0 day) were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/projects/TCGA-HNSC). The lists of neural-related genes were extracted and downloaded from the QIAGEN Ingenuity Pathway Analysis (IPA) (Trial version) databases.

2.2 Identification of neural-related hub genes

Based on the RNA-seq data of HNSCC samples (tumors vs. normal samples) obtained from TCGA, lists of differentially expressed genes (DEGs) (P < 0.05, | log₂(folder change) | > 1.5, false discovery rate (FDR) < 0.01) were identified by using the edgeR package [38]. After intersecting these genes with the neural-related gene lists, the neural-related DEGs were obtained and analyzed their enrichment in IPA neural-related pathways. Both Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were conducted using Annotation, Visualization, and Integrated Discovery databank (DAVID, version 6.8, http://david.abcc.ncifcrf.gov/home.jsp) [12, 13] to functionally annotate the neural-related DEGs.

Weighted gene co-expression network analysis (WGCNA) was performed to identify neural-related hub genes [51]. The similarity matrix was created by computing the Pearson correlation coefficient between the expression data of each two genes, and it was transformed to an adjacency matrix with a soft threshold of β = 5 (Supplementary Fig. 2A, 2B). Then, a topological matrix was obtained, and the dynamic pruning tree was constructed to identify the modules (Supplementary Fig. 2C). Finally, based on the genes from significantly related modules (the blue and yellow modules), GO and KEGG pathway enrichment analyses were conducted (Supplementary Fig. 2D, 2E, 2F). Meanwhile, the genes in the network with both Gene Significance > 0.2 and Gene Module Membership > 0.6 were identified as the hub genes.

2.3 Construction and validation of the NRGRS

The survival-associated genes were screened and obtained by univariate Cox regression analysis among neural-related hub genes. Then these genes were used to develop a possible risk score using the Lasso-Cox regression algorithm [47]. The genes and their constants for construction of the risk score were defined by minimum criteria, and the perfect penalty parameter λ related to the minimum 10-fold cross validation was chosen. The risk score was calculated through the following formula:

NRGRS = \(\sum _{i=1}^{n}(Coe{f}_{i}\times {x}_{i})\),

in which Coef_i is the coefficient, and x_i is the gene expression value of each selected gene. This formula was established based on the TCGA cohort, and then validated in the GEO cohort (GSE27020) consisting of 109 primary laryngeal carcinoma patients [8]. The patients were stratified into high-NRGRS and low-NRGRS groups based on the median value of the NRGRS in each cohort.

The predictive power of the NRGRS was evaluated by Kaplan–Meier (K–M) survival curves with log-rank tests and the receiver operating characteristic (ROC) curve in both the TCGA and GEO cohorts. To validate the independent prognostic value of NRGRS, univariate and multivariate Cox regression analyses were performed in the TCGA cohort.

Comprehensive analysis of molecular and immune characteristics and ICI therapy in different NRGRS subgroups

The gene set enrichment analysis (GSEA) method was performed on KEGG and HALLMARK gene sets to determine the signaling pathways enrichment in NRGRS-defined subgroups of the TCGA cohort [45]. In the gene mutation analysis, the quantity and quality of gene mutations were analyzed in two NRGRS subgroups by using the Maftools package of R [30]. Spearman’s correlation analysis was performed between NRGRS score and PD-L1 expression and total mutation burden (TMB).

To identify immune characteristics of the TCGA cohort, the expression data were imported into CIBERSORT (https://cibersort.stanford.edu/) to estimate the relative proportion of 22 types of immune cells [31]. The relative proportions of 22 types of immune cells and clinicopathologic factors were compared between the two NRGRS subgroups.

To explore the prognostic value of NRGRS in patients receiving immunotherapy, the prognostic ability of NRGRS was compared with other prognostic models, such as tumor immune dysfunction and exclusion (TIDE) model [17], tumor inflammation signature (TIS) [2] and immune-related gene prognostic index (IRGPI) [5]. TIDE score was calculated online (http://tide.dfci.harvard.edu/) [17], and TIS and IRGPI score was calculated according to the prognostic formula for the signature genes [2, 5].

2.5 Statistical analysis

Univariate Cox regression was performed to obtain survival-associated genes from neural-related hub genes, and to select survival-associated clinicopathologic factors and those factors with P < 0.1 would be included in further multivariate Cox regression analysis. After NRGRS defined subgrouping, the predictive performance of the risk score for survival was evaluated by time-dependent ROC curves. Survival analysis of NRGRS-defined subgroups was performed by Kaplan-Meier survival analysis with the log-rank test. Student t-test was performed to compare continuous variables between two groups. Categorical data were tested using the Chi-square test. Immune cell infiltration, TIDE, TIS and IRGPI scores between groups were compared by the Mann-Whitney Test. Spearman’s correlation was also estimated to explore the correlation between NRGRS and TMB or PD-L1 expression. A two-sided P < 0.05 was considered significant.

3.1 Neural-related DEGs in HNSCC tumors compared to normal samples

In differential expression analysis, a total of 3073 DEGs were obtained, of which 1463 genes were upregulated, and 1610 genes were downregulated in the tumor samples compared with para-tumor normal samples (Supplementary Fig. 1A). In QIAGEN Ingenuity Pathway Analysis (IPA), there were 49 neural-related pathways and 1971 genes within these pathways (Supplementary Table 1). By intersecting these genes with the lists of DEGs in HNSCC tumors compared to normal samples, 314 neural-related DEGs were obtained, of which 157 genes were upregulated and 157 were downregulated in tumor samples compared with normal samples (Supplementary Fig. 1B). According to IPA, the neural-related DEGs were significantly enriched in synaptogenesis signaling, neuroinflammation, CREB signaling, axonal guidance signaling, β-adrenergic signaling pathways, etc. (Fig. 2A). In GO and KEGG pathway analysis, the neural-related pathways were chemical synaptic transmission, postsynaptic membrane, dendrite and neuron projection, neuroactive ligand-receptor interaction, retrograde endocannabinoid signaling, GABAergic synapse and adrenergic signaling (Fig. 2B, Fig. 2C, details in Supplementary Table 2).

3.2 Neural-related hub genes

Through WGCNA analysis, 37 neural-related hub genes were obtained from the candidate genes (n = 314). According to somatic mutation analysis, most of the neural-related hub genes had missense mutations, and the mutation rates of ERBB4, NRG3 and SCN7A were greater than 3% (Supplementary Fig. 3A). The heatmap displaying the interaction of somatic mutation between hub genes is shown in Supplementary Fig. 3B. The somatic mutation of RORC was co-occurrent with the mutation of ADGRF1, NPY5R, ATP2A1 and ATP2A3 mutation, and the somatic mutation of ITGA5 was co-occurrent with GFRA1, SORBS1 and BMP1 mutations (Supplementary Fig. 3B).

Through univariate Cox regression analysis, there were 11 genes (ATP2A3, GNA14, GNG7, GPR12, ITGA5, NPY1R, NPY5R, NRG3, PYGM, TAC1, TGFBR3) related to survival among all neural-related hub genes (P < 0.15, Supplementary Fig. 4). Eight of these genes (GPR12, NPY5R and TAC1 were excluded since over half of the TCGA patients exhibited extremely low expression of these genes) were selected for Lasso-Cox regression analysis (Fig. 3A, 3B), and the risk score formula was finally developed as

NRGRS =\(0.171 \times ITGA5+0.084 \times PYGM-0.216 \times GNG7-0.119 \times ATP2A3.\)

The survival analysis for the four signature genes is shown in Supplementary Fig. 5. Univariate Cox regression analysis showed that age, gender and NRGRS were significantly associated with the prognosis of HNSCC (Fig. 3C). Select the factors with P < 0.1 for multivariate Cox regression analysis, which confirmed that NRGRS was an independent prognostic factor after adjusting other clinicopathologic factors (Fig. 3C).

3.3 Survival outcomes in NRGRS-defined subgroups

Table 1

Clinical characteristics of the HNSCC patients used in this study
		TCGA cohort	GEO cohort
No. of patients		501	109
Age (median, range)		61 (19–90)	64 (41–88)
Gender (%)
	Female	26.7	4.6
	Male	73.3	95.4
Smoking (%)		61.5	99.1
Alcohol (%)		66.3	53.2
Clinical stage (%)
	Stage I/II or unknown	25.4	27.4
	Stage III	20.4	33
	Stage IV	54.1	39.4
Survival status
	OS/DFS years (median)	1.75	2.92
	Alive (%)	56.3	68.8

Table 2

Clinicopathological factors in NRGRS subgroups of the TCGA cohort
		NRGRS-low subgroup	NRGRS-high subgroup	P-value
No. of patients		251	250
Age (median, range)		61 (24–90)	60 (19–88)	0.105
Gender (%)
	Female	24.8	28.5	0.348
	Male	75.2	71.5	0.348
Smoking (%)		64.8	58.2	0.132
Alcohol (%)		65.2	67.5	0.592
Clinical stage (%)
	Stage I/II or unknown	21.6	29.3	0.130
	Stage III	20.8	20.1
	Stage IV	57.6	50.6
Primary site (%)
	Oral cavity	43.2	59.8	< 0.001
	pharynx	56.8	40.2	< 0.001
Survival status
	OS years (median)	1.87	1.59	0.254
	Alive (%)	66.0	46.6	< 0.001

Taking the median NRGRS as the cut-off value in both the TCGA cohort and GEO cohort, the clinical background information of the patients from two cohorts and the correlation between NRGRS subgroups and clinicopathological factors in the TCGA cohort are shown in Tables 1 and 2, respectively. The NRGRS plot and survival status plot of all patients are shown in Fig. 4A, 4B. There were a significantly high proportion of dead patients in the NRGRS-high subgroup compared with those in the NRGRS-low subgroup in both TCGA (Chi-square test, P < 0.001, Fig. 4A) and GEO cohort (P = 0.003, Fig. 4B). NRGRS-high patients had a worse OS in the TCGA cohort (log-rank test, P < 0.001, Fig. 4C) and worse disease-free survival (DFS) in the GEO cohort (P = 0.003, Fig. 4D) than NRGRS-low patients. The predictive performance of the risk score for survival was evaluated by time-dependent ROC curves. The area under the curve (AUC) reached 0.628 at one year, 0.651 at three years, and 0.589 at five years in the TCGA cohort (Fig. 4E), while the AUC reached 0.773 at one year, 0.665 at three years, and 0.667 at five years in the GEO cohort (Fig. 4F).

3.4 Molecular characteristics of different NRGRS subgroups

GSEA was performed to determine the gene sets enriched in different NRGRS subgroups. The gene sets of the NRGRS-high samples were enriched in Wnt-beta-catenin signaling, neurotrophin signaling and metabolism-related pathways (Fig. 5A), while the gene sets of the NRGRS-low samples were enriched in DNA repair-related pathways and immune response-related pathways (Fig. 5B; P < 0.05, FDR < 0.25; details in Supplementary Table 3).

Next, we analyzed gene mutations to gain further biological insight into the NRGRS-defined subgroups. Across all samples from the TCGA cohort, missense mutations were the most common mutation type, followed by nonsense and frameshift deletions, while TP53 was the most frequent mutant gene (67%), followed by TTN (38%) and FAT1 (22%) (Supplementary Fig. 6A). Furthermore, significantly higher mutation counts were found in the NRGRS-high subgroup than in the NRGRS-low subgroup (student t-test, P = 0.011). The top 20 genes with the highest mutation rates in the NRGRS subgroups were identified (Fig. 5C). TP53 had an extremely high mutation rate in the NRGRS-high group compared to the NRGRS-low group (78% vs. 58%). Besides, the mutation of CDKN2A, FAT1 and NOTCH1 genes were more frequent in the NRGRS-high subgroup, while the mutation of CSMD3, MUC16, PIK3CA and SYNE1 genes were more common in the NRGRS-low subgroup (Fig. 5C). Spearman’s correlation analysis showed that NRGRS was not significantly correlated with TMB (r = − 0.054, P = 0.232, Supplementary Fig. 6B) or PD-L1 expression (r = 0.084, P = 0.061, Supplementary Fig. 6C).

3.5 Immune characteristics and the benefit of ICI therapy of different NRGRS subgroups

To analyze the composition of immune cells in different NRGRS subgroups, the Mann-Whitney test was conducted to compare the distribution of immune cells between NRGRS-high and NRGRS-low subgroups. Naïve B cells, memory B cell, plasma B cell, CD8 + T cells, activated memory CD4 + T cells, T helper, regulatory T (Treg) cell, monocytes and activated mast cells were more abundant in the NRGRS-low subgroup, while naive CD4 + T cells, resting NK cells, M0 and M2 macrophages and resting mast cells were highly infiltrated in the NRGRS-high subgroup (Fig. 6A).

TIDE, TIS and IRGPI scores were used to assess the potential clinical efficacy of immunotherapy in different NRGRS subgroups. A lower TIDE prediction score represented a lower potential for immune evasion and an association with a better outcome and higher response to immunotherapy [17], while higher TIS and IRGPI scores showed more benefit from immunotherapy [2, 5]. In our results, the NRGRS-high subgroup had a higher TIDE score and lower TIS and IRGPI scores than the NRGRS-low subgroup, implying that NRGRS-high patients had a lower likelihood of getting benefits from ICI therapy and a worse prognosis than the NRGRS-low patients (Fig. 6B).

HNSCC is a highly innervated solid tumor with neural detection in most cases [42]. Although it has been widely accepted that the intra-tumoral nerve density is positively correlated with poor clinical outcomes in HNSCC [40], the unresolved issues remain on the potential underlying mechanisms and the target therapeutic approach. This study constructed the risk score NRGRS for evaluation of HNSCC clinical outcome based on neural-related hub genes. The comprehensive analysis illustrated that the patients in the NRGRS-high subgroup had worse OS, more immunosuppressive TME and less benefit from ICI therapy than those in the NRGRS-low subgroup. NRGRS was proved to be a promising prognostic model for HNSCC. Our study provided evidence for the correlation between the alteration of the neural-related gene transcriptome and immune activity as well as immunotherapy efficacy in HNSCC.

4.1 Neural-related signature genes of NRGRS

The prognostic model proposed in the present study was composed of 4 neural-related hub genes, namely, ITGA5, PYGM, GNG7 and ATP2A3. Integrin α5 (ITGA5) primarily binds to integrin β1 to form α5β1-integrin [9], which is highly expressed on the developing axons, the regenerating neurites and Schwann cells within the injured nerves, and performs vital functions in neural development and regeneration [24]. Meanwhile, ITGA5 overexpressed in multiple types of cancer and is essential for aggressive behaviors of tumors, including HNSCC [53]. Muscle glycogen phosphorylase or myophosphorylase (PYGM) is considered as a key enzyme in glycogen degradation existing in the muscular system and astrocytes in the central nervous system [14, 33]. Innervation can cause an increase of PYGM expression on human muscle fibers in vitro [29]. PYGM has also been recognized as a prognostic biomarker in several cancer types since glycogen metabolism may exert an essential effect on cancer progression [20, 52]. In the calculation formula of NRGRS, the coefficients of ITGA5 and PYGM were positive, which suggested that ITGA5 and PYGM were positively related to NRGRS and provided evidence for the possible correlation between innervation and poor prognosis.

Sarcoplasmic/endoplasmic reticulum calcium ATPase 3 (ATP2A3) serves as an intracellular pump to sequester calcium into the endoplasmic reticulum, which also widely expresses on neurons [34]. Downregulation of ATP2A3 contributes to the transaction of uterine cervical preneoplastic lesions to neoplastic keratinocytes [39], while upregulation of ATP2A3 acts antitumor effect in breast cancer [16]. G-protein gamma 7 (GNG7) is involved in neurotransmitter receptor pathways, such as dopamine receptor signaling and β-adrenergic receptor signaling [36, 48, 49]. In HNSCC, a decrease of GNG7 protein expression induced by GNG7 promoter methylation is frequent and associated with advanced T stage [11]. According to the NRGRS formula, ATP2A3 and GNG7 had a negative relationship with NRGRS, suggesting that neural-related signaling involving ATP2A3 and GNG7 may participate in anti-tumor activities. The above four neural-related signature genes were reported to be associated with both the nervous system and cancer development, highlighting the role of the peripheral nervous system in cancer progression.

According to the analysis of the somatic mutation interaction between hub genes, the mutation of some genes was found to be co-occurrent with the neural-related signature genes, including RORC and ERBB4 (co-occurrent with ATP2A3), GFRA1, CHN1, and BMP1 (co-occurrent with ITGA5). Receptor protein-tyrosine kinase ERBB4, glycosylphosphatidylinositol-linked glial cell derived neurotrophic factor (GDNF) receptor alpha 1 (GFRA1) and chimerin 1 (CHN1) not only play well-defined roles in the development and maintenance of the nervous system as well as neural signal transduction, but also perform critical regulatory functions in the tumor microenvironment, tumor invasion, and metastasis [21, 25, 41]. As for retinoic acid-related orphan receptor C (RORC) and bone morphogenetic protein 1 (BMP1), they exert critical impact on tumorigenesis in various types of malignant tumors by interacting with the downstream signaling [3, 35]. These findings act a supplement to the essential role of NRGRS signature genes in nervous system and cancer progression.

4.2 Neural-related signaling in NRGRS-defined subgroups

According to the gene enrichment analysis results, neural-related DEGs were highly enriched in synaptogenesis/axonal guidance signaling and neurotransmitter-related signaling pathways, including adrenergic, GABAergic, and endocannabinoid signaling and opioid signaling. As for the GSEA enrichment result for NRGRS-defined subgroups, neurotrophin signaling was enriched in the NRGRS-high subgroup.

Synaptogenesis and axonal guidance signaling enrichment illustrate the neurogenesis involvement in HNSCC, which is consistent with clinical and research findings. Coculture of primary sensory neurons with the HNSCC cell line increased neuritogenesis of neurons in vitro [6]. The axon guidance regulator semaphorin 7A increased in oral cancer and was associated with larger tumor size, more metastatic regional lymph nodes, more advanced clinical stages and poorer 5-year survival [32]. Neurotrophins participate in the development, differentiation, survival, and function of the nervous system [33]. Nerve growth factor (NGF) and its receptor p75 neurotrophin receptor (p75NTR), brain-derived neurotrophic factor (BDNF) and its receptor tropomyosin receptor kinase (Trk) B were highly expressed in OSCC tumors and linked to tumor invasiveness and unfavorable outcome [34–37].

In terms of neurotransmitter-related signaling, GABA, cannabinoid, and noradrenaline, are released from the axon terminals and activate the related neurotransmitter receptors at the plasma membrane of cells in TME [38]. Glutamate decarboxylase 1 (the key enzyme of GABA synthesis), cannabinoid receptor 2 and β2-adrenergic receptor (β2-AR) significantly upregulated in OSCC and were highly associated with tumor invasion and metastasis [39–42]. Meanwhile, a recent study has identified that mutated/deficient TP53 leads to neurogenesis and an adrenergic-transdifferentiation of tumor-infiltrated nerves in HNSCC [6]. In our study, we found that nearly 4/5 of samples had TP53 mutation in the NRGRS-high group, and the mutation proportion was much higher than that in all samples of the TCGA cohort (78% vs. 67%). These results emphasize the reciprocal relationship between adrenergic signaling and HNSCC development. Combined with the above clinical and research findings on neural signaling in cancer, our results help consolidate the correlation between nerve-cancer crosstalk and poor prognosis in HNSCC, and suggest that targeting neural-related pathways might be a promising therapeutic approach.

4.3 NRGRS affects the immune characteristics and immunotherapy benefit

To verify the conjecture that the immunosuppressive TME induced by neural infiltration might be a potential explanation for the negative correlation between intra-tumor nerve density and prognosis, we explored the composition of immune cells in NRGRS-defined subgroups. Higher proportion of CD8 + T cells, activated memory CD4 + T cells, T helper cells, Tregs and activated mast cells were found in the NRGRS-low subgroup, whereas M0, M2 macrophages and resting mast cells were more abundant in the NRGRS-high subgroup. It has been widely recognized that CD8 + T cells and activated memory CD4 + T cells are indicators for an immunoactive microenvironment and a more favorable prognosis [10]. In HNSCC, M1 macrophages have been proven proinflammatory and antitumoral, while M2 macrophages are considered to act immunosuppressive behaviors and promote tumor growth and invasion [6, 28]. T helper cells function antitumor effects indirectly through secreting interferon-γ (IFN-γ) and tumor necrosis factor-α (TNG-α) as well as activating and recruiting antigen-specific effector cells [22, 26]. In HNSCC, high Tregs and mast cells infiltration seemed related to better survival [19, 27, 28]. Our results support these conclusions from previous studies, which implied that the NRGRS-high subgroup was characterized by active tumor progression and immunosuppression.

Concerning the effect of neural-related transcriptome alteration on immune characteristics in HNSCC, we then hypothesized that NRGRS-based differences in the immune-related TME might reflect different immune benefits from ICI therapy. Therefore, we investigated the relationship between NRGRS and known predictive biomarkers for immunotherapy. The TIDE prediction score is a predictor for ICI response based on two immune escape mechanisms, namely T-cell dysfunction in tumors with abundant cytotoxic T lymphocytes (CTL) and T-cell exclusion in tumors infiltrated with a low level of CTL [17]. In our study, patients in the NRGRS-high subgroup had higher TIDE scores and lower CD8 + T cell infiltration with higher T-cell exclusion scores than those in the NRGRS-low subgroup, suggesting their lower ICI responses might attribute to immune evasion via T-cell exclusion [17]. TIS score is a predictor for PD-1/PD-L1 blockade immunotherapy response consisting of 18 signature genes related to IFN-γ [2], and IRGPI score is a prognostic score based on three immune-related hub genes which also shows predictive value ICI therapy benefits in HNSCC [5]. In our study, NRGRS-low patients had significantly higher TIS and IRGPI scores than NRGRS-low patients, suggesting that NRGRS-high patients were more likely to benefit from ICI therapy and indicated more favorable clinical outcomes. As a result, tumors with higher NRGRS tend to have an immunosuppressive TME, which correlates with poorer immunotherapy responses, thus strengthening the prognostic significance of NRGRS in ICI therapeutic benefits.

Though NRGRS has a good consistency with other immunotherapy prognostic models, we get negative results when examining the correlation of NRGRS with TMB and PD-L1 expression. Recently, TMB has been proposed as an essential biomarker for ICI patients selection in clinical trials in several types of cancer [4]. However, TMB estimation seems impractical since it can be confounded by DNA input, sequencing depth, fixation artifacts, variant allele frequency cut points and genome coverage, etc [44]. Here, we found that the NRGRS had no significant correlation with TMB, suggesting thatTMB could not explain the effect of NRGRS in immunotherapy prognosis and further research was needed to address the possible mechanisms. Generally, PD-L1⁺ tumors are more likely to show responses to PD-1/PD-L1 blockade therapies compared to those PD-L1⁻ tumors [32]. However, it has been reported that higher PD-L1 expression on exosomes suppresses CD8 + effector T cell functions and contributes to HNSCC progression [46]. Moreover, the inconsistent results were also reported by a clinical trial that the expression of PD-L1 failed to affect the HNSCC patients benefiting from nivolumab treatment [7]. Our meta-analysis also showed that PD-L1 detected by immunohistochemistry was not a robust prognostic factor for survival in HNSCC [50]. Therefore, the relationship between PD-L1 expression and immunotherapy response needs further verification.

4.4 Limitations

There were several limitations of this study. Firstly, NRGRS was developed and validated using retrospective data from the TCGA and GEO databases, and its clinical utility awaits further validation in prospective research. Secondly, we found a strong connection between NRGRS and immune activity of HNSCC, but the actual mechanism by which neural-related transcriptome or neural signaling modulates immune-related TME is still unclear and needs to be addressed experimentally. Also, the direct evidence for the corelations between the NRGRS genes’ alterations and neural tumor infiltration is lacking. Thirdly, the predictive role of NRGRS on immunotherapy benefits was determined by compassion with other prognostic biomarkers rather than clinical cohorts, and its clinical significance requires additional validation.

In this study, we proposed NRGRS as a prognostic model and explored the molecular and immune characteristics in the NRGRS-defined subgroups. NRGRS is the first prognostic biomarker constructed based on neural-related genes and proved to be independently associated with patients’ survival in both the derivation and validation cohorts, providing insight into the prediction of HNSCC prognosis and evidence for nerve infiltration-induced poor clinical outcome. Additionally, it is the first time to investigate the potential association between neural-related gene transcriptome alteration and immune activity in HNSCC, providing a possible explanation for the poor prognosis of innervated or PNI + tumours and a promising direction for the design of adjuvant therapeutics based on neural-related pathways for immunotherapy. Looking forward, studies elucidating the specific mechanisms that how intra-tumoral nerves remodel antitumor immunity are encouraged. This will pave the way for developing new medications or strategies targeting neural signaling for the management of HNSCC, especially helping improve patients’ benefits from immunotherapy.

Acknowledgments

This study was supported by the Hong Kong Research Grant Council General Research Fund (No. 17120718) and HKU seed fund for basic research (No. 201910159037).

Conflict of Interests

The authors declare that they have no conflict of interest.

Author Contribution

All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Zhuo-ying TAO, Wei-fa Yang, Wang-yong ZHU and Lei-lei WANG. The first draft of the manuscript was written by Zhuo-ying TAO and revised by Yu-xiong SU and Xin-yuan GUAN. Kar Yan Li provided statistical consultation. All authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Funding

This study was supported by the Hong Kong Research Grant Council General Research Fund (No. 17120718) and HKU seed fund for basic research (No. 201910159037).

Ethics Approval and Consent to Participate

Not Applicable.

Amit M, Takahashi H, Dragomir MP, Lindemann A, Gleber-Netto FO, Pickering CR, Anfossi S, Osman AA, Cai Y, Wang R, Knutsen E, Shimizu M, Ivan C, Rao X, Wang J, Silverman DA, Tam S, Zhao M, Caulin C, Zinger A, Tasciotti E, Dougherty PM, El-Naggar A, Calin GA, Myers JN (2020) Loss of p53 drives neuron reprogramming in head and neck cancer. Nature 578:449–454. http://doi.org/10.1038/s41586-020-1996-3
Ayers M, Lunceford J, Nebozhyn M, Murphy E, Loboda A, Kaufman DR, Albright A, Cheng JD, Kang SP, Shankaran V, Piha-Paul SA, Yearley J, Seiwert TY, Ribas A, Mcclanahan TK (2017) IFN-gamma-related mRNA profile predicts clinical response to PD-1 blockade. J Clin Invest 127:2930–2940. http://doi.org/10.1172/JCI91190
Cao D, Qi Z, Pang Y, Li H, Xie H, Wu J, Huang Y, Zhu Y, Shen Y, Zhu Y, Dai B, Hu X, Ye D, Wang Z (2019) Retinoic Acid-Related Orphan Receptor C Regulates Proliferation, Glycolysis, and Chemoresistance via the PD-L1/ITGB6/STAT3 Signaling Axis in Bladder Cancer. Cancer Res 79:2604–2618. http://doi.org/10.1158/0008-5472.CAN-18-3842
Chan TA, Yarchoan M, Jaffee E, Swanton C, Quezada SA, Stenzinger A, Peters S (2019) Development of tumor mutation burden as an immunotherapy biomarker: utility for the oncology clinic. Annals of Oncology 30:44–56. http://doi.org/10.1093/annonc/mdy495
Chen Y, Li ZY, Zhou GQ, Sun Y (2021) An Immune-Related Gene Prognostic Index for Head and Neck Squamous Cell Carcinoma. Clin Cancer Res 27:330–341. http://doi.org/10.1158/1078-0432.CCR-20-2166
Evrard D, Szturz P, Tijeras-Raballand A, Astorgues-Xerri L, Abitbol C, Paradis V, Raymond E, Albert S, Barry B, Faivre S (2019) Macrophages in the microenvironment of head and neck cancer: potential targets for cancer therapy. Oral Oncology 88:29–38. http://doi.org/10.1016/j.oraloncology.2018.10.040
Ferris RL, Blumenschein G, Fayette J, Guigay J, Colevas AD, Licitra L, Harrington KJ, Kasper S, Vokes EE, Even C, Worden F, Saba NF, Docampo LCI, Haddad R, Rordorf T, Kiyota N, Tahara M, Lynch M, Jayaprakash V, Li L, Gillison ML (2018) Nivolumab vs investigator’s choice in recurrent or metastatic squamous cell carcinoma of the head and neck: 2-year long-term survival update of CheckMate 141 with analyses by tumor PD-L1 expression. Oral Oncology 81:45–51. http://doi.org/10.1016/j.oraloncology.2018.04.008
Fountzilas E, Kotoula V, Angouridakis N, Karasmanis I, Wirtz RM, Eleftheraki AG, Veltrup E, Markou K, Nikolaou A, Pectasides D, Fountzilas G (2013) Identification and Validation of a Multigene Predictor of Recurrence in Primary Laryngeal Cancer. PLoS ONE 8:e70429. http://doi.org/10.1371/journal.pone.0070429
Fujita M, Zhu K, Fujita CK, Zhao M, Lam KS, Kurth MJ, Takada YK, Takada Y (2015) Proinflammatory Secreted Phospholipase A2 Type IIA (sPLA-IIA) Induces Integrin Activation through Direct Binding to a Newly Identified Binding Site (Site 2) in Integrins αvβ3, α4β1, and α5β1. Journal of Biological Chemistry 290:259–271. http://doi.org/10.1074/jbc.m114.579946
Gentles AJ, Newman AM, Liu CL, Bratman SV, Feng W, Kim D, Nair VS, Xu Y, Khuong A, Hoang CD, Diehn M, West RB, Plevritis SK, Alizadeh AA (2015) The prognostic landscape of genes and infiltrating immune cells across human cancers. Nature Medicine 21:938–945. http://doi.org/10.1038/nm.3909
Hartmann S, Szaumkessel M, Salaverria I, Simon R, Sauter G, Kiwerska K, Gawecki W, Bodnar M, Marszalek A, Richter J, Brauze D, Zemke N, Jarmuz M, Hansmann M-L, Siebert R, Szyfter K, Giefing M (2012) Loss of protein expression and recurrent DNA hypermethylation of the GNG7 gene in squamous cell carcinoma of the head and neck. Journal of Applied Genetics 53:167–174. http://doi.org/10.1007/s13353-011-0079-4
Huang Da W, Sherman BT, Lempicki RA (2009) Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res 37:1–13. http://doi.org/10.1093/nar/gkn923
Huang Da W, Sherman BT, Lempicki RA (2009) Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 4:44–57. http://doi.org/10.1038/nprot.2008.211
Ignacio PC, Baldwin BA, Vijayan VK, Tait RC, Gorin FA (1990) Brain isozyme of glycogen phosphorylase: immunohistologicall localization within the central nervous system. Brain Research 529:42–49. http://doi.org/10.1016/0006-8993(90)90809-p
Iyer A, Chapoval S (2018) Neuroimmune Semaphorin 4A in Cancer Angiogenesis and Inflammation: A Promoter or a Suppressor? International Journal of Molecular Sciences 20:124. http://doi.org/10.3390/ijms20010124
Izquierdo-Torres E, Rodriguez G, Meneses-Morales I, Zarain-Herzberg A (2017) ATP2A3 gene as an important player for resveratrol anticancer activity in breast cancer cells. Mol Carcinog 56:1703–1711. http://doi.org/10.1002/mc.22625
Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, Li Z, Traugh N, Bu X, Li B, Liu J, Freeman GJ, Brown MA, Wucherpfennig KW, Liu XS (2018) Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nature Medicine 24:1550–1558. http://doi.org/10.1038/s41591-018-0136-1
Jimenez-Andrade JM, Ghilardi JR, Castaneda-Corral G, Kuskowski MA, Mantyh PW (2011) Preventive or late administration of anti-NGF therapy attenuates tumor-induced nerve sprouting, neuroma formation, and cancer pain. Pain 152:2564–2574. http://doi.org/10.1016/j.pain.2011.07.020
Jin Y, Qin X (2020) Profiles of immune cell infiltration and their clinical significance in head and neck squamous cell carcinoma. International Immunopharmacology 82:106364. http://doi.org/10.1016/j.intimp.2020.106364
Jin Y, Yang Y (2019) Bioinformatics-based discovery of PYGM and TNNC2 as potential biomarkers of head and neck squamous cell carcinoma. Bioscience Reports 39:BSR20191612. http://doi.org/10.1042/bsr20191612
Kim M, Kim DJ (2018) GFRA1: A Novel Molecular Target for the Prevention of Osteosarcoma Chemoresistance. Int J Mol Sci 19http://doi.org/10.3390/ijms19041078
Knutson KL, Disis ML (2005) Tumor antigen-specific T helper cells in cancer immunity and immunotherapy. Cancer Immunology, Immunotherapy 54:721–728. http://doi.org/10.1007/s00262-004-0653-2
Lee TL, Chiu PH, Li WY, Yang MH, Wei PY, Chu PY, Wang YF, Tai SK (2019) Nerve-tumour interaction enhances the aggressiveness of oral squamous cell carcinoma. Clin Otolaryngol 44:1087–1095. http://doi.org/10.1111/coa.13452
Lefcort F, Venstrom K, Mcdonald JA, Reichardt LF (1992) Regulation of expression of fibronectin and its receptor, alpha 5 beta 1, during development and regeneration of peripheral nerve. Development 116:767–782. http://doi.org/10.1242/dev.116.3.767
Li JP, Zeng SH, Zhang YH, Liu YJ (2021) Bioinformatics-based analysis of the association between the A1-chimaerin (CHN1) gene and gastric cancer. Bioengineered 12:2874–2889. http://doi.org/10.1080/21655979.2021.1940621
Li T, Wu B, Yang T, Zhang L, Jin K (2020) The outstanding antitumor capacity of CD4(+) T helper lymphocytes. Biochim Biophys Acta Rev Cancer 1874:188439. http://doi.org/10.1016/j.bbcan.2020.188439
Liang B, Tao Y, Wang T (2020) Profiles of immune cell infiltration in head and neck squamous carcinoma. Bioscience Reports 40:2. http://doi.org/10.1042/bsr20192724
Lukesova E, Boucek J, Rotnaglova E, Salakova M, Koslabova E, Grega M, Eckschlager T, Rihova B, Prochazka B, Klozar J, Tachezy R (2014) High Level of Tregs Is a Positive Prognostic Marker in Patients with HPV-Positive Oral and Oropharyngeal Squamous Cell Carcinomas. BioMed Research International 2014:1–11. http://doi.org/10.1155/2014/303929
Martinuzzi A, Askanas V, Kobayashi T, Engel WK (1988) Asynchronous regulation of muscle specific isozymes of creatine kinase, glycogen phosphorylase, lactic dehydrogenase and phosphoglycerate mutase in innervated and non-innervated cultured human muscle. Neuroscience Letters 89:216–222. http://doi.org/10.1016/0304-3940(88)90384-9
Mayakonda A, Lin D-C, Assenov Y, Plass C, Koeffler HP (2018) Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Research 28:1747–1756. http://doi.org/10.1101/gr.239244.118
Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA (2015) Robust enumeration of cell subsets from tissue expression profiles. Nature Methods 12:453–457. http://doi.org/10.1038/nmeth.3337
Patel SP, Kurzrock R (2015) PD-L1 Expression as a Predictive Biomarker in Cancer Immunotherapy. Molecular Cancer Therapeutics 14:847–856. http://doi.org/10.1158/1535-7163.mct-14-0983
Pfeiffer-Guglielmi B, Francke M, Reichenbach A, Hamprecht B (2007) Glycogen phosphorylase isozymes and energy metabolism in the rat peripheral nervous system—An immunocytochemical study. Brain Research 1136:20–27. http://doi.org/10.1016/j.brainres.2006.12.037
Pottorf WJ, De Leon DD, Hessinger DA, Buchholz JN (2001) Function of SERCA mediated calcium uptake and expression of SERCA3 in cerebral cortex from young and old rats. Brain Res 914:57–65. http://doi.org/10.1016/s0006-8993(01)02773-1
Rafi JH, Jafar T, Pathan MT, Reza R, Islam S, Sourna IJ, Alam R, Samad A, Ahammad F (2021) High expression of bone morphogenetic protein 1 (BMP1) is associated with a poor survival rate in human gastric cancer, a dataset approaches. Genomics 113:1141–1154. http://doi.org/10.1016/j.ygeno.2020.11.012
Ray K, Kunsch C, Bonner LM, Robishaw JD (1995) Isolation of cDNA Clones Encoding Eight Different Human G Protein γ Subunits, Including Three Novel Forms Designated the γ4, γ10, and γ11 Subunits. Journal of Biological Chemistry 270:21765–21771. http://doi.org/10.1074/jbc.270.37.21765
Riley RS, June CH, Langer R, Mitchell MJ (2019) Delivery technologies for cancer immunotherapy. Nat Rev Drug Discov 18:175–196. http://doi.org/10.1038/s41573-018-0006-z
Robinson MD, Mccarthy DJ, Smyth GK (2010) edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26:139–140. http://doi.org/10.1093/bioinformatics/btp616
Rotondo JC, Bosi S, Bassi C, Ferracin M, Lanza G, Gafa R, Magri E, Selvatici R, Torresani S, Marci R, Garutti P, Negrini M, Tognon M, Martini F (2015) Gene expression changes in progression of cervical neoplasia revealed by microarray analysis of cervical neoplastic keratinocytes. J Cell Physiol 230:806–812. http://doi.org/10.1002/jcp.24808
Schmitd LB, Scanlon CS, D’silva NJ (2018) Perineural Invasion in Head and Neck Cancer. Journal of Dental Research 97:742–750. http://doi.org/10.1177/0022034518756297
Segers VFM, Dugaucquier L, Feyen E, Shakeri H, De Keulenaer GW (2020) The role of ErbB4 in cancer. Cell Oncol (Dordr) 43:335–352. http://doi.org/10.1007/s13402-020-00499-4
Silverman DA, Martinez VK, Dougherty PM, Myers JN, Calin GA, Amit M (2021) Cancer-Associated Neurogenesis and Nerve-Cancer Cross-talk. Cancer Research 81:1431–1440. http://doi.org/10.1158/0008-5472.can-20-2793
Sim F, Leidner R, Bell RB (2019) Immunotherapy for Head and Neck Cancer. Oral Maxillofac Surg Clin North Am 31:85–100. http://doi.org/10.1016/j.coms.2018.09.002
Stenzinger A, Endris V, Budczies J, Merkelbach-Bruse S, Kazdal D, Dietmaier W, Pfarr N, Siebolts U, Hummel M, Herold S, Andreas J, Zoche M, Tögel L, Rempel E, Maas J, Merino D, Stewart M, Zaoui K, Schlesner M, Glimm H, Fröhling S, Allen J, Horst D, Baretton G, Wickenhauser C, Tiemann M, Evert M, Moch H, Kirchner T, Büttner R, Schirmacher P, Jung A, Haller F, Weichert W, Dietel M (2020) Harmonization and Standardization of Panel-Based Tumor Mutational Burden Measurement: Real-World Results and Recommendations of the Quality in Pathology Study. Journal of Thoracic Oncology 15:1177–1189. http://doi.org/10.1016/j.jtho.2020.01.023
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP (2005) Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proceedings of the National Academy of Sciences 102:15545–15550. http://doi.org/10.1073/pnas.0506580102
Theodoraki M-N, Yerneni SS, Hoffmann TK, Gooding WE, Whiteside TL (2018) Clinical Significance of PD-L1 + Exosomes in Plasma of Head and Neck Cancer Patients. Clinical Cancer Research 24:896–905. http://doi.org/10.1158/1078-0432.ccr-17-2664
Tibshirani R (1997) The lasso method for variable selection in the Cox model. Stat Med 16:385–395. http://doi.org/10.1002/(sici)1097-0258(19970228)16:4<385::aid-sim380>3.0.co;2-3
Wang Q, Jolly JP, Surmeier JD, Mullah BM, Lidow MS, Bergson CM, Robishaw JD (2001) Differential Dependence of the D1 and D5Dopamine Receptors on the G Protein γ7 Subunit for Activation of Adenylylcyclase. Journal of Biological Chemistry 276:39386–39393. http://doi.org/10.1074/jbc.m104981200
Wang Q, Mullah BK, Robishaw JD (1999) Ribozyme Approach Identifies a Functional Association between the G Protein β1γ7 Subunits in the β-Adrenergic Receptor Signaling Pathway. Journal of Biological Chemistry 274:17365–17371. http://doi.org/10.1074/jbc.274.24.17365
Yang X, Tian X, Wu K, Liu W, Li S, Zhang Z, Zhang C (2018) Prognostic impact of perineural invasion in early stage oral tongue squamous cell carcinoma: Results from a prospective randomized trial. Surg Oncol 27:123–128. http://doi.org/10.1016/j.suronc.2018.02.005
Zhang B, Horvath S (2005) A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol 4:Article17. http://doi.org/10.2202/1544-6115.1128
Zhang W, Liu S, Zhan H, Yan Z, Zhang G (2018) Transcriptome sequencing identifies key pathways and genes involved in gastric adenocarcinoma. Molecular Medicine Reportshttp://doi.org/10.3892/mmr.2018.9370
Zhu H, Wang G, Zhu H, Xu A (2021) ITGA5 is a prognostic biomarker and correlated with immune infiltration in gastrointestinal tumors. BMC Cancer 21http://doi.org/10.1186/s12885-021-07996-1

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

A Neural-Related Gene Risk Score for Head and Neck Squamous Cell Carcinoma

Status:

Version 1

Abstract

Figures

1 Introduction

2 Materials And Methods

3 Results

4 Discussion

Conclusion

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1