Exploring NK‐Cell molecules that impact the immune response and microenvironment in head and neck squamous cell carcinoma

Abstract NK cells play a role in various cancers, but their role in head and neck squamous cell carcinoma (HNSCC) still needs to be explored. All public data are obtained from the Cancer Genome Atlas Program (TCGA) database. All analysis was performed using specific packages in R software. In our study, we quantified the immune microenvironment of HNSCC through multiple algorithms. Next, we identified NK cell‐associated genes by quantifying NK cells, including SSNA1, TRIR, PAXX, DPP7, WDR34, EZR, PHLDA1 and ELOVL1. Then, we explored the single‐cell expression pattern of these genes in the HNSCC microenvironment. Univariate Cox regression analysis indicated that the EZR, PHLDA1 and ELOVL1 were related to the prognosis of HNSCC patients. Following this, we selected EZR for further analysis. Our results showed that the patients with high EZR expression might have a poor prognosis and worse clinical features. Biological enrichment analysis showed that EZR is associated with many oncogenic pathways and a higher tumour stemness index. Meanwhile, we found that EZR can remodel the immune microenvironment of HNSCC. Moreover, we noticed that EZR could affect the immunotherapy and specific drug sensitivity, making it an underlying clinical target. In summary, our results can improve the understanding of NK cell in HNSCC. Meanwhile, we identified EZR as the underlying clinical target of HNSCC.

For advanced or recurrent diseases, comprehensive treatment strategies are more common. 4But even after comprehensive treatment, the five-year survival rate of advanced HNSCC is still not high. 5In recent years, with an in-depth understanding of tumour biology, targeted therapy and immunotherapy have gradually become research hotspots in the treatment of HNSCC. 6Some new therapeutic drugs, such as immune checkpoint inhibitors against PD-1/PD-L1, have shown promise in the treatment of HNSCC, bringing new treatment options for patients. 7 cells, or natural killer cells, are a key component of the innate immune system.Such cells have a specific killing effect on non-self and altered cells in the body, such as virus-infected cells or tumour cells. 8Unlike T cells and B cells, NK cells do not need to be stimulated by specific antigens to display their cell-killing activity. 9In cancer, NK cells play a crucial role.They can identify and destroy early tumour cells, thereby preventing further development and metastasis of tumours. 10However, many established tumours can evade NK cell surveillance and attack through various mechanisms, such as secreting immunosuppressive factors or changing their antigenic properties. 11erefore, tumour immune escape from NK cells is a challenge in cancer therapy. 12In recent years, exploiting the antitumor activity of NK cells as a cancer treatment strategy has attracted considerable attention. 13For example, NK cell immunotherapy, which aims to enhance the ability of NK cells to attack tumours, is gradually becoming a promising treatment method.
In our study, we quantified the immune microenvironment of HNSCC through multiple algorithms.Next, we identified NK cell-associated genes by quantifying NK cells, including SSNA1, TRIR, PAXX, DPP7, WDR34, EZR, PHLDA1 and ELOVL1.Then, we explored the single-cell expression pattern of these genes in the HNSCC microenvironment.Univariate Cox regression analysis indicated that the EZR, PHLDA1 and ELOVL1 were related to the prognosis of HNSCC patients.Following this, we selected EZR for further analysis.Our results showed that the patients with high EZR expression might have a poor prognosis and worse clinical features.Biological relation analysis showed that EZR is associated with many oncogenic pathways and a higher tumour stemness index.Meanwhile, we found that EZR can remodel the immune microenvironment of HNSCC.
Moreover, we noticed that EZR could affect the immunotherapy and specific drug sensitivity, making it an underlying clinical target.

