Research Paper|Volume 16, Issue 2|pp 985—1001

Prognostic significance of senescence-related tumor microenvironment genes in head and neck squamous cell carcinoma

Young Chan Lee^1,^2,³, Yonghyun Nam³, Minjeong Kim², Su Il Kim¹, Jung-Woo Lee⁴, Young-Gyu Eun⁵, Dokyoon Kim^3,⁶

¹Department of Otolaryngology-Head and Neck Surgery, Kyung Hee University School of Medicine, Kyung Hee University Hospital at Gangdong, Seoul, Republic of Korea
²Department of Medicine (AgeTech-Service Convergence Major) College of Medicine, Kyung Hee University, Seoul, Republic of Korea
³Department of Biostatistics, Epidemiology and Informatics, The Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA 19104, USA
⁴Department of Oral and Maxillofacial Surgery, School of Dentistry, Kyung Hee University, Seoul, Republic of Korea
⁵Department of Otolaryngology-Head and Neck Surgery, Kyung Hee University School of Medicine, Kyung Hee University Medical Center, Seoul, Republic of Korea
⁶Institute for Biomedical Informatics, University of Pennsylvania, Philadelphia, PA 19104, USA

* Equal contribution

Received: August 24, 2023Accepted: November 6, 2023Published: December 21, 2023

https://doi.org/10.18632/aging.205346

Copyright: © 2024 Lee et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

The impact of the senescence related microenvironment on cancer prognosis and therapeutic response remains poorly understood. In this study, we investigated the prognostic significance of senescence related tumor microenvironment genes (PSTGs) and their potential implications for immunotherapy response. Using the Cancer Genome Atlas- head and neck squamous cell carcinoma (HNSC) data, we identified two subtypes based on the expression of PSTGs, acquired from tumor-associated senescence genes, tumor microenvironment (TME)-related genes, and immune-related genes, using consensus clustering. Using the LASSO, we constructed a risk model consisting of senescence related TME core genes (STCGs). The two subtypes exhibited significant differences in prognosis, genetic alterations, methylation patterns, and enriched pathways, and immune infiltration. Our risk model stratified patients into high-risk and low-risk groups and validated in independent cohorts. The high-risk group showed poorer prognosis and immune inactivation, suggesting reduced responsiveness to immunotherapy. Additionally, we observed a significant enrichment of STCGs in stromal cells using single-cell RNA transcriptome data. Our findings highlight the importance of the senescence related TME in HNSC prognosis and response to immunotherapy. This study contributes to a deeper understanding of the complex interplay between senescence and the TME, with potential implications for precision medicine and personalized treatment approaches in HNSC.

Introduction

Head and neck cancer (HNC) is the sixth most common cancer, with more than 68,000 patients reported to be newly diagnosed in the United States in 2021 [1]. The 5-year survival rate for HNC ranges from 30% to 70%, contingent upon the tumor's stage and site [2]. This heterogeneity arises from malignant cells harboring diverse genotypes, phenotypes, and interactions within the tumor microenvironment (TME) within each individual tumor, significantly contributing to tumorigenesis and malignant progression, presenting a major obstacle for cancer therapeutics [3]. Cellular senescence, characterized by a state of cell-cycle arrest, can be induced by various profound internal or external stresses, including oncogenic activation or DNA damage from chemotherapeutic agents [4]. The crucial role of senescent cells as a pivotal component of the tumor microenvironment (TME) has been highlighted in Hallmarks of Cancer [3]. Over time, numerous researchers have regarded cellular senescence as a mechanism against malignancy, leading to the conversion of cancer cells into senescent cells [5]. However, emerging evidence in recent years has revealed the dualistic roles of senescent cells, which can either impede or promote tumor development and malignant progression, depending on the specific conditions under which they are induced [6, 7]. In particular, senescent cells secrete several cytokines and growth factors known as the senescence-associated secretory phenotype (SASP). It has been reported that the SASP fosters a close relationship between senescent cells and the TME, inducing immunosuppression and inflammation to promote tumor growth and potentially influencing the response to immunotherapy [8–10]. Moreover, the targeting of senescence process has emerged as a promising strategy in cancer therapeutics [11]. Nevertheless, our current understanding of the interaction between senescence and TME remains limited. Specifically, it remains unclear whether the senescence related TME characteristics observed in patients with head and neck squamous cell carcinoma (HNSC) can serve as predictive biomarkers for clinical prognosis and therapeutic response.

In this study, we identified prognostic senescence-related TME genes (PSTGs) through a comprehensive analysis involving gene-gene network, differential expression analysis, and Cox regression analysis. Furthermore, we conducted clustering analysis using the expression of these genes. Subsequently, we developed a novel risk score model based on senescence related TME core genes (STCGs) to predict patient prognosis and response to immunotherapy. Importantly, the predictive capability of this model was validated in independent cohorts. Additionally, we explored the enrichment patterns of these core genes at the HNSC single-cell level.

Results

The identification of senescence related TME subtypes and characterization

Based on the transcriptome expression levels of PSTGs, the Consensus Clustering Method, an unsupervised ensemble clustering algorithm, was employed to cluster the samples in TCGA-HNSC cohort. The consensus matrix CDF curve revealed a distinct flattened segment at K=2, indicating optimal subgrouping (Figure 1C, 1D). Additionally, selecting K=2 for consensus clustering analysis minimized interference between subtypes. A total of 500 HNSC patients in the TCGA database were divided into two subtypes, named subtype 1 and subtype 2.

To better understand the characteristics of the two subtypes, the expression levels of PSTGs and various clinicopathologic features were compared between the two subtypes, as illustrated in Figure 2A. The two senescence-related TME subtypes differed in their clinical and molecular characteristics. Subtype 2 demonstrated a higher prevalence of HPV-positive and oropharyngeal cancer cases, whereas subtype 1 was characterized by a higher proportion of advanced T stage, N stage and overall stages.

Figure 2.

