Abstract Background Head and neck squamous cell carcinoma(HNSCC) is the sixth most common malignancy worldwide, with more than 890,000 new cases and 450,000 deaths annually. Its major risk factors include smoking, alcohol abuse, aging, and poor oral hygiene. Due to the lack of early and effective detection and screening methods, many patients are diagnosed at advanced stages with a five-year survival rate of less than 50%. In this study, we deeply explored the expression of Aging-related genes(ARGs) in HNSCC and analyzed their prognostic significance using single-cell sequencing and spatial transcriptomics analysis. This research aims to provide new theoretical support and directions for personalized treatment. Annually, more than 890,000 new cases of head and neck squamous cell carcinoma (HNSCC) are diagnosed globally, leading to 450,000 deaths, making it the sixth most common malignancy worldwide. The primary risk factors for HNSCC include smoking, alcohol abuse, aging, and poor oral hygiene. Many patients are diagnosed at advanced stages due to the absence of early and effective detection and screening methods, resulting in a five-year survival rate of less than 50%. In this research, single cell sequencing and spatial transcriptome analysis were used to investigate the expression of Aging-related genes (ARGs) in HNSCC and to analyse their prognostic significance. This research aims to provide new theoretical support and directions for personalized treatment. Methods In this study, we investigated the association between HNSCC and AGRs by utilizing the [48]GSE139324 series in the GEO database alongside the TCGA database, combined with single-cell sequencing and spatial transcriptomics analysis. The data were analyzed using Seurat and tSNE tools to reveal intercellular communication networks. For the spatial transcriptome data, SCTransform and RunPCA were applied to examine the metabolic activities of the cells. Gene expression differences were determined through spacerxr and RCTD tools, while the limma package was employed to identify differentially expressed genes and to predict recurrence rates using Cox regression analysis and column line plots. These findings underscore the potential importance of molecular classification, prognostic assessment, and personalized treatment of HNSCC. Results This study utilized HNSCC single-cell sequencing data to highlight the significance of ARGs in the onset and prognosis of HNSCC. It revealed that the proportion of monocytes and macrophages increased, while the proportion of B cells decreased. Notably, high expression of the APOE gene in monocytes was closely associated with patient prognosis. Additionally, a Cox regression model was developed based on GSTP1 and age to provide personalized prediction tools for clinical use in predicting patient survival. Conclusions We utilized single-cell sequencing and spatial transcriptomics to explore the cellular characteristics of HNSCC and its interaction with the tumor microenvironment. Our findings reveal that HNSCC tissues show increased mononuclear cells and demonstrate enhanced activity in ARGs, thereby advancing our understanding of HNSCC development mechanisms. Supplementary Information The online version contains supplementary material available at 10.1007/s12672-024-01672-z. Keywords: Spatial transcriptomics, Single-cell RNA sequening, Head and neck cancer, Pathogenesis, Aging-related genes, Integrative analysis, Immunotherapy, Head and neck squamous cell carcinoma, Immune microenvironment, Prognosis, Single-cell analysis Introduction Head and neck squamous cell carcinoma (HNSCC) is a group of malignant tumors originating from epithelial cells lining the mucous membranes of the upper airways and food passages (oral cavity, oropharynx, larynx, or hypopharynx), as well as a heterogeneous collection of malignant tumors of the salivary glands, upper gastrointestinal tract, and thyroid gland [[49]1–[50]3].HNSCC as the most common and deadly malignant tumor of the head and neck [[51]4] and also the sixth most common malignant tumor in the world. According to the latest global cancer statistics (GLOBOCAN 2022), there are more than 890,000 new cases and 450,000 deaths annually from head and neck cancers, including tumors of the larynx, lips, mouth, nasopharynx, oropharynx, and hypopharynx [[52]2, [53]5].Its classical causative risk factors include smoking and alcohol abuse, in addition to aging, poor oral hygiene and dietary deficiencies in vegetables as risk factors. It is interesting to note that betel nut product chewing is a known risk factor for the high prevalence of oral cancer in certain geographical areas [[54]6, [55]7].Because there is no effective early detection or screening method for HNSCC, careful physical examination is an important tool for early diagnosis, but its sensitivity is not always ideal [[56]8].Most patients have advanced disease at the time of diagnosis and have no clinical history of precancerous lesions [[57]9].This usually means that the disease has become locally aggressive and/or has spread to distant organs, leading to increased difficulty in treatment while predicting a poor prognosis for the patient, with the patient's five-year survival rate being less than 50% [[58]10].Targeted immunotherapy has provided survival benefit to patients with HNSCC, but fewer than 20% of patients achieve a durable response, making it critical to deepen our understanding of its pathogenesis in order to develop more effective therapeutic regimens [[59]11]. Aging is a multifactorial, time-dependent biological process and a time-dependent functional decline in most organisms [[60]12].Aging remains the most important risk factor for a variety of cancers, and its prevalence increases dramatically with age, peaking at around 85 years of age [[61]13].Cancer and aging may have a share a common origin, and can be considered as the same two different manifestations of the underlying process, the accumulation of cellular damage [[62]14]. In 2013 Guido Kroemer et al. proposed nine hallmarks of aging, four of which, genomic instability, epigenetic alterations, chronic inflammation and dysbiosis, are similar to those of Hanahan and Weinberg [[63]15] in 2011 who proposed specific cancer hallmarks. However, the two senescence hallmarks of telomere attrition and stem cell exhaustion during aging may play an inhibitory role in tumorigenesis [[64]16]. The relationship between aging and cancer is complex, and the closely linked for the developmental significance of cancer [[65]17], but this link is not fully clear in the mechanism of HNSCC occurrence and development. Future studies need to synthesize the multifaceted features of aging, as well as understand the interactions between HNSCC and aging, and provide new strategies for prevention and treatment. Spatial transcriptomics, as an emerging bioinformatics approach, breaks through the limitations of traditional transcriptomics and permits to explore the heterogeneity of gene expression in cells while maintaining their original spatial location information [[66]18]. The present study synthesized the use of multi-omics analyses, including single-cell sequencing and spatial transcriptomics, to utilize ARGs to predict the prognosis of patients with HNSCC and to reveal the HNSCC-senescence the association between HNSCC and aging was revealed. Based on the comprehensive multi-omics analysis, we gained an in-depth understanding of the thematic mechanisms of ARGs in HNSCC, with the aim of providing new insights into the molecular classification, prognostic assessment and personalized treatment of head and neck cancer. Method Raw data sources In a study aimed at performing single-cell sequencing analysis of head and neck squamous cell carcinoma (HNSCC), we accessed the [67]GSE139324 series [[68]19, [69]20], a gene expression dataset related to HNSCC, from the Gene Expression Omnibus (GEO) database ([70]https://www.ncbi.nlm.nih.gov/geo/) [[71]21]. From this dataset, we selected gene expression data from six samples for analysis, comprising three normal tissue samples from healthy individuals and three diseased tissue samples from HNSCC patients. This comparative analysis strategy enabled us to explore gene expression changes in cells during the disease state and to understand how these changes affect cell function and tumor development. The HNSCC transcriptome data were obtained from The Cancer Genome Atlas (TCGA) database ([72]https://portal.gdc.cancer.gov/) [[73]22, [74]23], which included a total of 548 samples from HNSCC patients, with 44 paraneoplastic tissue samples serving as controls for normal tissue. Additionally, we retrieved ARGs from the Human Aging Genomic Resources (HAGR) ([75]https://www.genomics.senescence.info/) [[76]24], a database focusing on ARGs (Supplementary Table 1). Single-cell sequencing data processing In our research, the “Seurat” [[77]25] software package is one of the most widely used and important analysis tools in the field of single-cell research. Its functions include data import, filtering, normalization, feature selection, scaling, dimensionality reduction, and clustering. After obtaining the samples from the GEO database, we first normalized the data using “Seurat” in the R package to perform dimensionality reduction and clustering. To highlight gene variability, we filtered the data again and selected 2000 highly variable genes (Top 2000 genes with differences between HNSCC samples and healthy samples). Subsequently, we used “RunPCA” for principal component analysis to reduce dimensionality. In terms of data integration and clustering, tSNE helps reveal the structure of cell communities in single-cell RNA sequencing data by showing the distribution of different cell populations in the reduced dimensionality space. By examining tSNE plots, potential cell subsets, state transitions, or continuums of transcriptional extent can be identified, enhancing our understanding of the data's structure and patterns. We visualized these structures using tSNE and annotated cell types based on literature. To increase the depth and breadth of single-cell RNA sequencing data analysis, we employed five different scoring algorithms (AddModuleScore, GSVA, AUCell, Ucells, and singscore) to provide a comprehensive and multidimensional interpretation of the data, thus improving the study's credibility and accuracy [[78]26, [79]27]. After scoring, we continued to identify genes that differed significantly between tissue types. Cell trajectory analysis involves tracking the dynamic changes of cells across various developmental or disease states using single-cell transcriptome data. The key steps include identifying distinct cellular states, typically visualized through dimensionality reduction techniques such as t-SNE or UMAP. Subsequently, the Monocle algorithm is employed to infer the trajectories of cells transitioning from one state to another. Finally, these trajectories are categorized into different phases to reflect the varying states of cells during development or disease progression. Additionally, “CellCall” is a crucial toolkit for identifying intercellular communication networks and internal regulatory signals by combining ligand/receptor expression with downstream transcription factor (TF) activity of certain ligand-receptor (LR) pairs [[80]28, [81]29]. We utilized the “CellCall” R software package to reveal integrated intercellular communication networks. Spatial transcriptome sequencing data processing First, we processed and analyzed the spatial transcriptome data using the “Seurat” R package. Subsequently, we normalized and denoised the data with the “SCTransform” program [[82]30, [83]31] and filtered out genes with non-zero expression in more than 10 spatial locations. Next, we downscaled the data using the “RunPCA” method. Finally, we utilized the “SpatialDimPlot” function to visualize cell clustering and dimensionality reduction analysis [[84]32]. Additionally, we utilized the “scMetabolism” [[85]33] R software package along with the AUCell method to estimate cellular metabolic activity. The VISION algorithm was employed to synthesize multiple single-cell data types for comprehensive cellular status assessment. These in-depth metabolic analyses revealed a wide variety of cellular functions and significant differences, thereby providing new insights into the pathogenesis of tumors. To enhance the depth and accuracy of single-cell analysis, we implemented a series of complementary steps using the “Scanpy” [[86]34] suite in Python. These steps included data preprocessing, dimensionality reduction, image processing, cluster analysis, visualization of results, exploration of cell developmental trajectories, and identification of associated genes. Subsequently, we integrated single-cell RNA sequencing data with spatial location information using the “stLearn” [[87]35] library in Python to perform spatial image reconstruction, cell type identification, spatial gene expression pattern analysis, and the discovery of cellular communication signals within tissues. From the collected data, we accurately categorized cell subpopulations, providing valuable insights for the in-depth exploration of cell functions and their interactions. Combining spatial transcriptome data and single-cell sequencing data for deconvolution analysis Based on bulk RNA-seq data, deconvolution techniques can estimate the proportion of each cell type in a tissue sample, thereby revealing cellular heterogeneity. Single-cell sequencing technology further enhances these findings by showing how specific cell types and their subpopulations respond to environmental factors. Spatial transcriptomics adds another dimension by analyzing RNA-seq data at the spatial level, providing a precise geographic map of cell distribution. To further improve the accuracy of deconvolution analysis, we integrated single-cell and spatial transcriptomics datasets.This synergistic analysis not only deepened our understanding of cellular heterogeneity but also provided unprecedented accuracy in mapping cellular organization. Using “spacexr” ([88]https://fsf.org/) [[89]36] an R software package for RCTD deconvolution analysis, we delved deeper into the variability of gene expression in tumor cells. This approach enabled us to construct a detailed, annotated, and benchmarked single-cell transcriptomics dataset, which also helped us to accurately distinguish and classify cell populations within tumors. Subsequently, we transformed the spatial transcriptome data into SpatialRNA structures for subsequent RCTD deconvolution analysis. This step facilitated the application of the least squares method to accurately infer gene expression distribution and the proportion of each cellular component in the spatial context. Furthermore, RCTD deconvolution analysis, facilitated by the “mistyR” [[90]37] R package, forms the foundation for investigating densely spatial interactions among cell clusters in tissues. With “mistyR”, we meticulously explored potential interactions among cell subpopulations in tissues based on their spatial distribution and known gene expression patterns. This analytical approach not only enhances our understanding of cellular interactions through spatial proximity but also unveils novel pathways for investigating complex mechanisms of intercellular communication. HNSCC monocyte senescence-related key genes combined with bulk data for prognostic and immunotherapy analysis This study explores the differential expression of ARGs between normal and tumor tissues and comprehensively assesses their prognostic implications in HNSCC. The aim of this research is to elucidate the role of these genes in tumor development and investigate their potential as therapeutic targets. We analyzed HNSCC patient samples with complete clinical and transcriptomic data sourced from the TCGA database. Using the R package “limma” [[91]38, [92]39], specifically designed for microarray data, we identified significantly differentially expressed genes between HNSCC and normal samples (adjusted P < 0.05 and logFC > 1). Subsequently, we conducted univariate and multivariate Cox regression analyses on these genes, and presented the results using forest plots generated by the “forestplot” package [[93]40, [94]41]. Based on the findings from multivariate Cox proportional hazards analysis, we constructed column plots to predict 1-, 3-, and 5-year recurrence rates in HNSCC, incorporating prognostic variables with adjusted P < 0.05 using the “rms” package ([95]https://hbiostat.org/r/rms/). Integrating multiple prognostic factors into a scoring system can enhance the precision of clinical decision-making and optimize patient management strategies. Statistical analysis Statistical analysis was conducted using RStudio (R version 4.3.2) and its associated libraries, along with the PyCharm integrated development environment for Python. For comparisons of continuous variables, we employed the nonparametric Wilcoxon rank sum test, and assessed correlations between variables using Pearson correlation analysis. In this study, a significance level of P < 0.05 was established to determine statistical significance. Results were considered statistically significant only if they met this criterion. Result HNSCC single cell sequencing standard preprocessing and cell annotation Figure [96]1 outlines the specific workflow of our study. We analyzed data from six single-cell samples in the [97]GSE139324 dataset, comprising three samples from HNSCC patients and three healthy control samples. Post-”harmony” integration revealed a tendency for HNSCC samples to cluster together, distinct from the CT samples. PCA plot distributions further highlighted the similarities and differences between these groups (Fig. [98]2A). To achieve more precise results in cell type identification and functional analysis, we selected the top 2000 highly variable genes (Top 2000 genes with differences between HNSCC samples and healthy samples) for subsequent analysis using the FindVariableFeatures function in the Seurat package. Nineteen unique cell clusters were then identified using UMAP visualization (Fig. [99]2B). Subsequently, we utilized classical marker genes, combined with UMAP clustering and marker gene expression bubble maps, to accurately annotate the cell clusters manually (Fig. [100]2C). Furthermore, we visualized the distribution of cell subpopulations in 2D and 3D tSNE space to illustrate their heterogeneity ((Fig. [101]2D, [102]E). Our analysis demonstrated a significant remodeling of the immune cell composition in the tumor microenvironment of HNSCC patients, as evidenced by a significant increase in the proportion of monocytes and macrophages and a significant decrease in the proportion of B cells and CD4 + T cells, especially CD8 + T cells, in HNSCC patients compared with normal tissues. Meanwhile, the increase in monocytes and macrophages suggested that tumor-associated macrophages (TAMs) played a key role in the pro-inflammatory and immunosuppressive effects of the tumor [[103]42]. The results of this study not only provide a basis for further understanding of the immune escape mechanism of HNSCC, but also provide a new potential direction for the development of immunotherapeutic strategies targeting macrophages or monocytes. (Fig. [104]2F). Lastly, we performed pathway enrichment analysis on the heatmaps of each subpopulation using the “ClusterGVis” R software package. This analysis uncovered significant enrichment in gene expression related to B-cell-mediated immune responses, cell killing, regulation of natural killer cell activation, and immune-response activation signaling pathways (Fig. [105]2G). These findings underscore the potential importance of immune cells in the development of HNSCC. Fig. 1. [106]Fig. 1 [107]Open in a new tab Flowchart of the main design of this study. First, raw data were obtained from the HAGR and GEO databases and cells were categorized using a dimensionality-reducing clustering method to identify nine major single-cell types. Subsequently, single-cell functional analysis was performed to explore the biological functions of each cell type. In addition, bulk analysis was performed to combine the overall data for prognostic analysis, to demonstrate the survival prediction model and its associated results, and then to assess the correlation between the expression level of APOE + monocytes and patient survival using Kaplan–Meier survival curves, as well as to analyze the correlation between the expression level of APOE_monocyte and specific immune-related factors. Further, the expression profiles of monocyte senescence-related genes were analyzed by spatial transcriptome sequencing, followed by in-depth analysis of cell–cell interactions. Finally, RCTD deconvolution technique was applied to analyze spatial cell developmental trajectories and investigate the role of cellular communication in prognosis Fig. 2. [108]Fig. 2 [109]Open in a new tab Integration of single-cell sequencing and senescence-associated genomic screening techniques to identify HNSCC-associated signature genes. A The left panel shows the gene screening of the single-cell expression matrix and the linear downscaling achieved through PCA clustering of the samples. The right panel shows the distribution of the top 20 principal components (PCs). B Nonlinear dimensionality reduction using UMAP was employed to cluster all single-cell data into 19 different cell clusters. C To illustrate the heatmap of the expression levels and average expression values of key genes in different cell populations, we manually annotated the immune cell subpopulations based on known marker genes to identify cell subpopulations. D We meticulously manually annotated the twenty cell clusters, and these subpopulations were then categorized into nine different cell types and depicted in a three-dimensional spatial context using tSNE. E This manual labeling process was further extended by performing dimensionality reduction clustering using UMAP, which ultimately categorized these 19 cell subpopulations into nine cell types. F The distribution of cell proportions of the nine identified cell types in various samples is shown. G Enrichment analysis was performed on the nine cell types obtained from the annotation. The leftmost line graph depicts the transition from CD4 + T cells to macrophages arranged from top to bottom. The accompanying heatmap provides the gene expression patterns associated with each cell type, while the rightmost section elucidates the functional pathways enriched in each cell type, including their corresponding highly expressed genes Assessment of senescence characteristics of HNSCC cells using single-cell scoring of ARGs Aging refers to the progressive decline in the function of organs and tissues over time [[110]43]. The aging process is evident at both the molecular and cellular levels. Cellular senescence is defined as the irreversible cessation of cell division that occurs after a limited number of divisions. Aging-related genes (ARGs) play a regulatory role in this process and are closely associated with the development of various cancers [[111]44]. To further assess the senescence characteristics in cell or tissue samples and to elucidate the potential mechanisms of ARGs in HNSCC. First, based on the set of ARGs and applying five well-recognized algorithms—AUCell, UCell, singscore, AddModuleScore, and a composite model synthesizing these algorithms—were used -to assess the senescence-related characteristics of immune cell subpopulations. By visualizing and analyzing UMAP plots, we observed that monocyte presented significantly higher scores in the HNSCC group when comparing HNSCC with CT. (Fig. [112]3A, [113]B) Combining the results of the five scoring algorithms scoring bubble plots and violin plots analysis, especially in monocytes, the transcriptional activity of ARGs was significantly higher in HNSCC (Fig. [114]3C, [115]D). Notably, the transcriptional profiles of the monocytes were lower than those of the HNSCC subgroups, except for the monocytes in the disease group with enhanced transcription (Fig. [116]3E). To further investigate the role and mechanism of monocytes in HNSCC, we obtained a monocyte-focused single-cell profile in by isolating other immune subgroups. Then, differential analysis was performed to obtain ARGs with differential expression (logFC > 1, P < 0.05). The violin plot showed the top nine differential genes with the largest differences, and the results showed enhanced expression in APOE, C1QA, CEBPB, and LMNA, and attenuated expression in HSPA1B, HSPA1A, and PLCG2 in monocytes from the HNSCC group (Fig. [117]3F). Finally, the FEATURE PLOT shows the distribution of the expression of these nine genes in various immune subpopulations (Fig. [118]3G). Fig. 3. [119]Fig. 3 [120]Open in a new tab Cell type identification and gene expression analysis of single-cell RNA sequencing data. A UMAP downscaled visualization plots of CT and HNSCC tissues illustrate the distribution of different cell types. B Cells were relabeled and visualized using a scoring metric called Scoring. C Multiple violin plots depict the distribution of scores of different cell types under the five single-cell scoring metrics: AUCell, UCell, singscore, ssgsea, and Add, as well as the composite score. D Correlation analysis is visualized by bubble plots. These plots depict positive and negative correlations between the five single-cell scoring methods and the overall scores for the different cell types. The size of the bubbles corresponds to the strength of the correlation, transitioning from blue (negative) to red (positive). E ARGs were introduced, highlighting their expression differences in different cell types in the HNSCC and normal groups. The normal group is shown in blue and the HNSCC group is shown in yellow. F Violin plots illustrate ARGs, showcasing their significant expression differences in monocytes of HNSCC group and normal group. *P < 0.05,**P < 0.01,***P < 0.001. G A densitogram illustrates ARGs expression with high relative cell density shown in dark Pseudotemporal analysis and intercellular communication analysis Cell trajectory analysis provided valuable insights into developmental trajectories, tumor immune cell dynamics, and cell differentiation relationships at single-cell resolution [[121]45]. The heatmap demonstrated the expression of ARGs at various time points throughout the process, particularly highlighting the high expression of genes such as HIF1A, GRN, and SOD2 during the neoplastic phase of the tumor (Fig. [122]4A). To further explore the sequential changes of cells during development or biological processes and to determine the developmental trajectories or state transitions of cells, we performed pseudotemporal and cellular trajectory analysis on monocyte subpopulations. We identified two subpopulations (Fig. [123]4B), and as time progressed, the cluster predominant in normal samples (cluster 11) evolved into a cluster abundant in HNSCC samples (cluster 17) (Fig. [124]4C, [125]D). By analyzing the interactions between monocytes and other cells, we classified the immune cells into two subgroups based on the expression of ARGs: the group with high expression of ARGs (aginghigh) and the group with low expression of ARGs (aginglow). The bubble diagram shows that the TNF signaling pathway exhibits higher activity in the aging^high^ group (Fig. [126]4E). This may indicate that increased TNF signaling pathway activity is associated with an enhanced inflammatory response in the tumor microenvironment [[127]46]. Inflammation is a crucial component of the tumor microenvironment, closely correlated with tumorigenesis, progression, and metastasis [[128]47]. Using the circle plot (Fig. [129]4F), we examined the communication between immune cells in these two subpopulations. Additionally, we investigated ligand-receptor interactions through heatmap analysis (Fig. [130]4G) and found significant differences in IL15-IL2RA, IL15-IL2RB, and IL15-IL2RG production by Treg cells. Finally, to explore the ligand-receptor-transcription factor connection, we constructed Sankey diagrams (Fig. [131]4H). Fig. 4. [132]Fig. 4 [133]Open in a new tab Pseudotemporal analysis and intercellular communication analysis. A The developmental heatmap of ARGs in monocytes shows time progression from left to right and expression levels from low (blue) to high (red). B A comparison plot of non-linear descending clustering of UMAP in monocytes from HNSCC group and normal group is presented. C The plot of tumor monocyte developmental trajectories allows for the observation of developmental trajectories of different monocyte populations in the tumor by combining panels A and C. D The developmental map of monocyte trajectories in normal and tumor groups combines panels B and C to show the developmental trajectories of different monocyte populations in these groups. E A correlation bubble diagram depicts the enrichment analysis of different signaling pathways in different cell populations; the size of the bubbles represents the magnitude of the correlation, and the color from blue to red indicates the correlation from negative to positive. F Correlation circle plots illustrate interactions between different cell types. G A heat map shows cell–cell ligand-receptor correlations. H A Sankey diagram visualizes ligand-receptor-transcription factor interactions between monocytes and Treg cells with higher expression of ARGs Spatial transcriptome sequencing unites genes characterizing monocyte senescence To delve into the spatial transcriptome characterization of HNSCC patients, we acquired spatial transcriptome data from a head and neck cancer patient, derived from a classic literature source [[134]48]. To ensure the integrity and reliability of the data, we implemented a series of stringent quality control procedures. Through these QC steps, we excluded interference from mitochondrial and ribosomal genes, ensuring the accuracy of the subsequent analysis. Next, we employed the SCTransform method to normalize and depth-correct the data, reducing the effect of technical differences. This processing flow revealed thirteen spatially distributed cell clusters, and through the UMAP projection technique, we obtained thirteen cell clusters and projected them onto spatial slices, allowing us to clearly observe their distribution (Fig. [135]5A, [136]B). We investigated the spatial distribution of oxidative phosphorylation processes in HNSCC (Fig. [137]5C). At the same time, we utilized “scMetabolism”, an R software package, to analyze the activity of metabolic pathways in the different cell subtypes, examining the impact of cellular senescence on the pathological process of HNSCC (Fig. [138]5D). Through the analysis of metabolic pathways in different cell subtypes, we found that cellular senescence not only affects the metabolic activity of tumor cells, but also plays an important role in the metabolic regulation in immune cells, especially monocytes. Specifically, the senescence process may impair immune surveillance and promote tumor progression by altering the metabolic status of monocytes and affecting their function in the tumor microenvironment. Finally, the ARGs characteristics of HNSCC monocytes identified in the single-cell analyses were integrated into these spatial clusters to form a clear bubble map (Fig. [139]5E), thus revealing the expression patterns of ARGs in different spatial clusters. Fig. 5. [140]Fig. 5 [141]Open in a new tab Spatial transcriptome sequencing associated with monocyte-associated glycosyltransferase genes. A Spatial transcriptome data were divided into 13 cell clusters after UMAP nonlinear dimensionality reduction clustering. B Projections of spatial transcriptome tissue sections into 13 cell clusters were made after dimensionality reduction clustering. C A visualized expression map of the oxygenate phosphorylation metabolic pathway was created in spatial transcriptome tissue sections, with selection criteria referenced from Figure D. D Correlation bubble plots show the expression of different functional pathways in the 13 Seurat clusters. E The bubble sizes of monocyte ARGs expression in the 13 cell clusters represent the correlation size, with color ranging from blue (negative) to red (positive) RCTD deconvolution analysis and spatial cell development analysis The RCTD deconvolution technique and single-cell data analysis were integrated to provide comprehensive insights into the spatial localization of genes, enabling the inference of cellular subpopulation distributions in heterogeneous samples. Initially, RCTD deconvolution analysis was applied to each tissue section location, enabling quantification of ARGs expression levels in various monocytes. Based on these findings, monocytes were classified into two categories: aginghigh and aginglow. Following this, we utilized the “Stlearn” and “Scanpy” Python packages to conduct Louvain clustering and data normalization on the spatial transcriptome data, identifying 12 distinct cell clusters (Fig. [142]6A). Fig. 6. [143]Fig. 6 [144]Open in a new tab Deconvolution of RCTD using spatial cell development analysis. A Louvain clustering was employed to classify the spatial transcriptome data into 12 distinct cell clusters. B RCTD deconvolution was utilized to map and correlate cell types from single-cell data to their corresponding cell types in the spatial transcriptome tissue sections. C A clear transition from cluster 0 to cluster 1 was observed and mapped within the spatial transcriptome tissue sections. D The Leiden clustering method further divided the spatial transcriptome data into 16 cell clusters, providing a more refined classification. E A correlation network diagram was constructed to illustrate the interactions among different cell types. F A heatmap was generated to display the correlations between various cell types. G The custom R package “mistyR” was employed to analyze and visualize the enhanced statistics and correlation magnitudes of different cell types. H A histogram illustrates the contributions of three custom functions, intra, para_15, and juxta_5, to the measure of cell importance. I–N Three subsequent sets of heatmaps and network diagrams demonstrate the correlation between cell types, as calculated from the intra, para_15, and juxta_5 microenvironments Image analysis indicated that cell cluster 0 predominantly contained aginglow cells, while cell cluster 1 was predominated by aginghigh cells, thereby further validating our classification of cell subpopulations. Furthermore, through analysis of the spatial developmental trajectories of cells, we observed a significant migration from cell cluster 1 to cell cluster 0 within the tumor region (Fig. [145]6B, [146]C). To comprehensively understand the spatial organization of cells, we utilized the Leiden clustering method to delineate 16 more detailed cell populations (Fig. [147]6D). Analysis of intercellular network topology revealed a complex network of interactions among cells (Fig. [148]6E). Notably, the heatmap depicting intercellular correlation strength showed a positive correlation between aginglow cells and B cells, suggesting enhanced B cell activity in HNSCC samples, which may potentially inhibit the development of HNSCC (Fig. [149]6F).Additionally, using the “mistyR” R package, we investigated the spatial interactions between cells and meticulously analyzed the positions of various cell populations (Fig. [150]6G). The histograms depicting cell type contributions clearly illustrate the cellular composition across distinct clustering categories, underscoring the pivotal role of intra and para_15 in elucidating cellular population dynamics and spatial relationships (Fig. [151]6H).Utilizing the intra, juxta_5, and para_15 analysis functions within “mistyR”, we enhanced our comprehension of intricate cell interactions and inter-cellular correlations (Fig. [152]6I–[153]N). SPOTlight deconvolution analysis using spatial pseudotemporal analysis The distribution of tumors and normal tissues in the sections is depicted in Fig. [154]7A. To accurately delineate the spatial expression trajectories of various cellular entities, we employed SPOTlight deconvolution, a technique that enables precise annotation of single-cell data integrated with spatial datasets. This approach revealed the expression dynamics of a diverse array of cell types, including APOE- and APOE + monocytes, B cells, CD4 + T cells, CD8A T cells, endothelial cells, macrophages, mast cells, NK cells, and regulatory T cells (Fig. [155]7B–K). Integration of these datasets highlighted the significantly elevated expression of APOE + monocytes in the tumor group, with more nuanced expression changes observed in other cell types, consistent with our previous single-cell findings. These results underscore the pivotal role of APOE in HNSCC. Fig. 7. [156]Fig. 7 [157]Open in a new tab SPOTlight deconvolution analysis with spatial pseudotemporal analysis. A Louvain clustering demonstrated on spatial transcriptome tissue sections identifies ten cell populations. B–K SPOTlight deconvolution of spatial transcriptome tissue sections visualizes the expression of different cell types derived from single-cell data. L A scatterplot demonstrates the relationship between the average expression level of genes (X-axis) and the empirical dispersion (Y-axis) to identify genes with significant variability in gene expression experiments. M Developmental trajectories of five cell populations expressing high levels of APOE in monocytes are shown. N A proposed temporal progression of the five cell populations is depicted, with the gradual passage of time indicated from black to blue. O Changes in APOE gene expression in cell cluster 0, cell cluster 1, cell cluster 2, cell cluster 4, and cell cluster 7 are observed with pseudotemporal progression Upon thorough analysis, we selected the first 1000 cells with high information content (high variability) for a rigorous time-series approach, aiming to derive accurate and biologically meaningful trajectories (Fig. [158]7L). Specifically, we conducted pseudo-timeseries analysis on tumor tissues encompassing clusters 0, 1, 2, 4, and 7. This examination unveiled a notable transformation process: as pseudotime advanced, cluster 7 bifurcated at a critical node, giving rise to clusters 0 and 4, alongside the emergence of cluster 1 (Fig. [159]7M, [160]N). Significantly, our meticulous investigation into the dynamics of APOE gene expression during pseudotemporal progression revealed a marked increase as cells transitioned from cluster 7 to cluster 1, highlighting the pivotal role of APOE in tumor cell clusters (Fig. [161]7O). Collectively, these findings underscore the pivotal role of cluster 1 in primary tumor formation and suggest that APOE may play a significant role in this tumorigenic mechanism, presenting itself as a potential therapeutic target for future development strategies against HNSCC. Analysis of space cell communications Utilizing the RCTD deconvolution analysis technique, we effectively integrated single-cell datasets with spatial transcriptomics information. This approach accurately mapped the characteristics and spatial distributions of APOE + monocytes, APOE-monocytes, B cells, CD4 + T cells, CD8A T cells, regulatory T cells, and endothelial cells (Fig. [162]8A). Fig. 8. [163]Fig. 8 [164]Open in a new tab Spatial cell interaction analysis. A Visualization and analysis of APOE + monocytes and other cell types through the deconvolution of single-cell data integrated with spatial transcriptome data. B Identification of the top 50 ligand-receptor pairs in APOE + monocytes. C Expression analysis of the highest ranked ligand-receptor pair, S100A9_TLR4, under different conditions, as shown in Figure B. D A network diagram illustrating cell–cell ligand-receptor interactions. E A circle diagram depicting cell–cell interactions Next, we thoroughly analyzed the HNSCC tissue sections using the “Stlearn” package in Python. Specifically, our focus was on identifying the ligand-receptor signals within the tissue spots and assessing their statistical significance through permutation tests. These methods effectively identified the top 50 ligand-receptor pairs and provided each pair with a spatial score to indicate its spatial association (Fig. [165]8B). Through comprehensive analysis of the ligand-receptor interaction between S100A9 and TLR4, we verified through spatial enrichment assessment that this interaction was significantly enriched in the tumor region compared to the surrounding non-tumor region. Specifically, the interaction between S100A9 and TLR4 was notably concentrated in the spatial vicinity of tumor cells, highlighting its pivotal role in intercellular spatial communication. These findings suggest that the S100A9-TLR4 interaction may regulate the development and progression of tumors, potentially influencing tumor growth and metastasis. This discovery provides critical insights for future investigations into the mechanisms underlying spatial cellular interactions within tumors (Fig. [166]8C). Through additional network diagram and chordogram analyses, our focus was on examining the elevated APOE expression levels specifically within tumor tissues. We paid particular attention to the ligand-receptor interaction between S100A9 and TLR4, employing a receptor tandem model to explore this interaction mechanism in depth. Notably, we discovered that during the progression of HNSCC and its response to therapy, APOE + monocytes and CD4 + T cells play a crucial role in regulating and modulating the functions and behaviors of other cell types (Fig. [167]8D, [168]E). Through this comprehensive analysis, we significantly enhanced our understanding of the internal organization of the tumor ecosystem and the spatial interactions among cells, revealing complex and well-organized patterns. APOE + monocytes combined with bulk data for prognostic and immunotherapy target analysis For the investigation of the prognostic impact of APOE + monocytes in HNSCC patients, we utilized the survival + survminer software package to determine the optimal cut-off value for the expression level of APOE + monocytes, which was statistically significant (P < 0.05). Based on this cut-off value, we categorized the patients in our study into high and low expression groups. Through Kaplan–Meier survival curve analysis, we observed that the low-expression group of APOE + monocytes in HNSCC patients exhibited a significantly higher survival rate than the high-expression group (Fig. [169]9A). This finding aligns with our previous hypothesis that high APOE expression is associated with poorer prognosis in HNSCC, highlighting its potential as both a biomarker and therapeutic target. Fig. 9. [170]Fig. 9 [171]Open in a new tab APOE-based prognosis and immunotherapy target analysis. A Kaplan–Meier (K-M) curves illustrating the survival probability of patients with high and low expression of APOE-positive monocytes over time. B Box plots depicting the efficacy of immunotherapy in patients with varying levels of APOE-positive monocyte expression. C–L Scatterplots showing the correlation between APOE-positive monocytes and tumor microenvironment markers To explore the role of APOE in immune checkpoint therapy for HNSCC, we employed the TIDE algorithm, a tool used to predict response to immune checkpoint inhibitors (e.g., PD-1/PD-L1 inhibitors) in cancer patients [[172]49]. Our analysis revealed that high APOE expression in APOE + monocytes showed significant differences in TIDE scores between sham and true factions. These results indicated a poor outcome for the high-expression group when treated with immune checkpoint therapy, consistent with Fig. [173]9A, further emphasizing the strong association between high APOE expression and poor prognosis in HNSCC patients (Fig. [174]9B). Further correlation analysis revealed associations between APOE + monocytes and key markers of the tumor microenvironment. Notably, APOE + monocytes were positively correlated with markers like Dysfunction, CD274, CAF, Merck18, CD8, and IFNG (Fig. [175]9C–H). Conversely, they exhibited negative correlations with Exclusion, MDSC, MSI.Expr.Sig, and TAM.M2 (Fig. [176]9I-[177]L). Collectively, the intricate interactions between APOE + monocytes and the tumor microenvironment not only significantly influence tumor development, but also profoundly impact immunotherapy efficacy and overall survival. These findings not only enhance our insights into HNSCC pathogenesis, but also hold potential clinical value for refining immunotherapeutic strategies and improving patient outcomes. Furthermore, these results underscore the critical role of further investigating APOE-positive monocytes in tumor biology, opening new avenues for future therapeutic developments. Differential expression of ARGs in HNSCC and its prognostic analysis In this study, we conducted a comprehensive evaluation of the differential expression of ARGs in head and neck squamous cell carcinoma (HNSCC) and investigated its correlation with prognosis. By analyzing the differences in the expression of ARGs in normal and tumor samples derived from single-cell data, we identified 12 differentially expressed ARGs using the R package “limma” (Fig. [178]10A). Further univariate and multivariate Cox regression analyses were performed to reveal significant predictors of HNSCC prognosis. The results indicated that GSTP1 and age are significant independent predictors of HNSCC prognosis, with high expression of GSTP1 being correlated with poorer prognosis (Fig. [179]10B, [180]C). In addition, the KM curves further confirmed that the survival probabilities predicted by the column line graphs coincided with the actual observed survival data (Fig. [181]10D). Based on these findings, we constructed a column chart using the “rms” package, which included variables such as GSTP1 expression, age, and pN stage, allowing us to predict the overall survival of HNSCC patients at 1, 3, and 5 years (Fig. [182]10E). The consistency index (C-index) of the column-line plot was 0.614, and the calibration curve demonstrated a good predictive effect (Fig. [183]10F), indicating that the model had good discriminatory ability. Our findings highlight the potential of genomics in cancer prognostic assessment and provide clinicians with a practical tool to assess patient prognosis in an individualized manner while potentially guiding treatment decisions. Fig. 10. [184]Fig. 10 [185]Open in a new tab Differential expression of ARGs in HNSCC and their prognostic significance. A A volcano plot illustrating the changes in gene expression. B One-factor Cox regression analysis results. C Multifactor Cox regression analysis results. D Comparison of overall survival between patients with high and low GSTP1 gene expression. E Column line plots of survival probability prediction models at 1, 3, and 5 years. F Assessment of the prediction models' accuracy in predicting 1-, 3-, and 5-year survival Discussion ARGs play a role in regulating several aspects of the cellular life cycle, DNA repair, protein homeostasis, cellular metabolism, and immune response. However, the relationship between cellular senescence and tumors is complex [[186]50]. Cellular senescence is a natural tumor suppressor mechanism that limits the proliferation of damaged or mutated cells [[187]51, [188]52]. This mechanism works by activating tumor suppressor pathways such as p53, p16^INK4a, and Rb, which help prevent the development of precancerous cells [[189]53–[190]55]. Despite this protective role, senescent cells accumulate and secrete a variety of inflammatory factors, growth factors, and proteases, a phenomenon known as SASP. SASP can promote the formation of the tumor microenvironment and facilitate the growth, invasion, and metastasis of cancer cells [[191]56]. In addition, SASP may affect immune surveillance functions, making it more difficult for the immune system to clear cancer cells [[192]57]. In conclusion, cellular senescence is a complex biological process that serves as an anticancer mechanism to limit the proliferation of damaged cells, but can also promote tumor development through pro-inflammatory signals and other mechanisms. Through our analysis of cells with high and low levels of ARGs levels, we identified ARGs-related marker genes, which, in combination with nomogram, can predict patient survival more accurately than traditional methods and provide personalized, quantitative prognostic assessment. In the tumor microenvironment of head and neck squamous cell carcinoma (HNSCC), the expression levels of ARGs in various cell types undergo significant changes compared to normal tissues. These alterations may affect the ability of immune cells, such as T cells and monocytes, to combat tumors. For monocytes, the upregulation of HSPA1B and HSPA1A suggests that these cells may be in a highly stressed state in the tumor microenvironment, enabling them to adapt to stressors such as hypoxia and nutrient deficiencies [[193]58]. Additionally, the upregulation of PLCG2 may enhance cellular activation and function, promoting inflammatory responses, cell migration, and phagocytosis [[194]59]. Conversely, the downregulation of APOE may impair the lipid metabolism and anti-inflammatory function of monocytes, thereby enhancing the inflammatory state in the tumor microenvironment and favoring tumor growth [[195]60, [196]61]. Similarly, the downregulation of C1QA reduces the activity of the complement system, which is essential for clearing pathogens and damaged cells, potentially attenuating the immune clearance of tumor cells [[197]62]. In addition, the downregulation of CEBPB may affect the differentiation and function of monocytes, further impairing their ability to regulate inflammatory responses and immune surveillance. The downregulation of LMNA impacts the nuclear structure and stability of monocytes and may indirectly affect their migratory ability and control of tumor growth [[198]63]. In the tumor microenvironment of HNSCC, monocytes with high and low expression of ARGs show significant differences in cellular communication.IL-15 reveals important immunoregulatory changes in the tumor microenvironment through its specific receptor IL-15RA and the IL-2RB and IL-2RG receptors shared with IL-2. Although IL-15 shares some of its receptors with IL-2, it exerts unique functions through IL-15RA, which is particularly important in the survival and activation of NK cells and CD8 + T cells. Upregulation of this pathway may both enhance the anti-tumor response of NK cells and effector T cells and be utilized by tumor cells to promote their survival and immune escape [[199]64]. Furthermore, the downregulation of communication of monocytes with CD8 + T cells, macrophages, endothelial cells, and Tregs by highly expressed ARGs, especially in the chemokine signaling pathway, suggests that the aging process of the tumor microenvironment may weaken immune surveillance and anti-tumor responses. Such changes may be triggered by inflammatory, metabolic changes or immunosuppressive signals released by tumor cells, further promoting immune escape and tumor growth. Thus, upregulation of IL-15 and IL-2 pathways, while having the potential to enhance immune responses, may be equally hijacked by tumor cells in the complex tumor microenvironment, attenuating immune surveillance and promoting tumor progression. Therapeutic strategies targeting HNSCC should consider the dual role of these signals to avoid immunosuppression and improve efficacy. stLearn, a toolkit for spatial transcriptome data analysis, enables spatial localization and developmental trajectory analysis of cells [[200]65]. Clusters of cells with high expression of ARGs were observed in tumor sections, and developmental trajectory analysis revealed these cells infiltrating into surrounding tissues. Following RCTD inverse convolution analysis, which determines the distribution and proportions of different cell types in the tissue, Mistyr can further analyze the interactions among these cell types. In tumor tissues, interactions between monocytes and other cells primarily occur intraview, meaning within the same view or within the same cell type. This suggests that in the tumor microenvironment, monocytes may be regulated more through internal signaling mechanisms rather than direct interactions with other cell types. Additionally, tumor cells may inhibit direct interaction with immune cells through various mechanisms to promote immune escape. The reduced interaction between monocytes and other cells may be one outcome of tumor cells evading surveillance by the host immune system [[201]66]. The significant expression of S100A9/A8 and TLR4 in tumor tissues suggests that these proteins may play important roles in regulating the inflammatory environment of tumors, immune responses, and tumor growth and metastasis [[202]67]. APOE (apolipoprotein E) is a protein primarily expressed in the liver and brain, playing a crucial role in lipid metabolism, cholesterol transport, cellular repair, and immune regulation [[203]68]. An imbalance in lipid metabolism is a key feature of aging and is closely associated with many age-related diseases. Our study found that APOE expression in monocytes was closely related to the prognosis of patients with squamous cell carcinoma of the head and neck. Notably, patients with high APOE expression in monocytes had worse prognoses than those with low expression. Furthermore, high APOE expression also correlated with poorer immunotherapy outcomes. To assess the potential of ARGs as prognostic biomarkers for patients with head and neck squamous cell carcinoma, we utilized Nomogram model, which combines molecular markers with clinical features, provides a personalized, quantifiable prognostic assessment tool that more accurately predicts survival probabilities for individual patients than traditional methods. This model assists physicians and patients in developing more targeted treatment plans. By predicting disease progression and treatment response based on patient-specific genetic and clinical characteristics, this prognostic model helps personalize the treatment of patients with squamous cell carcinoma of the head and neck, providing the most appropriate treatment plans. Conclusion In conclusion, in this study, we employed a combination of single-cell sequencing and spatial transcriptomics to explore the cellular and molecular characteristics of head and neck squamous cell carcinoma (HNSCC) and its relationship with the tumor microenvironment. We observed that HNSCC tissues exhibited a decrease in B cells and an increase in monocytes and macrophages, with monocytes showing significantly enhanced ARGs activity. Pseudotemporal analysis elucidated the time-dependent differentiation process among these immune cells. Spatial transcriptomics further revealed a notable increase in APOE + monocytes within the tumor region, closely correlated with tumor growth and metastasis. These findings not only deepen our understanding of HNSCC pathogenesis but also offer crucial theoretical insights and novel perspectives for developing personalized therapeutic strategies. Limitations This study systematically examined the role of aging-related genes (ARGs) in head and neck squamous cell carcinoma (HNSCC) and offered key insights. Nevertheless, several potential limitations should be noted. First, while the sample size was sufficient to support the preliminary analyses, future studies should focus on increasing the sample size to validate and extend these findings across a broader patient population. This will allow for a more comprehensive understanding of gene expression differences and their potential association with disease progression. In addition, while single-cell and spatial transcriptomics technologies provide cutting-edge methods for resolving cellular heterogeneity within the tumor microenvironment, current sequencing technologies have not yet achieved true single-cell resolution, and imaging-based approaches face challenges in high-throughput gene detection. Finally, although tools such as Seurat, tSNE, and SCTransform demonstrated excellent performance in data processing, their normalization and downscaling steps may introduce analytical bias. Supplementary Information [204]Additional file 1. ^(24.7KB, xlsx) Acknowledgements