Identification of several senescence‐associated genes signature in head and neck squamous cell carcinoma

Abstract Background As one of the core aging processes, cellular senescence is associated with tumorigenesis, growth, and immune modulation in cancers. Nevertheless, the prognosis of senescence‐associated genes (SAGs) signature in head and neck squamous cell carcinoma (HNSCC) remains to be further evaluated. Methods The transcriptome and corresponding clinical datasets of SAGs in patients with HNSCC were downloaded from public databases. A new prognostic SAGs signature was established with least absolute shrinkage and selection operator discussion. Patients with HNSCC were fallen into two risk groups based on each sample's risk mark and the cutoff point. The survival analysis was extended to determine the predictive accuracy of the SAGs signature. Furthermore, the evaluation of SAGs signature was made according to clinicopathological characteristics, survival state, the infiltration of inflammatory cells, and efficacy of immunotherapy. Results 41 SAGs were recognized and adopted to establish the forecast signature. The survival analysis indicated that patients with HNSCC in the high‐senescent score group had significantly reduced overall survival compared with those in the low‐senescent score group. It was certified that the risk score of SAGs signature was a separate predicting agent for HNSCC applying Cox regression analysis. According to functional analysis, some immune‐associated pathways were increased in the low‐senescent score group significantly. High‐senescent score group was correlated with poor clinicopathological characteristics, given less the infiltration of inflammatory cells state and worse immunotherapeutic effect. Conclusion A new SAG signature predicting result and response to immunotherapy of HNSCC was identified. Cellular senescence may be a hidden target for HNSCC.


| INTRODUC TI ON
Head and neck squamous cell carcinoma (HNSCC) ranks the sixth in terms of frequent cancer globally; there are about 0.6 million new cases diagnosed each year. 1 The primary treatment of HNSCC rely on the clinical phase including operation, radiotherapy, chemotherapy and immunotherapy. Treatments for a persistently high mortality rate in advanced HNSCC patients are limited; the 5-year total survival rate is only 63%. 2 For improving the prognosis results of HNSCC patients, steady prognostic signatures shall be identified.
Cellular senescence, one of the core aging processes, has the characteristics such as an eternal cell cycle arrest. 3 It is considered to be a response to injury inside and outside the cell. Nevertheless, the relationship of cell senescence and tumor, particularly in HNSCC, is complicated and poorly interpreted as yet. It has been deemed that cellular senescence is a double-edged sword during tumorigenesis and growth. In one aspect, aging cells go into eternal cell cycle arrest and be cleared by the human immune system, which inhibits cancer growth in prophase. 4 On the other hand, when immune system cannot effectively remove aging cells, this accumulation of senescence cells facilitates the senescence-related secretory phenotype (SASP), 5 which changes immune microenvironment and remodels tissue by releasing growth elements, cytokines, enzymes, and extracellular matrix (ECM) components, leading to tumorigenesis and development. 3 In the current studies, both an interior persistent DDR and extrinsic stress promote the progress of tumorigenesis and tumor cell evasion from immunosurveillance, bringing worse prognosis in HNSCC patients ( Figure 1A) 6,7 .
For assessing the relationship between cellular senescence and prognosis in HNSCC in an integrated way, we constructed a new risk model on basis of differentially expressed and prognostic senescence-associated genes (SAGs) and dug their hidden value as sibylline biomarkers for prognosis and immunotherapy response from databases. Next, the SAGs scoring model can forecast total survival of HNSCC patients and guide clinical treatment by calculating the total risk mark. In the end, hidden measures for treating HNSCC were identified by analyzing the association between cellular senescence, prognosis, the immune microenvironment, and the response to immunotherapy.