Distinct characteristics of the senescence related TME subtypes. (A) Heatmap showing the expression of 91 PSTSs of the two subtypes. Red represents high expression and blue represents low expression. Subtype, subsite, stage, T, N and HPV status were used as sample annotations. (B) differentially enriched pathway (subtype 2 vs subtype 1). (C) Forest plot showing somatic mutation in the two subtypes. (D) Barplot showing somatic mutation in the two subtypes. (E) Top 10 most significantly differentially methylated CpGs under comparing subtype 1 and 2. (F) Differences in SASP and GO_senescence ssGSEA score between the two subtypes in the TCGA–HNSC cohort.

We conducted an analysis to investigate the presence of differential enrichment in gene ontology (GO) and pathways between the two subtypes based on the DEGs (Figure 2B). The enriched GO terms encompassed biological processes such as immune response and defense response, cellular components including intrinsic and integral to the plasma membrane, and molecular functions such as cytokine binding and MHC class II receptor activity. Additionally, the enriched pathways included crosstalk between dendritic cells and natural killer cells, as well as altered T cell and B cell signaling in rheumatoid arthritis.

We investigated the difference in genomic alteration between the subtypes (Figure 2C, 2D). In subtype 1, the most frequently observed mutation was in TP53, found in 84% of the samples, followed by TTN, CSMD3, and FAT1. On the other hand, subtype 2 exhibited common mutations in TP53, TTN, FAT1, CDKN2A, and PIK3CA. While CDKN2A and FAT1 mutations were primarily nonsense mutations, the majority of mutations were missense mutations. The TMB was found to be higher in subtype 1 compared to subtype 2 (Subtype1; median 2/MB, Subtype2; median 1.92/MB). TP53 and NSD1 mutations were more frequent in subtype 1, whereas CASP8 mutations were more prevalent in subtype 2.

We identified the top 10 most significantly differentially methylated CpGs (p-value ≤0.05; FC >1; Δβ-value >0) when comparing subtype 1 and 2 (Figure 2E). Nine CpG sites (cg23967461, cg20926024, cg24749947, cg06121450, cg14873488, cg25467652, cg15995296, cg01943873, and cg06422467) were hypomethylated in subtype 1 compared to subtype 2, and one CpG site (cg15375424) was hypermethylated in subtype 1. These CpG sites of deregulation of DNA methylation were in ACVR1, SIL1, TMCO1, MIA2, AGTRAP, CORO1B, HCG20 and IRF1. Metascape analysis showed that the genes suppressed by DNA methylation in subtype 1 were enriched in the ZNF528 target genes among the Transcription Factor Target ontology. SASP and GO_senescence ssGSEA scores were calculated for each sample in TCGA-HNSC to represent the senescence state using the ssGSEA method from the GSVA package, according to the SenMayo gene set and GOBP_MULTICELLULAR_ORGANISM_AGING [12]. It was observed that the SASP ssGSEA score was significantly higher in subtype 2 than in subtype 1 (p < 0.05), and the GO_senescence ssGSEA score of the patients in subtype 1 was significantly higher than that of the patients in subtype 1 (p < 0.05) (Figure 2F).

The Kaplan-Meier survival curves revealed significant differences in overall survival (OS) (p = 0.012) between the two subtypes, with subtype 2 demonstrating a more favorable prognosis compared to subtype 1. Furthermore, the significance of the prognostic difference between the clusters was enhanced in the HPV-positive cohort (p = 0.0047), while it diminished in the HPV-negative cohort (p = 0.25). The difference in disease-specific survival (DSS) between the two subtypes exhibited a consistent pattern comparable to that observed for overall survival (Figure 3A, 3B). Remarkably, within the subgroup subjected to radiation treatment, a substantial escalation in prognostic divergence between the two subtypes was observed (OS; p=0.0025, DSS; p=0.00055) (Figure 3C).

Figure 3.

Prognostic difference between the two subtypes. (A) KM curves indicating prognostic differences in overall survival between the two subtypes in the TCGA-HNSC cohort (Left; total, Middle; HPV positive cohort, Right; HPV negative). (B) KM curves indicating prognostic differences in disease specific survival between the two subtypes in the TCGA-HNSC cohort (Left; total, Middle; HPV positive cohort, Right; HPV negative). (C) KM curves indicating prognostic differences in overall survival and disease specific survival between the two subtypes in the TCGA-HNSC cohort subgroup who received radiotherapy.

We utilized immune cell signatures to evaluate the infiltration of immune cells between two subtypes, aiming to assess their distinct immune characteristics. Through the application of CIBERSORT, we discovered significant differences in 16 immune cell types including B cell, CD4 memory activated T cell, macrophages, and dendritic cells between the subtypes (Figure 4A).

Figure 4.

Difference of immune characteristics between the two subtypes. (A) Differences in immune cell scores between two subtypes in the TCGA–HNSC cohort. (B) Differences in ESTIMATE immune infiltration between two subtypes in the TCGA–HNSC cohort. (C) Differences in TIDE analysis score between two subtypes in the TCGA–HNSC cohort. (D) Differences in immunotherapy response prediction between two subtypes in the TCGA–HNSC cohort.

To evaluate the characteristics of senescence-related TME subtypes on immunity, we performed the ESTIMATE algorithm. When comparing the immune scores (the proportion of immune cells), stromal score (the proportion of stromal cells) and ESTIMATE score (the proportion of nontumor components) between the two subtypes, we observed that the stromal score, immune score, and ESTIMATE score were significantly higher in subtype 2 (Figure 4B).

TIDE analysis predicted distinct responses to immune checkpoint inhibitors between the two subtypes. Subtype 1 exhibited higher TIDE and Exclusion scores compared to subtype 2, indicating a greater likelihood of immunotherapy evasion. Furthermore, subtype 1 displayed elevated proportions of MDSC, TAM.M2, and CAF, suggesting increased T cell exclusion (Figure 4C). These contrasting immune characteristics between the subtypes led us to hypothesize that they would elicit different responses to immunotherapy. Consistently, the proportion of responders to immunotherapy was significantly higher in subtype 2 (Responder in subtype 1 (24.7%) vs subtype 2(34.1%), p<0.001) (Figure 4D).