| Screening of NK cell marker genes
Data of RNA sequences and related clinical information of 528 HNSC patients were gathered from the TCGA database (accessed time: 2023/7/11, https:// tcga-data.nci.nih.gov/ tcga/ ), and the high-throughput sequencing data were transformed into TPM values. 14Based on multiple immune infiltration algorithms including CIBERSORT, EPIC, MCPCOUNTER, QUANTISEQ, TIMER and XCELL algorithms from the online tool (www.home-for-resea rchers.[17][18][19][20] After observing the distribution of the data, EPIC, QUANTISEQ and XCELL algorithms were screened for the differential expression analysis, and the co-expressed genes including upregulated and downregulated genes were also retained as the NK cell-related genes.

| Related single-cell analysis and prognostic analysis
First, correlation analysis was performed on these NK cell-related genes including SSNA1, TRIR, PAXX, DPP7, WDR34, EZR, PHLDA1 and ELOVL1.Then, differentially expressed genes (DEGs) analysis and univariate prognostic analysis were also conducted, followed by single-cell analysis of these genes using the online tool (http:// tisch.comp-genom ics.org/ home/ ). 21Lastly, Kaplan-Meier (KM) curves were analysed and differences in clinical parameters including age, gender, grade and stage were also analysed.

| Functional enrichment analysis
Gene set variation analysis (GSVA) was utilized to identify the potential molecular mechanisms between high-and low-expression groups based on the Hallmark gene set from the molecular signature database (https:// www.gsea-msigdb.org/ gsea/ msigdb). 22Gene ontology (GO) enrichment analysis was performed to explore the different biological processes (BP) and molecular functions (MF) of the differently expressed genes between two groups. 23

| Genomic feature analysis
Differential analyses of tumour mutation burden (TMB) and tumour stemness index (mRNAsi) were performed to examine the genomic landscape of molecules involved in two groups.With the input of "SNP6" files, we calculated arm-and focal-level somatic copy number alterations (SCNAs) in the tumour using GISTIC2, which were downloaded from the genomic data commons data portal (https:// portal.gdc.cancer.gov/ ). 24

| Immune cell infiltration analysis
Three scores including the immune score, the stromal score and the estimate immune score calculated based on the "estimate" R package were used for the evaluation of immune cells infiltration.Furthermore, immunocyte infiltration analysis according to the 42-immunocytes project was also performed between these two groups.

| ICB therapy response and drug sensitivity prediction
The underlying ICB (PD-L1, CTLA4, HAVCR2, LAG3 and PD1) response for HNSC patients was predicted with tumour immune dysfunction and exclusion (TIDE) analysis. 25Besides, the submap algorithm was also performed to reveal which group of patients may respond to anti-PD-1 and anti-CTLA4 therapy.The potential chemotherapy response of the chemotherapy drugs for patients was evaluated using the Genomics of Drug Sensitivity in Cancer (GDSC) database. 26

| Statistical analysis
A version of R software 4.2.2 was used to perform all the R packages.Comparisons between two groups were analysed with Wilcoxon rank-sum tests, and continuous variables were analysed using Wilcoxon rank-sum tests.In all cases, p values were set up as two-sided, and p values less than 0.05 were considered statistically significant.

| Identification of NK cell marker genes
The whole flow chart of our study was shown in Figure S1.The  was relatively high (Figure 1B).Differential expression analysis revealed a notable difference in NK cell abundance between tumour and normal samples (Figure 1C).According to the Venn diagram, five common up-regulated genes (Figure 1D) and three common down-regulated genes (Figure 1E) were identified across these three algorithms.

| Relationship between these genes and determination of the hub NK cell marker gene
We first examined the pairwise correlation between genes using Pearson correlation analysis across HNSCC samples (Figure 2A).significance of these genes, and EZR was identified for prognosis in HNSC patients with the lowest p-value (Figure 2B).The differential expression was assessed for each gene and the EZR has the lowest p value (Figure 2C).Following this, the expression distribution of these genes in multiple cell subgroups was analysed both in the GSE103322 cohort and GSE139324 cohort (Figure 2D,E).

| Analysis of clinical characteristics and functional enrichment analysis
After prognostic analysis, EZR served as the key gene strongly associated with the prognosis of HNSCC patients.We found that EZR expressed highly in most tumours (Figure 3A).HNSCC patients were stratified into high-and low-expression groups using the median value of EZR expression as a cutoff.KM survival curves based on overall survival, disease-specific survival and progress-free interval were all performed and a significant difference in overall survival time was observed between the two groups (Figure 3B-D).To unravel the profile of clinical manifestation for these two groups, clinical characteristics analysis was then performed (Figure 3E).
Next, to explore the potential biological functions and signalling pathways between the two groups, GSVA analysis was conducted.
In the high-expression group, xenobiotic metabolism, oxidative phosphorylation, DNA repair and E2F targets were enriched.In the low-expression group, inflammatory response, complement, epithelial-mesenchymal transition (EMT) and protein secretion were enriched (Figure 4A).In terms of BP, regulation of vascular endothelial cell proliferation and regulation of drug response were enriched in the high expression group while the sensory perception of smell was enriched in the low expression group (Figure 4B,C).
Besides, the enriched MF terms in the high expression group were mRNA base-pairing post-transcriptional repressor activity, translation repressor activity and translation regulator activity.The enriched MF terms in the low expression group were olfactory receptor activity, odorant binding, and mRNA base-pairing posttranscriptional repressor activity (Figure 4D,E).

| Comparison of genomic characteristics between groups
The TMB levels in the HNSCC cohort were relatively high in 33 different types of tumours (Figure 5A).However, in the HNSCC cohort, no significant difference in TMB level was detected between the two groups (Figure 5B).Then, genes associated with differential mutants (the number of mutations more than 20) were identified between two groups (Figure 5C), and the co-occurrence and exclusive relationship between these differential mutant genes was shown in Figure 5D.In addition, an overview of the clinical characteristics and mRNAsi distribution for each sample can be found in Figure 5E, and no significant difference in mRNAsi level was observed from differential expression analysis (Figure 5F).Chromosomal aberrations of each HNSC patient including copy number percentage and copy number gistic score were assessed separately (Figure 6A,B).Also, we found that there is no significant difference in amplification and deletion frequencies between the two groups (Figure 6C,D).

| Features of immune infiltration and prediction of ICB therapy
According to ESTIMATE analysis, patients in the low-expression group had higher immune, stromal and estimate scores than those in the high-expression group (Figure 7A-C).Then, scores of immune function and immune cell abundance were calculated using the CIBERSORT and ssGSEA algorithms.The results showed that

F I G U R E 6 Analysis of chromosomal aberrations of between two groups. (A) The copy number percentage of each HNSC patient. (B)
The copy number gistic score of each HNSC patient.(C,D) No significant difference of amplification and deletion frequencies was observed between two groups.immune terms associated with the T-cell were enriched in the highexpression group (Figure 7D,E).Next, we compared the levels of immune checkpoint expression between the two groups and found that high-expression patients were more highly expressed (Figure 7F).
The distribution of TIDE scores was different among HNSCC patients (Figure 8A).The TIDE, disfunction and exclusion scores were lower in the high-expression group (Figure 8B).According to TIDE analysis, patients in the high-expression group may be more sensitive to ICB therapy (Figure 8C).Submap analysis also suggested that patients in the high-expression group would tend to respond to anti-CTLA4 therapy (Figure 8D).Finally, the drug sensitivity analysis revealed that patients in the high expression group may tend to respond to Navitoclax_1011, PLX − 4720_1036, YK − 4 − 279_1239, AZ960_1250, Pevonedistat_1529, Docetaxel_1819, AZD6738_1917 and Telomerase Inhibitor IX_1930 (Figure 8E).

| DISCUSS ION
HNSCC is the most common head and neck tumour with relatively high morbidity and mortality worldwide. 27 especially for patients with advanced or recurrent disease. 28In recent years, targeted therapy has become a new direction in the treatment of HNSCC. 29By targeting specific biomarkers or signalling pathways, targeted drugs aim to reduce damage to normal cells, thereby improving treatment outcomes and patients' quality of life.
However, further research is needed on its long-term effects and safety.
In our study, we quantified the immune microenvironment of HNSCC through multiple algorithms.Next, we identified NK cell-associated genes by quantifying NK cells, including SSNA1, TRIR, PAXX, DPP7, WDR34, EZR, PHLDA1 and ELOVL1.Then, we explored the single-cell expression pattern of these genes in the HNSCC microenvironment.Univariate Cox regression analysis indicated that the EZR, PHLDA1 and ELOVL1 were related to the prognosis of HNSCC patients.Following this, we selected EZR for further analysis.Our results showed that the patients with high EZR expression might have a poor prognosis and worse clinical features.
Biological enrichment analysis showed that EZR is associated with many oncogenic pathways and a higher tumour stemness index.Meanwhile, we found that EZR can remodel the immune microenvironment of HNSCC.Moreover, we noticed that EZR could affect the immunotherapy and specific drug sensitivity, making it an underlying clinical target.
Our results showed that the patients with high EZR expression might have a higher activity of xenobiotic metabolism, oxidative phosphorylation, G2M checkpoint, PI3K/AKT/mTOR signalling, fatty acid metabolism and bile acid metabolism.Oxidative phosphorylation is the main energy-producing process of the cell, involving the generation of ATP. 30 In cancer, the role of oxidative phosphorylation is complex.While many cancer cells prefer glycolysis for energy production due to the Warburg effect, some tumour cells still rely on or enhance oxidative phosphorylation. 31Enhanced oxidative phosphorylation provides cancer cells with the necessary energy and biosynthesis to support their rapid proliferation and growth. 32In specific cases, enhanced oxidative phosphorylation in cancer cells is associated with chemotherapy resistance and cancer stem cell properties.
The G2M checkpoint is a critical stage in the cell cycle responsible for ensuring that DNA has been replicated intact before cells enter mitosis. 33During cancer development, this checkpoint is often disturbed.Cancer cells may bypass the impaired G2M checkpoint, causing them to continue dividing even in the presence of DNA damage, which further leads to the accumulation of genetic mutations. 34The PI3K/AKT/mTOR signalling pathway plays a key role in many biological processes, including cell growth, proliferation, metabolism and survival. 35In cancer, this pathway is often activated to promote tumour formation and progression.Abnormal PI3K activity can lead to overactivation of AKT, which in turn activates mTOR, which promotes protein synthesis, cell proliferation and inhibits apoptosis.In addition, this pathway is also closely related to cancer drug resistance, metastasis and cancer stem cell properties. 36The role of lipid metabolism in cancer has received extensive attention in recent years. 37Cancer cells often reprogram their lipid metabolism pathways to meet their demands for rapid proliferation and growth.This not only provides the cancer cell with a source of energy but also with the necessary components of the cell membrane and signalling molecules. 38Our results provide indications for the mechanism of action of EZR in HNSCC.
Meanwhile, we found that EZR was associated with a higher tumour stemness index and can remodel the tumour microenvironment of HNSCC.Cancer stem cells are a class of cells with self-renewal and multilineage differentiation capabilities and are regarded as key factors in cancer development, recurrence and metastasis. 39ncer stem cells play a crucial role in HNSCC.These cells exhibit strong resistance to traditional treatments such as chemotherapy and radiotherapy, making the treatment of HNSCC a great challenge. 40Therefore, in-depth understanding and targeting of cancer stem cells in HNSCC have become the focus of current research.The tumour microenvironment refers to the environment composed of cells and matrix around tumour cells, including immune cells, fibroblasts, vascular cells, and extracellular matrix. 41This environment plays a critical role in cancer progression, immune evasion, and therapy resistance.In HNSCC, the role of the tumour microenvironment is particularly pronounced.The tumour microenvironment of HNSCC is often accompanied by chronic inflammation, which promotes tumour development and metastasis. 42In addition, cells in the microenvironment can promote resistance of HNSCC cells to treatment.
Although our analysis is based on high-quality data and code, some limitations still need to be noted.Firstly, due to the systematic bias of bioinformatics itself, the results of its analysis cannot fully reflect the true situation within the organism.Secondly, potential biases between samples still affect the stability of the conclusion, such as lineage bias.Finally, although we identified EZR as the research gene, further experimental validation is still needed.

ACK N O WLE D G E M ENTS
None.

FU N D I N G I N FO R M ATI O N
None.

CO N FLI C T O F I NTE R E S T S TATE M E NT
None.
heatmap of Figure 1A was generated to visualize the overall distribution of the level of NK cells quantified by CIBERSORT, EPIC, MCPCOUNTER, QUANTISEQ, TIMER and XCELL algorithms in each patient.According to our data, we found that the level of NK cell abundance in EPIC, QUANTISEQ and XCELL algorithms F I G U R E 1 Identification of genes encoding markers for NK cells.(A) Quantitative analysis of the distribution of NK cells in each patient according to the CIBERSORT, EPIC, MCPCOUNTER, QUANTISEQ, TIMER, and XCELL algorithms.(B) Data distribution of the level of NK cell abundance.(C) Differential expression analysis of NK cell abundance in EPIC, XCELL and QUANTISEQ algorithms.(D,E) Wayne diagram of the differential up-regulated (D) and down-regulated (E) genes NK cell genes among three algorithms.F I G U R E 2 Analysis of the relationship between these genes and hub NK cell marker gene.(A) Correlation analyses between 8 genes associated with NK cell marker.(B) Univariate analysis of these genes and EZR was identified as the hub gene.(C) Differential expression analysis of these 8 NK cell marker genes.(D,E) Both GSE103322 and GSE139324 cohorts were used to analyse the expression distribution of these genes in multiple cell subgroups.

F I G U R E 3
Clinical relevance analysis.(A) The levels of EZR expression among all tumours based on TCGA dataset.(B-D) KM analyses of HNSC patients based on overall survival time, disease specific survival time and progress free interval time.(E) A comparison of two groups' clinicopathological features including age, gender, grade and stage.

F I G U R E 4
GSVA and GO enrichment analysis between two groups.(A) GSVA showing enriched pathways in two groups.(B,D) GO analysis showing BP and MF terms were enriched in high expression group.(C,E) GO analysis also revealed different BP and MF terms enriched in low expression group.

F I G U R E 5
Analysis of differential genomic differences between two groups.(A) The levels of TMB value among all tumours based on TCGA dataset.(B) Differential expression analysis of TMB level between two groups in HNSC patients.(C) Identification of differential mutant genes between two group.(D) The relationship of co-occurrence and exclusive between these differential mutant genes.(E) Data distribution of the level of mRNAsi score in each HNSC patient.(F) Differential expression analysis of mRNAsi score between two groups.
Tobacco, alcohol intake and HPV infection are the main risk factors.Despite great advances in surgery, radiotherapy and chemotherapy techniques over the past few decades, the treatment of HNSCC still faces many challenges, F I G U R E 7 Analysis of immune cell infiltration landscape in HNSC patients.(A-C) Differential expression analysis of ESTIMATE, Immune, and Stroma scores between two groups.(D,E) The boxplots for the comparison of immune cells and immune functions between two groups.(F) The boxplots for the comparison of the immune checkpoints genes between two groups.

F I G U R E 8 | 11 of 12 SHAO
Prediction for immunotherapy response and identification of appropriate agents.(A) The distribution of TIDE scores in HNSC patients.(B) The comparison of TIDE, Dysfunction and Exclusion scores between two groups.(C) The proportion of responder and nonresponder to immunotherapy in two groups.(D) Submap analysis for two groups in HNSC also suggest that patients in high expression group would be respond to immunotherapy.(E) Box plots of estimated drug sensitivity for several chemotherapeutic agents in two groups.et al.