| Data sources and processing
CellAge (Avelar, et al., 2020) (https://genom ics.senes cence.info/ cells/) provided the list of genes with manually gathered information of human genes related to cellular senescence. 8

| Differential gene expression analysis
The differentially expressed genes (DEGs) between tumor samples (n = 499) and normal tissues (n = 44) was performed with R package "limma". Pearson associations >0.30 and p < 0.05 were adopted to set up the expression association network of SAGs in HNSCC.
F I G U R E 1 (A) Schematic diagram of the cellular senescence process. (B) The flow chart of this study. Abbreviations: pDDR, Persistent DNA damage response; SASP, senescence-associated secretory phenotype; HNSCC, Head and neck squamous cell carcinoma; TCGA, the cancer genome atlas; GEO, The Gene Expression Omnibus; SAGs, Senescence-Related Genes; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; TME, tumor microenvironment; GO, Gene Ontology; GSEA, Gene set enrichment analysis

| Establishment of a prognostic SAGs signature
Differentially expressed SAGs were screened by performing univariate Cox analysis applying R package "survival" based on the thresholds of p < 0.01. We identified and developed the optimal prognostic SAGs signature applying the "glmnet" R package by adopting the least absolute shrinkage and selection operator (LASSO) Cox regression discussion. The calculation of risk marks of HNSCC samples was made based on related regression parameters and expression level of each gene: Risk mark = ∑ n i=0 coefficient x different expression levels of SAGs. Furthermore, 499 HNSCC patients were fallen into high-and low-risk groups on basis of the average value of risk marks. The receiver operating characteristic (ROC) curve of 1-, 3-, and 5-years was constructed to assess the effectiveness of the SAGs signature model. The comparison of total survival diversities between two risk groups was made by performing Kaplan-Meier survival curves. If the SAGs signature was assessed as a hidden independent prognostic predictor of HNSCC patients applying univariate and multivariate Cox regression explorations.
A nomogram was established for survival prediction by including age, N stage, T stage, gender, grade, and risk score. The precision of the nomogram was assessed with the calibration curves. Chi-square tests and Wilcoxon signed-rank tests were adopted to assess clinicopathological elements.

| Functional enrichment and pathway analyses
On basis of DEGs, Gene Ontology (GO) analyses was made with a modified p < 0.05. Gene set enrichment analyses (GSEA, version 4.1.0) between two risk groups was utilized to find a possible molecular basis.

| Immune cell infiltration and immunotherapy
We calculated the immune mark, stromal mark and estimate mark of each HNSCC sample applying the ESTIMATE algorithm. 9 The cibersort method was adopted to find the fraction of immune-related cell types. 10 The different infiltrating levels of immune cell between two risk groups were calculated. The single-sample gene set enrichment analysis (ssGSEA) was adopted to obtain the immune-related pathway. The high-senescent score and low-senescent score groups were compared in immune checkpoint inhibitor (ICI)-associated genes. Immunotherapy response was evaluated to explore the value of ICI treatments for HNSCC by The Cancer Immunome Atlas (TCIA, https://tcia.at/) website 10 .

| Chemotherapeutic drug sensitivity analysis
In order to study the sensitivity of high-and low-senescent score groups to chemotherapy drugs, the "pRRophetic" package was used to evaluate HNSCC patients' half-maximum inhibitory concentration (IC50) of four conventional chemotherapy drugs (paclitaxel, gemcitabine, docetaxel, and cisplatin).

| Statistical analysis
All statistical analysis and figures were performed applying R language (version 4.1.0). The Chi-square test and Wilcoxon signedrank tests were employed for classified and constant variables data.
Pearson correlation analysis was adopted for association assessment. p < 0.05 meant statistical significance.

| Construction of a prognostic SAG signature
As shown in Figure 2A, 41 prognostic SAGs were recognized with univariate Cox regression analysis ( Table 1). The heatmap was set up in Figure 2B (red means high expression and green low expression). Figure 2C shows the expression association network of the SAGs in HNSCC (red displays positive association and blue negative association). Based on best value of λ and coefficient in LASSO cox regression analysis ( Figure 3A,B), the SAGs signature was constructed by choosing 23 prognostic SAGs (Table 2).

| Verification of the prognostic SAGs signature for HNSCC
We divided 499 patients with HNSCC into the high-senescent score group (n = 249) and the low-senescent score group (n = 250) applying the average risk mark as the cutoff value ( Figure 4A).
According to Figure  We constructed a nomogram to predict survival for HNSCC patients based on tumor grade, gender, age, tumor stage, and senescent risk score. As shown in Figure 4G, by calculating the total score of each feature of some patient marked in red, we can know that his 1-year, 3-year and 5-year survival rates are 57.9%, 21.7% and 13.3%, respectively. According to the calibration cures, the nomogram showed excellent prediction capacity by comparing with the ideal model ( Figure 4H). Then, the association between the clinicopathological elements and risk mark was analyzed in the heatmap showed that the grade (p < 0.001), clinical phase (p < 0.01), T stage (p < 0.001), and lymph node metastasis (p < 0.05) were greatly related to the risk score ( Figure 6).

| Functional enrichment analysis
According to Figure 7A, the results of GO and GSEA analysis indicated that the DEGs were mainly participated in the cellular or humoral immune response, complement activation, and immunoglobulin production. Uniformly, some immune-associated paths were greatly increased in the low-risk group, such as T-cell receptor signaling path, autoimmune thyroid disease, intestinal IgA generation and Fc epsilon Ri signaling path ( Figure 7B).

| The association between the SAGs signature and immune state
The tumor immune microenvironment of HNSCC was further evaluated. HNSCC patients in the low-senescent score group showed higher immune marks and estimate marks than the high-senescent  Figure 8D).

| Immunotherapy forecast of SAG signature for HNSCC
For researching the effect of immunotherapy, the relationship between risk marks and Immune checkpoint inhibitors (ICI) associated genes were analyzed. The expression level of ICI associated genes in the low-risk group were mostly greatly higher than the high-risk group ( Figure 9A, p < 0.001), indicating that the low-risk HNSCC patients better responded to ICI. This guess was confirmed with TCIA and revealed that the low-senescent score group better responded to PD-1 inhibiting agent alone (p = 0.0078), CTLA4 inhibiting agent alone (p = 0.02) or an integration of PD-1 and CTLA4 inhibiting agent (p = 0.0019) ( Figure 9B).

| Chemosensitivity analysis
By analyzing the IC50 values of four common chemotherapeutic drugs for HNSCC patients, we found that the IC50 value of gemcitabine ( Figure 10B, p < 0.001) was significantly reduced in the high-risk group, suggesting that patients at high risk of aging had a better response to the chemotherapy drug; However, IC50 values of paclitaxel ( Figure 10A), docetaxel ( Figure 10C), and cisplatin ( Figure 10D) did not differ significantly between high-and low-senescent score groups.

| DISCUSS ION
As the sixth most frequent cancer globally, head and neck squamous cell carcinoma (HNSCC) is driven by abnormal changes of genetic,  Cytokine Signaling in immune system according to gene database. [25][26][27][28][29][30] In brief, our findings suggest that these 23 prognostic SAGs signature significantly influence the survival time of cancer patients, which might offer hidden therapeutic objects for HNSCC. The To reveal the biological roles of SAGs signature, GO improvement analysis and GSEA analysis were made, showing that immuneassociated biological processes and pathways were increased in the low-senescent score group, which recruited more functional B F I G U R E 9 Immunotherapy prediction of SAGs-based prognostic model for HNSCC. (A) Immune checkpoint inhibitors related genes expression between highand low-risk patients. (B) Differences in immunophenoscores between patients in high-and low-risk groups received anti-PD1 alone, anti-CTLA4 alone, and combination therapy with anti-CTLA4 and anti-PD1 immunosuppressive cytokines were upregulated in patients with malignant melanoma at low risk of aging, which is consistent with our study 35,36 . In contrast, the high-senescent score group showed higher quantity of M2 macrophage cells, which was positively related to the risk mark. The research found that M2 macrophages are immune suppressive cells associated with angiogenesis and tissue remodeling, which enhance the immune suppression ability of tumors to promote cancer. 37 We summarize that cell senescence have systemic effects on tissue remodeling and immune microenvironment, which may mediate the pathogenesis of HNSCC.

TA B L E 2 Clinical characteristics of the HNSCC patients in this study
Current studies found that immunotherapy like pembrolizumab and nivolumab have been useful in HNSCC, inducing the human immune system to eliminate cancer. 11 In our study, we found that ICI associated genes (PD-L1, CTLA4, PD1) were greatly improved in the low-senescent score group, which had a better response to PD-1 inhibiting agent or CTLA4 inhibiting agent or the integration of PD-1 and CTLA4 inhibiting agent with TCIA. 38 We conclude that our SAGs signature might make contributions to developing immunotherapy measures for HNSCC patients.
Chemotherapy is one of the important treatment methods for advanced head and neck squamous cell carcinoma, and its clinical efficacy is worthy of affirmation. Some advanced patients will have significantly reduced tumor size after chemotherapy, and will get the opportunity for surgery again. Gadhikar et al. reported that cisplatin inhibits cell growth by inducing senescence of head and neck cancer cells. 39 Therefore, aging genes may be associated with chemotherapy resistance of squamous cell carcinoma of the head and neck. In order to further obtain the predictive value of aging genes in chemotherapy response of HNSCC patients, we calculated the sensitivity of four commonly used chemotherapy drugs to HNSCC patients according to the "pRRophetic" algorithm.
Interestingly, our study found that gemcitabine was more sensitive in patients with high-senescent risk, which provides a priority reference value for rational clinical selection of chemotherapy agents for HNSCC patients.
There are several limitations in this research. The conclusions would be more useful in case of experimental verification. On basis of the research outcomes, new research will perform HNSCC animal models to dig the association of 23 senescence-associated genes expression and immune microenvironment.

| CON CLUS IONS
A new 23 SAGs signature was established for predicting outcomes of HNSCC. SAGs potentially influences tissue remodeling and immune microenvironment for the occurrence and development of HNSCC.
Our outcomes offer insight to forecast results and recognize therapeutic targets for HNSCC patients.

CO N FLI C T O F I NTE R E S T
The author has declared that there is no competing interest associated with this study.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data used for our analysis in this study are openly available at (https://portal.gdc.cancer.gov/ and http://www.ncbi.nlm.nih.gov/ geo/).