Construction of risk scoring model based on senescence related TME status

STCGs were selected by performing the LASSO Cox regression algorithm to select penalty coefficient using 91 PSTGs. The minimum log(lambda) was set to the optimal value through 5-fold cross-validation (Figure 5A, 5B).

Figure 5.

Construction of risk scoring model based on senescence related TME status. (A) LASSO PLD. (B) LASSO c-index. (C) core genes, coefficient. (D) Heatmap showing the expression of 21 STCGs of the two risk groups. (E) KM curves indicating prognostic differences between the two risk groups (Left; HNSC-TCGA test cohort, Middle; GSE41613, Right; KHUMC cohort). (F) Receiver operating characteristic (ROC) curves predicting 1-year, 2-year, 3-year, 4-year OS. (G) Multivariable analysis with hazard ratio (HR) for OS represented in a Forest plot.

Finally, several risks (OLR1, VEGFC, ITGA5, P4HA1, TINAGL1, TBX3, FGF7, PPARG, EDA2R) and protective (CDKN2A, ADAM33, TNFRSF4, SOCS1, TNFRSF25, MYO1G, CD38, FCRLA, EGFL6, ICOS, COL8A2, LYZ) mRNAs were identified in patients with HNSC (Figure 5C).

The entire TCGA-HNSC set was split into a training set and a test set at a ratio of 7:3. We calculated risk scores based on the expression level of each gene and the coefficient in training set. The median risk score was used as the cutoff value to classify each individual in the TCGA-HNSC cohort as high-risk or low-risk. The heatmap in Figure shows the difference in expression of 21 genes between the two risk groups (Figure 5D).

Prognostic performance of this STCGs based risk model was tested in test set, independent cohorts. Based on the risk model, the test set, GEO cohort (GSE41613), and KHUMC cohort were stratified into high-risk and low-risk groups. Survival analysis demonstrated that patients in the high-risk group exhibited significantly lower overall survival compared to those in the low-risk group in the test set (high-risk (n=49) vs low-risk (n=48), log-rank test, p-value=0.026) and the GEO cohort (high-risk (n=33) vs low-risk (n=64), log-rank test, p-value=0.049). Furthermore, the low-risk group demonstrated a trend towards higher overall survival in the KHUMC cohort (independent cohort, Asian ancestry) (high-risk (n=29) vs low-risk (n=46), log-rank test, p-value=0.079), although it did not reach statistical significance (Figure 5E).

We also calculated the predictive ability for 1-year, 2-year, 3-year, 4-year OS using ROC analysis, and the area under the ROC curve was 0.696, 0.685, 0.703, and 0.686, respectively (Figure 5F). In multivariate Cox regression analysis, the risk score was the most important independent prognostic factor to predict OS (HR: 0.36, CI: 0.24- 0.54, p <0.001) (Figure 5G). Therefore, it was confirmed that senescence related TME risk score can be used as a prognostic biomarker for HNSC.

Immunotherapy response prediction of senescence related TME risk model and STCGs expression in single cell level

In previous studies, the IPS has demonstrated its predictive value for immune checkpoint inhibitor (ICI) response in melanoma patients, owing to its high immunogenic potential [13]. In our investigation, we observed distinct IPS patterns between the two risk groups among HNSC patients. Specifically, the low-risk group exhibited significantly increased IPS scores in CTLA4 negative/PD-1 negative, CTLA4 negative/PD-1 positive, CTLA4 positive/PD-1 negative, and CTLA4 positive/PD-1 positive. Furthermore, when utilizing the TIDE score to predict ICI response, we observed a significantly higher proportion of responders in the low-risk group (Figure 6A, 6B). These findings suggest that HNSC patients with a low risk of senescence EMT represent promising candidates for ICI therapy. Alluvial diagram showed the mutuality of the molecular subtypes [14], subtype and risk score group in TCGA-HNSC cohort. Mainly, basal and classical subtypes were associated with Subtype 1 and high-risk score group. The atypical subtype associated with HPV positive has a significant correlation with subtype 2 and low risk score group (Figure 6C).

Figure 6.

Immunotherapy response prediction of senescence related TME risk model and STCGs expression in single cell level. (A) Differences in immunophenoscore between two risk group in the TCGA–HNSC cohort. (B) Differences in immunotherapy response prediction between two risk group in the TCGA–HNSC cohort. (C) Alluvial plot showing the changes of molecular subtype, senescence related TME subtype and risk group. (D) Single cell profiling of senescence related TME core genes.

The TME, consisting of immune cells, inflammatory cells, and stromal cells, exerts a crucial role in the initiation, progression, metastasis, recurrence, and acquisition of drug resistance in tumor. Therefore, to further evaluate the senescence related TME characteristics at the single cell level, we investigated the single-cell transcriptome of GEO dataset from primary HNSC tissue. In UMAP and violin plots, STCGs were enriched in fibroblast, mono/macrophage, and T cells rather than cancer cells, suggesting that these cell types, and not cancer cells, contribute to the senescent features of HNSC. Interestingly, expression of risk genes was increased in fibroblast and endothelial cells, and protective genes were mainly expressed in T cells (Figure 6D). In short, the above results suggest a significant association between STCGs and HNSC, with a particularly notable relationship observed in stromal and immune cell populations.

Discussion

The long-lasting impact of senescent cells on tissue homeostasis has gained prominence, particularly with the identification of the SASP [15]. SASP serves a dual role by not only reinforcing cellular senescence through autocrine signaling but also mediating paracrine effects [16]. Through paracrine signaling, SASP factors have the capability to remodel tissues, impacting the proliferation and migration abilities of neighboring cells such as stromal cells, immune cells, and cancer cells [4, 17]. In addition, SASP factors possess the potential to stimulate angiogenesis and augment the immunosuppressive microenvironment [5]. Coppé et al. found that SASP factors selectively act on immune cells and stromal cells present in the TME, triggering paracrine senescence [18]. This process promotes cancer cell epithelial-mesenchymal transition (EMT) and enhances invasiveness. However, the knowledge regarding the link between cellular senescence in the stroma and TME has remained significantly limited thus far. Furthermore, the impact of cellular senescence associated TME on the efficacy of cancer treatments, including immunotherapy, and its potential as a prognostic indicator remains elusive. To the best of our knowledge, this is the first study to offer a comprehensive evaluation of the senescence related TME status by integrating senescence related TME genes through a gene-gene network and clustering. Furthermore, we have introduced a novel risk model that utilizes a selected gene set to predict prognosis and confirmed the expression of STCGs in immune cells at single-cell levels.

Initially, we identified a set of genes from the intersection in the lists of TAS genes, TME related genes, and immune-related genes. Leveraging these genes as seed nodes, we successfully derived interconnected gene lists within the same network through gene-gene interactions. This proposed model holds the potential to enhance the accuracy and efficacy of selected genes associated with senescence in the TME.

Through the application of consensus clustering using PSTGs, we identified two distinct subtypes. These subtypes exhibited significant differences in mutation profiles, methylation patterns, immune profiles reflecting the TME status, prognosis, and immunotherapy response. Subgroup 1 is more aggressive compared with tumor in subgroup 2. Interestingly, in the analysis of HNSC patients who received radiotherapy, the prognosis of patients in subgroup 1 was worse, and the difference in prognosis between subtypes increased more significantly. This observation suggests that the expression of radiation-induced senescence related genes may have implications for prognosis in the context of this study. The higher mutation rate of the TP53 gene may be implicated in subtype 1 concerning senescence related TME. The loss of p53 function promotes chromosomal instability, leading cells to undergo either senescence or apoptosis through direct and indirect mechanisms [19]. Numerous studies have indicated the significance of DNA methylation patterns in senescence as a pivotal factor influencing tumor behavior [20, 21]. In our results, subtype 1 patients had significantly hypomethylation in several CpG sites than subtype 2 patients and methylation-silenced genes in subtype 1 were enriched in transcription factor target, especially ZNF528 target genes. The data suggest that epigenetic silencing of ZNF528 could be an important factor in the determination of senescence related TME subtype in HNSC.

Furthermore, we observed the significant difference of immune cell profile and TME status between the two subtypes. The results showed that subtype 2 showed significant increases in the infiltration of immune cells such as the activated CD4+ T cells, B cell and macrophages. In addition, there was a significant difference in stromal score, immune score and ESTIAMTE score between two subgroups. The results of this study provide evidence supporting the association between senescence related TME subtypes and distinct TME features. In this context, analysis of functional differences in DEGs indicates that subtype 2 is closely related primarily to immune responses, defense responses, regulation of lymphocyte activation, and crosstalk between dendritic cells and natural killer cells. Interestingly, in the comparison of ssGSEA between the two subgroups, the SASP and senescence gene sets exhibited contradictory findings. It is well-established that senescent cells typically release SASP factors, however, we observed a decrease in SASP expression in subtype 1, characterized by high expression of senescence gene sets. This discrepancy can be attributed to the fact that SASP-associated genes may predominantly reflect immune cell activity within the tumor microenvironment (TME). Conversely, the upregulation of the senescence gene set likely arises from an augmented occurrence of cellular senescence or impaired immune-mediated clearance of senescent cells [22].

Furthermore, we established STCGs including twenty-one genes (OLR1, VEGFC, ITGA5, P4HA1, TINAGL1, TBX3, FGF7, PPARG, EDA2R, CDKN2A, ADAM33, TNFRSF4, SOCS1, TNFRSF25, MYO1G, CD38, FCRLA, EGFL6, ICOS, COL8A2, LYZ) were selected by performing the LASSO Cox regression algorithm for predicting the prognosis and therapy response of HNSC patients. Previous studies have provided some level of understanding regarding the biological functions of the genes encompassed in STCGs. Peroxisome proliferator-activated receptors (PPARs) are ligand-activated transcription factors that belong to the nuclear hormone receptor superfamily [23]. PPARG is a key regulator of lipid metabolism in diverse immune cells, thus exerting a significant influence on immune regulation [24, 25]. Notably, recent investigations have revealed an unfavorable prognostic association between overexpression of PPARG and certain cancer types (e.g., prostate cancer, esophageal cancer), particularly those linked to obesity [26, 27]. Given the highest coefficient of PPARG as a risk core gene in our findings, further exploration into the interplay between lipid metabolism and the senescence process in cancer is warranted [28]. Furthermore, consistent with our findings, previous studies have reported that OLR1 and FGF7 are associated with adverse cancer prognosis due to alterations in immune response and TME status [29–31].

In this study, we have validated the predictive capability of the risk score obtained from STCGs across two independent cohorts. This validation emphasizes the robustness and reliability of the model based on STCGs. Notably, a significant disparity in HPV status was observed between the two classified subtypes, leading to the anticipation that our risk model would solely exhibit efficacy in HPV-positive cancer. Surprisingly, in the validation dataset (GSE41613), our developed model demonstrated prognostic capability even in HPV-negative oral cavity cancer patients. Furthermore, in the dataset comprising Asian patients (KHUMC cohort) with different ancestral backgrounds than the TCGA cohort, it was observed that the risk score differentiated the prognosis, although statistical significance was not attained.

To evaluate the expression of STCGs in TME cells, we identified TME cell populations using scRNA-seq data. We found a particularly high expression of risk STCGs in fibroblast and endothelial cells. The bioactive Vascular endothelial growth factor (VEGF) secreted by senescent fibroblasts has a significant impact on tumor angiogenesis and the progression of cancer [32]. Additionally, studies have demonstrated that senescent endothelial cells contribute to the enhanced aggressiveness of breast cancer cells [33]. Thus, our findings strongly suggest that cellular senescence in these specific cell types within TME exhibits pro-tumor properties. Conversely, a notable upregulation in the expression of protective STCGs was observed specifically in T cells, indicating the potential of these immune cells to exert tumor control mechanisms via senescence.

Our study has several limitations. Firstly, despite utilizing public data for the development and validation of risk models associated with senescence, it is imperative to validate these models using prospective data from multicenter studies to enhance their applicability. Secondly, the transcriptional profiles analyzed in this study were generated on RNA sequencing, and the gene expression of the identified STCGs were further confirmed in single cell RNA sequencing in HNSCC tissues. However, several protein products in the risk model were not validated in either HNSCC cell lines or HNSC tumor tissues. Our future research will contain some experiments such as immunohistochemical testing or western blot for expression of STCG in HNSCC tissue. Thirdly, while a risk model based on STCGs has been established, there is a need to improve its diagnostic performance by integrating relevant clinical parameters. Finally, although hypotheses regarding the functions and mechanisms of senescence-associated genes within the TME have been proposed, further comprehensive investigations are essential to unravel the specific mechanisms underlying their actions.

In conclusion, this study comprehensively investigated the prognostic and immunological features of senescence related TME genes in HNSC. By leveraging these senescence related TME genes, we successfully developed a risk model to predict HNSC prognosis and immunotherapy response, which was robustly validated using external transcriptome datasets. These findings provided evidence for the role of senescence in the TME and highlighted the potential of senescence-related biomarkers as promising therapeutic targets.

Materials and Methods

Data acquisition and processing

RNA-seq derived gene expression data (N=520) which transformed into transcripts per kilobase million (TPM), somatic mutation data (N=511), Human Methylation 450 data (N=528) and clinicopathologic data for HNSC were downloaded from the TCGA database using the R package ‘TCGAbiolinks’ (version 3.15) [34] Microarray gene expression data (GSE41613, N=167) from GEO database were obtained and processed to normalized matrix data by GEOquery R package and used as a validation dataset [35].

A total of 1,889 Tumor associated senescence (TAS) were selected based on well-established databases and published literature, including MSigDB (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp), SenMayo gene set and, The Human Ageing Genomic Resources (HAGR) [12, 36]. The immune-related genes (N= 4,723) were downloaded from ImmPort database [37]. TME-related genes were obtained from several previously published studies and selected by removing duplicates [38]. The gene lists were expanded by searching for senescence related TME genes through network analysis with a gene-gene network since there may be genes that have not yet been investigated among the collected gene lists in this study. A gene-gene network is a weighted graph representing genes as nodes and connections between genes as edges. The edges were measured by the magnitude of the correlation coefficients, indicating the similarity between two genes in the network. Network-based label propagations were performed by employing graph-based semi-supervised learning (SSL) to investigate the potentially associated genes related to prognostic senescence, propagating label information on query nodes to unlabeled nodes along with edges [39].

KHUMC cohort

Between January 2011 and January 2019, we enrolled 72 patients diagnosed with HNSC at the Kyung Hee University Medical Center (KHUMC) who received curative treatment. The cohort, referred to as the KHUMC cohort, was prospectively followed up for over 5 years following treatment, during which their clinical data, including age, sex, smoking history, and treatment type, were obtained. Furthermore, survival and recurrence information were retrospectively assessed for each patient.

Tissue samples were collected immediately after surgery or biopsy from patients diagnosed with HNSC. Total RNA was extracted from these samples using The TRIzol^® Reagent. To generate cDNA, the extracted total RNA was reverse transcribed using the Tetro cDNA Synthesis Kit (Bioline, London, UK) according to the manufacturer’s protocol. Subsequently, the isolated total RNA samples were sent to Applied Biosystems Macrogen Korea for sequencing, pre-processing, and transcriptome analysis. The concentration of total RNA was determined using Quant-IT RiboGreen^® (Invitrogen, #R11490). Only high-quality RNA preparations, with an RNA Integrity Number (RIN) greater than 7.0, were selected for RNA library construction. Each sample was independently used to prepare a library with 1 μg of total RNA, employing the Illumina TruSeq^® Stranded mRNA Sample Prep Kit (Illumina Inc., #RS-122-2101). The libraries were quantified using KAPA Library Quantification Kits for Illumina Sequencing platforms, following the qPCR Quantification Protocol Guide (Kapa Biosystems, Wilmington, MA, USA, #KK4854), and their quality was assessed using the TapeStation D1000 ScreenTape (Agilent Technologies, Palo Alto, CA, USA, #5067-5582). Indexed libraries were then submitted to an Illumina NovaSeq (Illumina Inc.), and Macrogen Inc. (Seoul, Korea) performed paired-end (2 × 100 bp) sequencing.

Clustering analysis for senescence related TME subtypes and characterization

To identify genes with differential expression between tumor and normal samples, we employed the Limma R package. We further selected differentially expressed genes that showed prognostic significance using the univariate Cox regression analysis.

To cluster HNSC samples, we utilized ConsensusClusterPlus to construct a consistency matrix [40]. The expression data of genes linked to senescence related TME served as the basis for determining the molecular subtypes of the samples. The “hc” algorithm (hierarchical clustering) and “1-Pearson correlation” were chosen as the metric distance for conducting 500 bootstraps with each bootstrap iteration involved 80% of the patients in the training set. We performed a grid search ranging from 2 to 10 clusters to determine the optimal number of clusters. The optimal cluster was selected based on the cumulative distribution function (CDF) and consensus matrix.

To explore the relationship between senescence-related TME subtypes and senescence-associated biological functions across different samples, we conducted a single-sample gene set enrichment analysis (ssGESA) using the R package GSVA [41]. Subsequently, ssGSEA scores were calculated for each sample on various functions for comparison between the subtypes. To identify significant differences in gene mutations between the senescence related TME subtypes in the TCGA-HNSC dataset, we used the maftools R package (version 2.6.05). Based on the subtype, we categorized the original mutation annotation format (MAF) into two distinct groups. Tumor mutation burden (TMB) scores were computed for each patient within the two subtypes using somatic mutation data. For the analysis of methylation, TCGA-HNSC HumanMethylation450 data was used and preprocessed using the Limma R package (v3.46.0) with an FDR-corrected P-value of 0.01 and absolute log fold change > 0.5 to identify differentially methylated CpGs between subtypes. Genome annotation of these differentially methylated probes was based on the Illumina protocol. Metascape (https://metascape.org) analysis was performed to identify the pathways associated with a gene set.

To compare the difference in survival between the two subtypes, Kaplan-Meier survival analysis was performed using the ‘survival’ and ‘survminer’ R packages. The Kaplan-Meier curves were utilized to generate p-values and hazard ratios (HR) with 95% confidence intervals (CI) through log-rank tests. For visual representation of the p-values, HR, and 95% CI for each variable, a forest plot was constructed using the ‘forestplot’ R package.

Estimation of the immune cell landscape, immunophenoscore and prediction of immunotherapy responsiveness

CIBERSORT utilized a gene expression signature matrix derived from purified immune cell populations to deconvolute the composition of complex tissue samples [43]. It leverages support vector regression and an artificial immune system-like algorithm to estimate the relative proportions of different immune cell types in HNSC sample. The accuracy of immune cell fractionation was considered significant when the CIBERSORT output achieved a p-value of less than 0.05. For each HNSC sample, the stromal and immune scores were estimated by applying the ESTIMATE algorithm to the normalized expression matrix [44]. Data on individual immunophenoscore (IPS) for HNSC patients were obtained by downloading from the Cancer Immunome Atlas (https://tcia.at/home).

To validate the predicted treatment responsiveness, we employed the TIDE algorithm, which utilizes gene expression profiles as a computational method for predicting the efficacy of immune checkpoint blockade [45]. We also calculated the M2 tumor-associated macrophages (TAMs), tumor-associated fibroblasts (CAFs), and myeloid-derived suppressor cells (MDSCs), the dysfunction score of tumor-infiltrating cytotoxic T lymphocytes (CTLs) (T cell dysfunction), and the exclusion score of CTLs by immunosuppressive factors (T cell exclusion) through TIDE analysis.

Single-cell RNA-seq data

To understand the expression of STCGs in different cell types, we applied Tumor Immune Single Cell Hub (TISCH, http://tisch.comp-genomics.org/). TISCH is an online database dedicated to the TME and comprises a comprehensive collection of 76 tumor datasets encompassing 27 types of cancer, including single-cell transcriptome profiles comprising nearly 2 million cells. In our study, we focused on examining the expression patterns of the STCGs in HNSC sample. To accomplish this, we utilized GSE103322, a HNSC single-cell RNA sequencing dataset, which is part of the extensive data available in the TISCH database [46]. GSE103322 contained data of 5,902 cells derived from 18 patients with oral cavity tumor. We explored the expression of the STCGs in HNSC at single-cell level and identified the distribution of expression of senescence-related TME core genes in GSE103322. The expression of STCGs was collapsed by mean value, and the gene expression level displayed using UMAP and violin plots was quantified by the normalized values.

Supplementary Materials

Supplementary Table 1

Author Contributions

YCL conceived, designed and supervised the study, developed the methodology, acquired data, conducted data analysis and interpretation, and wrote the manuscript. YN acquired, analyzed and interpreted data, and wrote the manuscript. MK, SIK and JL were involved in the development of the methodology, acquired data, and reviewed the manuscript. YE and DK contributed to the conception and design of the study as well as to the development of the methodology, provided facilities, performed data analysis and interpretation, and participated in writing the manuscript. All authors read and approved the final manuscript.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

The experiments were undertaken with the understanding and written consent of each subject. The study methodologies were approved by our institutional review board (IRB: 2018-05-046-12) and were conducted according to the principles of the Declaration of Helsinki.

Funding

This work was supported by the BK21 plus program “AgeTech-Service Convergence Major” through the National Research Foundation (NRF) funded by the Ministry of Education of Korea [5120200313836] and the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) [RS-2023-00241230]. This work was also funded by National Library of Medicine (NLM) R01 LM012535.

References

1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. CA Cancer J Clin. 2021; 71:7–33. https://doi.org/10.3322/caac.21654 [PubMed]
2. Nederland IK. Netherlands Comprehensive Cancer Organisation. Neth Cancer Regist Age Gend. 2014.
3. Hanahan D. Hallmarks of Cancer: New Dimensions. Cancer Discov. 2022; 12:31–46. https://doi.org/10.1158/2159-8290.CD-21-1059 [PubMed]
4. Gorgoulis V, Adams PD, Alimonti A, Bennett DC, Bischof O, Bishop C, Campisi J, Collado M, Evangelou K, Ferbeyre G, Gil J, Hara E, Krizhanovsky V, et al. Cellular Senescence: Defining a Path Forward. Cell. 2019; 179:813–27. https://doi.org/10.1016/j.cell.2019.10.005 [PubMed]
5. Lee S, Schmitt CA. The dynamic nature of senescence in cancer. Nat Cell Biol. 2019; 21:94–101. https://doi.org/10.1038/s41556-018-0249-2 [PubMed]
6. Wang B, Kohli J, Demaria M. Senescent Cells in Cancer Therapy: Friends or Foes? Trends Cancer. 2020; 6:838–57. https://doi.org/10.1016/j.trecan.2020.05.004 [PubMed]
7. Kowald A, Passos JF, Kirkwood TB. On the evolution of cellular senescence. Aging Cell. 2020; 19:e13270. https://doi.org/10.1111/acel.13270 [PubMed]
8. Xu C, Xie N, Su Y, Sun Z, Liang Y, Zhang N, Liu D, Jia S, Xing X, Han L, Li G, Tong T, Chen J. HnRNP F/H associate with hTERC and telomerase holoenzyme to modulate telomerase function and promote cell proliferation. Cell Death Differ. 2020; 27:1998–2013. https://doi.org/10.1038/s41418-019-0483-6 [PubMed]
9. Jackson JG, Pant V, Li Q, Chang LL, Quintás-Cardama A, Garza D, Tavana O, Yang P, Manshouri T, Li Y, El-Naggar AK, Lozano G. p53-mediated senescence impairs the apoptotic response to chemotherapy and clinical outcome in breast cancer. Cancer Cell. 2012; 21:793–806. https://doi.org/10.1016/j.ccr.2012.04.027 [PubMed]
10. Coppé JP, Patil CK, Rodier F, Sun Y, Muñoz DP, Goldstein J, Nelson PS, Desprez PY, Campisi J. Senescence-associated secretory phenotypes reveal cell-nonautonomous functions of oncogenic RAS and the p53 tumor suppressor. PLoS Biol. 2008; 6:2853–68. https://doi.org/10.1371/journal.pbio.0060301 [PubMed]
11. Calcinotto A, Kohli J, Zagato E, Pellegrini L, Demaria M, Alimonti A. Cellular Senescence: Aging, Cancer, and Injury. Physiol Rev. 2019; 99:1047–78. https://doi.org/10.1152/physrev.00020.2018 [PubMed]
12. Saul D, Kosinsky RL, Atkinson EJ, Doolittle ML, Zhang X, LeBrasseur NK, Pignolo RJ, Robbins PD, Niedernhofer LJ, Ikeno Y, Jurk D, Passos JF, Hickson LJ, et al. A new gene set identifies senescent cells and predicts senescence-associated pathways across tissues. Nat Commun. 2022; 13:4827. https://doi.org/10.1038/s41467-022-32552-1 [PubMed]
13. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, Hackl H, Trajanoski Z. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep. 2017; 18:248–62. https://doi.org/10.1016/j.celrep.2016.12.019 [PubMed]
14. Walter V, Yin X, Wilkerson MD, Cabanski CR, Zhao N, Du Y, Ang MK, Hayward MC, Salazar AH, Hoadley KA, Fritchie K, Sailey CJ, Weissler MC, et al. Molecular subtypes in head and neck cancer exhibit distinct patterns of chromosomal gain and loss of canonical cancer genes. PLoS One. 2013; 8:e56823. https://doi.org/10.1371/journal.pone.0056823 [PubMed]
15. Kuilman T, Michaloglou C, Vredeveld LC, Douma S, van Doorn R, Desmet CJ, Aarden LA, Mooi WJ, Peeper DS. Oncogene-induced senescence relayed by an interleukin-dependent inflammatory network. Cell. 2008; 133:1019–31. https://doi.org/10.1016/j.cell.2008.03.039 [PubMed]
16. Acosta JC, O’Loghlen A, Banito A, Guijarro MV, Augert A, Raguz S, Fumagalli M, Da Costa M, Brown C, Popov N, Takatsu Y, Melamed J, d’Adda di Fagagna F, et al. Chemokine signaling via the CXCR2 receptor reinforces senescence. Cell. 2008; 133:1006–18. https://doi.org/10.1016/j.cell.2008.03.038 [PubMed]
17. Faget DV, Ren Q, Stewart SA. Unmasking senescence: context-dependent effects of SASP in cancer. Nat Rev Cancer. 2019; 19:439–53. https://doi.org/10.1038/s41568-019-0156-2 [PubMed]
18. Coppé JP, Desprez PY, Krtolica A, Campisi J. The senescence-associated secretory phenotype: the dark side of tumor suppression. Annu Rev Pathol. 2010; 5:99–118. https://doi.org/10.1146/annurev-pathol-121808-102144 [PubMed]
19. He Q, Au B, Kulkarni M, Shen Y, Lim KJ, Maimaiti J, Wong CK, Luijten MN, Chong HC, Lim EH, Rancati G, Sinha I, Fu Z, et al. Chromosomal instability-induced senescence potentiates cell non-autonomous tumourigenic effects. Oncogenesis. 2018; 7:62. https://doi.org/10.1038/s41389-018-0072-4 [PubMed]
20. Xie W, Kagiampakis I, Pan L, Zhang YW, Murphy L, Tao Y, Kong X, Kang B, Xia L, Carvalho FL, Sen S, Chiu Yen RW, Zahnow CA, et al. DNA Methylation Patterns Separate Senescence from Transformation Potential and Indicate Cancer Risk. Cancer Cell. 2018; 33:309–21.e5. https://doi.org/10.1016/j.ccell.2018.01.008 [PubMed]
21. Cruickshanks HA, McBryan T, Nelson DM, Vanderkraats ND, Shah PP, van Tuyn J, Singh Rai T, Brock C, Donahue G, Dunican DS, Drotar ME, Meehan RR, Edwards JR, et al. Senescent cells harbour features of the cancer epigenome. Nat Cell Biol. 2013; 15:1495–506. https://doi.org/10.1038/ncb2879 [PubMed]
22. Takasugi M, Yoshida Y, Hara E, Ohtani N. The role of cellular senescence and SASP in tumour microenvironment. FEBS J. 2023; 290:1348–61. https://doi.org/10.1111/febs.16381 [PubMed]
23. Mangelsdorf DJ, Thummel C, Beato M, Herrlich P, Schütz G, Umesono K, Blumberg B, Kastner P, Mark M, Chambon P, Evans RM. The nuclear receptor superfamily: the second decade. Cell. 1995; 83:835–9. https://doi.org/10.1016/0092-8674(95)90199-x [PubMed]
24. Banno A, Reddy AT, Lakshmi SP, Reddy RC. PPARs: Key Regulators of Airway Inflammation and Potential Therapeutic Targets in Asthma. Nucl Receptor Res. 2018; 5:101306. https://doi.org/10.11131/2018/101306 [PubMed]
25. Oh YM, Mahar M, Ewan EE, Leahy KM, Zhao G, Cavalli V. Epigenetic regulator UHRF1 inactivates REST and growth suppressor gene expression via DNA methylation to promote axon regeneration. Proc Natl Acad Sci USA. 2018; 115:E12417–26. https://doi.org/10.1073/pnas.1812518115 [PubMed]
26. Siltari A, Syvälä H, Lou YR, Gao Y, Murtola TJ. Role of Lipids and Lipid Metabolism in Prostate Cancer Progression and the Tumor’s Immune Environment. Cancers (Basel). 2022; 14:4293. https://doi.org/10.3390/cancers14174293 [PubMed]
27. Ma S, Zhou B, Yang Q, Pan Y, Yang W, Freedland SJ, Ding LW, Freeman MR, Breunig JJ, Bhowmick NA, Pan J, Koeffler HP, Lin DC. A Transcriptional Regulatory Loop of Master Regulator Transcription Factors, PPARG, and Fatty Acid Synthesis Promotes Esophageal Adenocarcinoma. Cancer Res. 2021; 81:1216–29. https://doi.org/10.1158/0008-5472.CAN-20-0652 [PubMed]
28. Trayssac M, Hannun YA, Obeid LM. Role of sphingolipids in senescence: implication in aging and age-related diseases. J Clin Invest. 2018; 128:2702–12. https://doi.org/10.1172/JCI97949 [PubMed]
29. Huang T, Wang L, Liu D, Li P, Xiong H, Zhuang L, Sun L, Yuan X, Qiu H. FGF7/FGFR2 signal promotes invasion and migration in human gastric cancer through upregulation of thrombospondin-1. Int J Oncol. 2017; 50:1501–12. https://doi.org/10.3892/ijo.2017.3927 [PubMed]
30. Shang A, Zhou C, Bian G, Chen W, Lu W, Wang W, Li D. miR-381-3p restrains cervical cancer progression by downregulating FGF7. J Cell Biochem. 2019; 120:778–89. https://doi.org/10.1002/jcb.27438 [PubMed]
31. González-Chavarría I, Fernandez E, Gutierrez N, González-Horta EE, Sandoval F, Cifuentes P, Castillo C, Cerro R, Sanchez O, Toledo JR. LOX-1 activation by oxLDL triggers an epithelial mesenchymal transition and promotes tumorigenic potential in prostate cancer cells. Cancer Lett. 2018; 414:34–43. https://doi.org/10.1016/j.canlet.2017.10.035 [PubMed]
32. Lee YI, Choi S, Roh WS, Lee JH, Kim TG. Cellular Senescence and Inflammaging in the Skin Microenvironment. Int J Mol Sci. 2021; 22:3849. https://doi.org/10.3390/ijms22083849 [PubMed]
33. Hwang HJ, Lee YR, Kang D, Lee HC, Seo HR, Ryu JK, Kim YN, Ko YG, Park HJ, Lee JS. Endothelial cells under therapy-induced senescence secrete CXCL11, which increases aggressiveness of breast cancer cells. Cancer Lett. 2020; 490:100–10. https://doi.org/10.1016/j.canlet.2020.06.019 [PubMed]
34. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, Sabedot TS, Malta TM, Pagnotta SM, Castiglioni I, Ceccarelli M, Bontempi G, Noushmehr H. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016; 44:e71. https://doi.org/10.1093/nar/gkv1507 [PubMed]
35. Davis S, Meltzer PS. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics. 2007; 23:1846–7. https://doi.org/10.1093/bioinformatics/btm254 [PubMed]
36. Tacutu R, Thornton D, Johnson E, Budovsky A, Barardo D, Craig T, Diana E, Lehmann G, Toren D, Wang J, Fraifeld VE, de Magalhães JP. Human Ageing Genomic Resources: new and updated databases. Nucleic Acids Res. 2018; 46:D1083–90. https://doi.org/10.1093/nar/gkx1042 [PubMed]
37. Bhattacharya S, Andorf S, Gomes L, Dunn P, Schaefer H, Pontius J, Berger P, Desborough V, Smith T, Campbell J, Thomson E, Monteiro R, Guimaraes P, et al. ImmPort: disseminating data to the public for the future of immunology. Immunol Res. 2014; 58:234–9. https://doi.org/10.1007/s12026-014-8516-1 [PubMed]
38. Wen Y, Guo G, Yang L, Chen L, Zhao D, He X, Zhang R, Huang Z, Wang G, Zhang L. A tumor microenvironment gene set-Based prognostic signature for non-small-cell lung cancer. Front Mol Biosci. 2022; 9:849108. https://doi.org/10.3389/fmolb.2022.849108 [PubMed]
39. Hu K, Xiang J, Yu YX, Tang L, Xiang Q, Li JM, Tang YH, Chen YJ, Zhang Y. Significance-based multi-scale method for network community detection and its application in disease-gene prediction. PLoS One. 2020; 15:e0227244. https://doi.org/10.1371/journal.pone.0227244 [PubMed]
40. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010; 26:1572–3. https://doi.org/10.1093/bioinformatics/btq170 [PubMed]
41. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14:7. https://doi.org/10.1186/1471-2105-14-7 [PubMed]
42. Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw. 2010; 33:1–22. https://doi.org/10.18637/jss.v033.i01 [PubMed]
43. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015; 12:453–7. https://doi.org/10.1038/nmeth.3337 [PubMed]
44. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
45. 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, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018; 24:1550–8. https://doi.org/10.1038/s41591-018-0136-1 [PubMed]
46. Sun D, Wang J, Han Y, Dong X, Ge J, Zheng R, Shi X, Wang B, Li Z, Ren P, Sun L, Yan Y, Zhang P, et al. TISCH: a comprehensive web resource enabling interactive single-cell transcriptome visualization of tumor microenvironment. Nucleic Acids Res. 2021; 49:D1420–30. https://doi.org/10.1093/nar/gkaa1020 [PubMed]