Integrating scRNA and bulk-RNA sequencing develops a cell senescence signature for analyzing tumor heterogeneity in clear cell renal cell carcinoma

Introduction Cellular senescence (CS) plays a critical role in cancer development, including clear cell renal cell carcinoma (ccRCC). Traditional RNA sequencing cannot detect precise molecular composition changes within tumors. This study aimed to analyze cellular senescence’s biochemical characteristics in ccRCC using single RNA sequencing (ScRNA-seq) and traditional RNA sequencing (Bulk RNA-seq). Methods Researchers analyzed the biochemical characteristics of cellular senescence in ccRCC using ScRNA-seq and Bulk RNA-seq. They combined these approaches to identify differences between malignant and non-malignant phenotypes in ccRCC across three senescence-related pathways. Genes from these pathways were used to identify molecular subtypes associated with senescence, and a new risk model was constructed. The function of the gene DUSP1 in ccRCC was validated through biological experiments. Results The combined analysis of ScRNA-seq and Bulk RNA-seq revealed significant differences between malignant and non-malignant phenotypes in ccRCC across three senescence-related pathways. Researchers identified genes from these pathways to identify molecular subtypes associated with senescence, constructing a new risk model. Different subgroups showed significant differences in prognosis level, clinical stage and grade, immune infiltration, immunotherapy, and drug sensitivity. Discussion Senescence signature markers are practical biomarkers and predictors of molecular typing in ccRCC. Differences in prognosis level, clinical stage and grade, immune infiltration, immunotherapy, and drug sensitivity between different subgroups indicate that this approach could provide valuable insights into senescence-related treatment options and prognostic assessment for patients with ccRCC. The function of the gene DUSP1 in ccRCC was validated through biological experiments, confirming its feasibility as a novel biomarker for ccRCC. These findings suggest that targeted therapies based on senescence-related mechanisms could be an effective treatment option for ccRCC.


Introduction
Renal cell carcinoma (RCC) is the most typical malignant tumor in the kidney from renal epithelial cells. RCC patients' five-year overall survival (OS) is below 20% (1). Clear cell RCC (ccRCC) is the most common histology subtype and contributes about 70% to all cases of RCC (2). Surgical excisions are the preferred method of ccRCC treatment. However, the outcome is disappointing and is expected to be recurrent (3). The pathogenesis of ccRCC is influenced by the tumor microenvironment (TME), which includes malignant tumor cells, tumor-associated macrophages (TAMs), CD8 T cells, and fibroblasts (4). Even though many biomarkers have been proposed to assess risk models of ccRCC, a significant insufficient value remains with these models (5). Therefore, developing a new predictive model based on the significant genes and pathways of ccRCC is critically essential.
Cell senescence (CS) is a permanent break in the cell cycle and can be caused by various physiological and pathological situations, including tissue injury, aging, and tumorigenesis (6). CS can suppress the unregulated proliferation of cancer cells, thus inhibiting tumor progression (7). Studies also indicated that CS was a robust prognostic biomarker for many cancers in which senescent cells in tumor tissues may facilitate the proliferation and invasions of adjacent pre-neoplastic cells (8). However, the biological mechanisms and prognostic roles of senescence-related genes remain unclear. There is currently little understanding of whether the characteristics of CS within the ccRCC samples can be used for treatment guidance and screening out the prognostic risk subgroups.
Traditional RNA sequencing (bulk RNA-seq) profiling is conducted on a mixed population of cells, which is insufficient for detecting particular cell types and unable to assess the complexity of intra-tumor heterogeneity (9). In contrast, the single-cell RNA-seq (scRNA-seq) technique has flourished recently, allowing researchers to intuitively identify specific gene characteristics at a genome-wide scale and investigate cellular heterogeneity. In this study, we integrated scRNA-seq and bulk RNA-seq data to analyze senescence-related pathways with "ccRCC characteristics" at multiple levels. Based on the genes in these distinct senescencerelated pathways, we constructed senescence-related subtypes and risk models and identified significant differences between different subtypes and risk groups in biology and clinical phenotypes. Additionally, we discovered a novel molecular marker, DUSP1, in ccRCC. These findings highlight the important value of cellular heterogeneity in ccRCC and lay the foundation for further development of clinically relevant applications.

Data acquisition
The scRNA-seq files were collected from GSE159115 via the GEO database. The dataset includes 14 samples. Among them, seven lesions were ccRCC tumor samples, one lesion was a Chromophobe RCC tumor sample, and six lesions were canceradjacent normal tissues. For ccRCC clinical phenotype data, 526 ccRCC samples and 72 cancer-adjacent normal samples were downloaded from The Cancer Genome Atlas (TCGA) dataset that matched survival information (survival time and survival status). TCGA ccRCC gene-level copy number variation (CNV) data of Masked Copy Number Segment type assessed using the GISTIC2 method was performed. The predictive value for the risk model was validated using RECA-EU data of 91 ccRCC samples downloaded from the ICGC Data Portal. Furthermore, we utilized the GSE167573 dataset downloaded from the GEO database and the E-MTAB-1980 dataset downloaded from the ArrayExpress database (https://www.ebi.ac.uk/arrayexpress/) to further validate the performance of our risk model. Finally, to validate the model's application in immunotherapy, we downloaded the datasets GSE78220 and GSE135222, along with their respective clinical information, from the GEO database. Additionally, we obtained the expression matrix and clinical data for IMvigor210 using the R package "IMvigor210CoreBiologies".
Single-cell data processing Data analysis was conducted using Seurat R package (version 3.6.3, https://satijalab.org/seurat/). At first, by setting the criteria that each gene expressed in no fewer than 3 cells and no fewer than 250 genes expressed in each cell, single-cell data were filtered, and 32352 cells were obtained. Then, according to the criteria that each cell expressed 100 to 5000 genes, 25% less mitochondrial content, and 100 to 50000 unique molecular identifier (UMI) counts. The proportions of mitochondria and rRNA were calculated using the PercentageFeatureSet function. The data of the 14 samples were individually normalized using log normalization. Hypervariable genes were screened with the FindVariableFeatures function (variant features were identified based on a variance stabilizing transformation [vst]). Subsequently, we used the canonical correlation analysis (CCA) method to identify the sample batch with the FindIntegrationAnchors function and integrated the 14 samples using the IntegrateData method genes were subjected to scaling using the ScaleData function, and anchors were obtained through principal component analysis (PCA) dimensionality reduction. The condition dim = 25 was set for cell clustering using FindNeighbors and FindClusters functions (Resolution = 0.2). We next performed t-distributed Stochastic Neighbor Embedding (tSNE) for dimension reduction on cells using the RunTSNE function and annotated subpopulations using several classical immune cell markers. Marker genes in subpopulations were screened using the FindAllMarkers function with the logarithm of fold change (log FC) = 0.5, minimal percent of the differentially expressed gene (Minpct) = 0.5, and adjusted p< 0.05 as screening thresholds.

Construction of molecular subtypes and risk model
We conducted consensus clustering on gene expression profiles using the ConsensusClusterPlus package. We simultaneously utilized the PAM algorithm and "Euclidean" metric distance and performed 500 bootstraps, with 80% of patients in the training set for each bootstrap. With the number of clusters set as 2 to 10, the consensus matrix and the consensus cumulative distribution function (CDF) were analyzed to calculate the optimal clusters. LASSO Cox regression was achieved using the R package glmnet. Then, the R package timeROC was utilized to perform the ROC analysis of RiskScore for prognostic classification.
Acquisition and quantification of gene sets 31 genes related to cell cycle progression (CCP) and 24 genes related to angiogenesis were obtained from previous studies (10). In addition, we downloaded 27 genes related to the G1/S phase from the KEGG website. To analyze differences in cell cycle scoring, tumor metastasis scoring, inflammatory response scoring, and telomerase scoring among different subtypes, several common gene sets, including HALLMARK_G2M_CHECKPOINT, HALLMARK_ EPITHELIAL_MESENCHYMAL_TRANSITION, HALLMARK_ INFLAMMATORY_RESPONSE, HALLMARK_HYPOXIA and REACTOME_TELOMERE_EXTENSION_BY_TELOMERASE were downloaded from the GSEA website (https://www.gseamsigdb.org/gsea/msigdb/). Finally, we obtained relevant gene sets for 10 tumor-related pathways previously reported in the study, involving different cancer aspects (11). All gene sets were calculated using ssGSEA analysis to obtain corresponding scores.

SA-b-galactosidase detection assay
The activity of SA-b-gal with 786-O cell was performed in accordance with the manufacturer guidelines (C0602, Beyotime Biotechnology, Shanghai, China). The stained cells were visualized under an inverted microscope.

Statistical analysis
Apart from the stated bioinformatic methods, R (version 4.1.0, www.r-project.org) and GraphPad Prism 8.0 were used to analyze this research. Spearman's rank correlation was utilized to evaluate the connection between two continuous variables. The relevance of two group divergences was tested using Student's t-test. p< 0.05 was statistically significant.

Single-cell clustering and dimension reduction analysis
After filtering, a total of 27,300 cells were obtained. As shown in Figure S1A, the UMI count was significantly related to the number of mRNAs, while the amount of UMI/mRNAs was insignificantly linked to the number of mitochondrial genes. Figures S1B, C represents the violin plot before and after quality control. Next, tSNE was performed on 27300 cells for dimension reduction using the RunTSNE function, and 12 subpopulations were identified and annotated using several classical immune cell markers ( Figure S2). Among these subpopulations, subpopulation 4 was T cell (CD2, CD3D, CD3E, CD3G); subpopulation 10 was B cell (CD79A); subpopulation 11 referred to Mast cell (TPSAB1 and CPA3), subpopulation 1 referred to Macrophage (CD163, CD68, CD14); subpopulations 3 and 9 referred to Fibroblast (ACTA2, PDGFRB, NOTCH3); subpopulations 0, 5, 6, and 8 referred to ccRCC (CA9); subpopulations 2 and 7 referred to Endothelial cell (KDR, PECAM1, PLVAP, PTPRB, VCAM1).
The distribution of the 14 samples was summarized in a tSNE plot ( Figure 1A); different subpopulations after clustering were presented in a tSNE plot ( Figure 1B); the distribution of the annotated cells was visualized in a tSNE plot ( Figure 1C). The numbers of different sample cells before and after filtering are statistically summarized in Table 1. Marker genes in 7 subpopulations were selected using the FindAllMarkers function with logFC = 0.5 and minpct = 0.5. After screening with a corrected p-value< 0.05, we only presented the expression patterns of the top 5 marker genes exhibiting the most significant contribution in the subpopulations ( Figure 1D). The results of the marker genes are described in Table scRNA_marker_gene.txt. Furthermore, we analyzed the percentages of these 7 subpopulations in each sample ( Figure 1E). Next, the CopyKat package was employed to predict CNV changes in cells in single-cell data to distinguish tumor cells from normal cells in each sample (although tumor and normal tissues were selected at the time of sampling, it cannot be guaranteed that tumor tissues did not contain normal cells). Among them were 1944 cancer cells and 25356 normal cells ( Figure 1F).

Cell senescence characteristics in singlecell level
As shown in Figure 2, fibroblast senescence-related pathways had higher scores in malignant cells than in non-malignant cells, and endothelial cells were only present in non-malignant cells.

Validation of cell senescence abnormalities based on bulk RNA-seq data
Our results demonstrate that malignant cells exhibit higher expression of cell senescence-associated pathways than nonmalignant cells at the single-cell level. We examined expression profiles in tumor and normal tissue samples using bulk RNA-seq data to further analyze these pathways. Through GSEA software analysis, we found that GOBP_REGULATION_OF_CELL_AGING, GOBP_NEGATIVE_REGULATION_OF_CELL_AGING, and KEGG_P53_SIGNALING_PATHWAY were significantly enriched in tumor tissues from the TCGA dataset ( Figure 3A).
Next, the scores of these senescence-associated pathways in tumor tissues and normal tissues in each sample in the TCGA dataset were subjected to a ssGSEA analysis. The significance of each cell senescence-associated pathway in tumor tissues and adjacent tumor tissues was calculated, and we found the enrichment scores of GOBP_REGULATION_OF_CELL_AGING, GOBP_NEGATIVE_REGULATION_OF_CELL_AGING, KEGG_P53_SIGNALING_PATHWAY was higher in tumor tissues than in tumor-adjacent tissues ( Figure 3B). The scores of the pathways mentioned above in the TCGA dataset are presented in Table tcga.cellage.score.txt.

Construction of cell senescencerelated subtypes
As described above, our analysis revealed three cell senescencerelated pathways, including (GOBP_REGULATION_OF_CELL_ AGING, GOBP_NEGATIVE_REGULATION_OF_CELL_AGING, and KEGG_P53_SIGNALING_PATHWAY) were significantly enriched in tumor tissues. Therefore, a univariate Cox analysis was implemented on the genes in these three pathways in the TCGA and ICGC datasets using the survival package and screened prognosisassociated genes with p< 0.05. The results showed 57 prognosisassociated genes in the TCGA dataset and 18 prognosis-associated genes in the ICGC dataset. Furthermore, we found that 7 of these genes were correlated with prognosis in both datasets. The results of univariate Cox analysis on the genes in three senescence-related pathways in TCGA and ICGC datasets are outlined in Table  tcga.cellage.cox.txt and icgc.cellage.cox.txt.
Next, based on seven prognostic genes, we performed clustering using the ConsensusClusterPlus package on 526 ccRCC samples from the TCGA dataset. CDF Delta area curves show that the clustering with Cluster = 3 was more stable (Figures 4A, B). Ultimately, we selected k = 3 and obtained three subtypes (cluster) ( Figure 4C). A subsequent analysis on the prognostic characteristics of these three subtypes was conducted, the results of which revealed that these three subtypes had notable differences in prognoses ( Figure 4D). Overall, the clust1 subtype exhibited the best prognosis, followed by the clust2 subtype (second) and clust3 subtype (worst).
Additionally, we used the same method to analyze the patients in the ICGC dataset and observed marked differences in prognoses among these three molecular subtypes ( Figure 4E), which coincided with the results of the TCGA dataset. The results above suggested the

Differential analysis of clinical phenotypes in cell senescence-related subtypes
For the patients in the TCGA dataset, the distribution of various clinical features in the three molecular subtypes was compared. The corresponding results demonstrated significant differences among the three subtypes concerning gender, T stage, N stage, Stage, Grade, and survival status of patients in the TCGA dataset ( Figure 5). Additionally, we performed a comparative analysis of various clinical features among the three molecular subtypes in 91 ccRCC patients from the RECA-EU dataset. The results of our analysis demonstrated significant differences in the T stage and nearly significant differences in the M stage among the three subtypes. However, no significant differences were observed in the N stage grouping. (Supplementary Table 1).

Differences in variations of cell senescence-related subtypes
We next integrated copy-number variants (CNVs) in TCGA-KIRC using the gistic2 software under a confidence level of 0.9, with hg38 as the reference genome. As presented in Figures 6A-C, differences were noted in CNVs among the three subtypes. Also, the maftools package was employed to analyze the single nucleotide variant (SNV) data downloaded from TCGA, from which the top 15 genes with the most variations were selected and visualized ( Figure 6D).

Biological characteristics of cell senescence-related subtypes
The CCP score in each sample from the TCGA dataset was calculated using the ssGSEA method. The results indicated that the clust3 subtype, which had the worst prognosis, exhibited a higher CCP score. ( Figure 7A). Previous research has demonstrated that tumor cells can suppress the induction of cell senescence in the cell cycle, and an essential characteristic of senescent cells is that upregulation of cyclin-dependent kinases such as INK4a and p21 can lead to cell cycle arrest (15).. It is noteworthy that the results also demonstrated an increase in G1/S phase-and G2 checkpointrelated scores in the clust3 subtype ( Figures 7B, C). These data indirectly illustrated that the cell cycle was not the only influencing factor for cell senescence, and other mechanisms in the body may act together with the cell cycle to regulate cellular senescence.
Likewise, telomerase inhibition induces cellular senescence (15). Tumor cells often activate the telomerase activity to prevent the loss of telomeres in the body. Telomere Extension by Telomerase was the primary function of this pathway. The results of this analysis indicate that the clust3 subtype with a worse prognosis had a higher score of telomere extension by telomerase ( Figure 7D). Nevertheless, apart  from transmitting a "please kill me" message, factors secreted by senescent cells can impact adjacent cells, thus hastening tumor migration and metastasis by inducing epithelial-mesenchymal transition (EMT). Meanwhile, senescent tumor cells can expedite the formation of blood vessels and lymphatic vessels by recruiting specific macrophages and also supply oxygen and nutrients for the growth of other tumor cells, thereby facilitating tumor growth and metastasis. The results showed a higher EMT score of cluster 3 ( Figure 7E). Furthermore, clustering analysis revealed that angiogenesis and hypoxia scores were significantly lower in the second cluster of patients ( Figures 7F, G). Next, the immune score and stromal score of samples from TCGA were estimated using the ESTIMATE method. According to these two scores, a higher degree of immune infiltration was noted in the clust3 subtype ( Figure 7H). The CIBERSORT method scored 22 immune cells from the TCGA dataset. It was found that cell senescence-associated subtypes significantly differed among some immune cells. Macrophages, in particular, exhibit a notable variance in their level of infiltration across distinct senescence subgroups ( Figure 7I). We scored the enrichment of 10 tumor-related pathways in samples from TCGA-KIRC and observed significant differences in all 10 pathways ( Figure 7J). Notably, the inflammation score of patients in the clust3 subtype was significantly higher than that of patients in the clust1 and clust2 subtypes, as shown in the clustering analysis results ( Figure 7K). These findings suggest a potential correlation between the molecular subtype of ccRCC and the TME, which may have implications for personalized treatment strategies in the future.
Construction and validation of cell senescence-related risk model As described above, we identified three different molecular subtypes through 7 essential genes and found differences in clinical phenotype, mutation, and immune characteristics among subtypes. The clust3 subtype showed the worst prognosis, followed by the clust2 subtype, while the clust1 subtype exhibited the optimal prognosis. Then, we conducted differential analyses of the clust1 vs. no_clust1 subtype, clust2 vs. no_clust2 subtype, and clust3 vs. no_clust3 subtype with the limma package, the results of which are summarized in   Comparison of different clinical features distribution among three molecular subtypes in TCGA dataset.  (Figures 8A-C). Next, a univariate Cox analysis of the 964 DEGs was realized using the coxph function of the survival package and identified 613 genes showing significant effects on prognosis (p< 0.05), consisting of 100 "risk" genes and 513 "protective" genes ( Figure 8D). The corresponding results are summarized in Table tcga.cox.txt. Furthermore, a LASSO regression analysis on the 613 essential genes was implemented to reduce the number of genes for the risk model. The trajectory of each independent variable was analyzed. It was suggested that number of independent variables with a coefficient close to 0 increased gradually as the lambda increased gradually ( Figure 8E). We employed 10-fold cross validation for model construction and analyzed the confidence interval for each lambda ( Figure 8F). It should be noted that the optimization model was developed using lambda = 0.0386. Therefore, we selected 21 genes at lambda = 0.0386 as target genes for subsequent analysis. Following stepwise regression, the number of genes was reduced from 21 to 10. Ultimately, a 10-gene signature including ITGA8, SEMA3G, DPYSL3, IFITM1, ZNF521, SOCS3, PCSK6, DUSP1, SLC44A4, and IL20RB was developed, and Senescore was calculated using the formula:

Prognostic analysis and validation of the risk model
Next, we calculated the risk score for each sample individually based on the expression profiles of samples in the training dataset from TCGA data. We separately analyzed the classification efficiency of this model to predict 1-year to 5-year prognoses. As depicted in Figure 9A, this model had a high area under the curve (AUC) value; At last, we performed Z-score normalization on the Riskscore and assigned the samples to a high-risk subgroup (Riskscore > 0 after Z-score normalization) and low-risk subgroup (Riskscore< 0 after Z-score normalization). Kaplan-Meier (KM) curves were plotted accordingly. The analysis indicated significant differences between the two groups in terms of overall survival, disease-specific survival, and progression-free survival in TCGA-KIRC cohort ( Figure 9B, Supplementary  Figures 2E-G). To validate the robustness of the model, the same method was applied for validation using the ICGC dataset. The AUC values of the risk model established with the 10 genes mentioned above are presented in Figure 9C. After the Z-score normalization of Riskscore. The samples were assigned to the highrisk subgroup and low-risk subgroup. The KM curves showed significant differences between these two groups in the ICGC dataset (p< 0.05). Furthermore, we validated the robustness of our model in two additional microarray datasets, as shown in Supplementary Figures 2A-D

Correlations between risk model and clinical characteristics
To ascertain the correlations between the RiskScore and the clinical characteristics of ccRCC, we analyzed the differences in the RiskScore among different TNM grades and Stages in the TCGA-KIRC dataset. The results exhibited that a higher clinical grade was associated with a higher risk score ( Figure 10).

Biological characteristics of cell senescence-related risk score
Based on the results mentioned above, cell senescenceassociated subtypes were associated with cell cycle, telomere extension by telomerase, hypoxia, angiogenesis, and immunity. Therefore, the correlations of those scores with the cell senescence-related risk score were analyzed using the rcorr function in the Hmisc package ( Figure 11).

Prediction of immunotherapeutic efficacy with the cell senescence-related risk score
To assess the relevance of the Senescore to immunotherapy, we evaluated the predictive capability of the Senescore for patient response to ICB therapy. In the anti-PD-L1 cohort (IMvigor210 cohort), a high Senescore was associated with a worse prognosis ( Figure 12C; log-rank test, p< 0.05). Additionally, we found the FIGURE 10 Senescore differences in clinical pathological features in TCGA dataset. ns, p≥0.05; *, p< 0.05; **, p<0.01; ***, p<0.001; ****, p<0.0001.

FIGURE 11
Association between the Senescore and biological feature scores. varied response of 348 patients in the IMvigor210 cohort to PD-L1 blockers, encompassing complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD). The patients with SD/PD had a higher Senescore versus those with CR/PR ( Figure 12A). Based on percentage statistics between the low-and high-Senescore groups, significantly better treatment outcomes were noted in patients with a low Senescore ( Figure 12B). We analyzed survival differences among all samples in the IMvigor210 cohort as well as the survival differences at different Stages. The results showed significant survival differences among Stage I + II samples ( Figure 12D), but insignificant survival differences between low-and high-Senescore groups in Stage III + IV samples ( Figure 12E). Particularly, the Senescore exhibited outstanding predictive performance in early-stage clinical samples.
The result from TIDE showed that a higher TIDE prediction score denoted a higher probability of immune escape, indicating a lower probability that the patient might benefit from immunotherapy. Furthermore, the Senescore and TIDE scores were higher in non-responders to immunotherapy, which indirectly suggested that patients with high Senescore were less prone to benefit from immunotherapy ( Figures 12F, G).
Additionally, we assessed the efficacy of Senescore on other immunotherapy cohorts. Our findings indicate that the GSE135222 cohort shows a different response when compared to the IMvigor210 cohort. Patients with high Senescore performances exhibited enhanced benefits from immunotherapy and sustained better prognoses. Nevertheless, in the GSE78220 immunotherapy cohort, the high and low Senescore divisions did not function as dependable predictive markers for patient outcomes (Supplementary Figure 3).

Correlations of cell senescence-related risk scores and drug sensitivity
This study also compared the response of high-and low-risk populations to conventional chemotherapeutic agents such as Erlotinib, MG-132, and Paclitaxel. The high Senescore group showed a higher sensitivity to the abovementioned agents ( Figure 13).

Senescore integrated with clinicopathological features for improved prognostic models and survival prediction
Through univariate and multivariate Cox regression analyses of the Senescore and clinicopathological characteristics, the Senescore was identified as the most significant prognostic factor (Figures 14A, B). To quantify patients' risk assessment and survival probability, we integrated the Senescore and other clinicopathological characteristics to establish a nomogram ( Figure 14C). Based on the results. The Senescore exhibited the most significant impact on survival prediction. Furthermore, a calibration curve was utilized for the predictive accuracy assessment. As presented in Figure 14D, the calibration curves for the prediction at the three calibration points (1, 3, and 5 years) almost coincided with the standard curves, suggesting the good predictive performance of the nomogram. Also, decision curve analysis (DCA) was carried out to evaluate the reliability of this model. The benefit of either Senescore or nomogram was remarkably higher than that of the extreme curve. Compared to other clinicopathological characteristics, nomogram, and Senescore exhibited more vital survival prediction ability ( Figure 14E).

Discussion
Based on single-cell data analysis, this study interpreted the abnormalities in cell senescence-associated pathways within the TME. Meanwhile, we screened cell senescence-associated pathways enriched in tumor tissues in bulk RNA-seq dataset by GSEA and constructed cell senescence-associated subtypes and risk models based on the genes in these pathways. Moreover, this study linked the efficacy of immunotherapy to cell senescence-associated risk scores. ccRCC is a highly heterogeneous renal tumor developed by different complex epigenetic driving mechanisms and molecular pathways. Due to its malignant progression and high recurrence rate, ccRCC is the deadliest type of renal tumor. Currently, postoperative recovery in patients with ccRCC remains unsatisfactory, and only a few patients may benefit from drug therapy targeting tyrosine kinase inhibitors (TKI) and anti-PD-1 antibodies (16). Biomarkers that can accurately predict prognosis and guide the treatment of ccRCC have not been fully identified and applied clinically. Immunosenescence refers to the age-related decline in immune system function, which can impair the body's ability to defend against infections, vaccines, and tumors (17). As a result, older individuals are more susceptible to infections and less capable of mounting an effective immune response to diseases. The increased risk of cancer in older adults is of particular concern, as tumor cells often exploit deficiencies in the immune system to evade detection and clearance, making treatment more challenging (18). Consequently, immunosenescence has emerged as a prominent topic in medical research. By exploring the underlying mechanisms of immunosenescence and its impact on human health, researchers hope to develop new strategies to enhance immune function and improve the health and well-being of older individuals. Cellular senescence is a permanent state of the arrest of the cell cycle, and senescent cells/genes have been observed to accumulate during aging and play a role in various tumors (19). Xu et al. constructed a senescence-associated prognostic model significantly correlated with lung adenocarcinoma diagnosis and prognosis (20). In addition, a previous study showed that senescence-associated protein P400 is a prognostic marker for renal cell carcinoma (21). However, senescence characteristics and senescenceassociated prognostic models of ccRCC are still rare. The specific molecular markers for ccRCC also need to be further elucidated.
Many studies have attempted to build models based on gene sequencing and clinical data to predict the prognosis of patients with ccRCC (22,23). However, the clinical application of these models has had little effect. In this study, we used RNA-seq combined with bulk RNA-Seq to analyze the cellular senescence characteristics of ccRCC. Our assessment has recognized IFITM1, SOCS3, DPYSL3, IL20RB, SLC44A4, SEMA3G, ITGA8, PCSK6, ZNF521, and DUSP1 as genes that exhibit noteworthy differential expression in association with senescence in patients with ccRCC. A prognostic model based on these 10 genes is established. Our senescence-related prognostic features showed good performance, with the AUC values of 1-year, 3-year, and 5-year OS predicted by the TCGA database being 0.84, 0.81, and 0.79, respectively. We also used the ICGC dataset to further validate and evaluate the prognosis. The ICGC database predicted the AUC value of 1-year, 3-year, and 5-year OS to be 0.6, 0.69, and 0.7, respectively. The results showed that the senescence-associated prognostic model we constructed could predict the survival rate of ccRCC. Regrettably, extending this gene signature to other types of Differential analysis of Senescore and IC50 drug sensitivity. ****, p<0.0001.
Gong et al. 10.3389/fimmu.2023.1199002 tumors is unfeasible. Currently, the markers of aging in tumors remain unclear and may vary depending on tumor type and subtype. Each tumor type has distinct biological and genetic characteristics, which could affect the expression of senescence markers in tumor cells. Therefore, a more thorough analysis and evaluation must first be conducted before investigating the applicability of aging markers in specific tumors. In the future, we hope to identify more generalized, broad-spectrum aging markers that can be utilized across different tumor types, facilitating more precise and convenient guidance for studying tumor subtypes. DUSP1 is a threonine-tyrosine bispecific phosphatase that targets negatively regulating the MAPK signaling pathway (24). Nevertheless, its role in tumorigenesis is controversial. DUSP1 is highly expressed in a range of tumors, including lung, breast, ovarian, gastric, and prostate cancers (25-29), while low expression in HCC (30). Several studies have demonstrated the involvement of DUSP1 in malignant tumor progression through various signaling pathways (31). However, the role of DUSP1 in ccRCC remains unknown. The results of this study indicate that DUSP1 is one of the essential regulatory genes of senescence characteristics in ccRCC. The results of specific staining experiments demonstrated that DUSP1 plays an essential role in the senescence of renal clear cell carcinoma cells by effectively inhibiting the generation of senescent cells. In vitro experiments, we found that DUSP1 knockdown significantly promoted the proliferation, migration, and invasion of renal carcinoma cells. Targeting DUSP1 may serve as a potential "senescence biomarker" for predicting clinical outcomes in patients with ccRCC.
The tumor microenvironment (TME) comprises various types of cells, including tumor cells, inflammatory cells, immune cells, stromal stem cells, endothelial cells, and tumor-associated fibroblasts, which play essential roles in the proliferation and drug resistance of tumor cells. Cellular senescence in TME can trigger immune cell infiltration and promote tumor progression (32,33). This study found that T cells, B cells, mast cells, macrophage, fibroblast, and endothelial cells were significant tumor-infiltrating cell clusters compared to adjacent normal tissues. We also found that fibroblast senescence-related pathway scores were higher in ccRCC than adjacent normal tissue cells. Tumor fibroblasts are the primary source of tumor extracellular matrix (ECM) (34). Shi et al. demonstrated that the recruitment of cancer-associated fibroblasts (CAFs) in the ccRCC microenvironment occurs through interaction with malignant proximal renal tubular epithelial cells (PTEC) (35). Peng et al. showed that infiltrating CAFs could reduce CD8 + T cell infiltration in the TME of ccRCC by secreting galactolectin 73 (Gal1) (36). The study found that CAFs can provide metabolic support to cancer cells by releasing alanine and deoxycytidine, thereby enhancing chemotherapy drug resistance (37,38). Extracellular vesicles from CAFs have also been shown to contain various surface proteins (CD105, Hsp70, TGF-b1, etc.) and metabolites (lactate, amino acids, lipids, etc.), which can affect tumor progression and drug resistance (39). CAFs may represent a novel therapeutic target to combat resistance to ccRCC treatment. Targeting CAFs with immunotherapy is also emerging as an effective treatment for ccRCC.
This study investigated the correlation between multiple molecular characteristics, such as telomerase, EMT, and angiogenesis, representing different physiological processes and cell senescence scores. Furthermore, the study explored the relationship between immune infiltration levels in the tumor microenvironment and cell aging scores. Although specific impacts of each characteristic on immune infiltration have been reported, for instance, the telomerase catalytic subunit (TERT) activates endogenous retroviruses to promote the formation of a tumor immune suppression microenvironment (40). At the same time, epithelial-mesenchymal transition (EMT) plays a vital role in immune evasion and tumor immune suppression. The EMT score may be a predictive biomarker for clinical response to immune checkpoint blockade. In addition, tumor angiogenesis is associated with the infiltration of different immune cell types in the TME, which could affect the response of tumors to immunotherapy. These results suggest that cell aging affects immune infiltration through multiple molecular mechanisms. However, further research is necessary to understand the precise impact of cell aging on immune infiltration. Through correlational analysis of these results, it is expected to link cell senescence, immune infiltration, and molecular features, providing a comprehensive understanding and assistance for ccRCC treatment.
Furthermore, this study explored the potential use of risk scores associated with cellular senescence in immune therapy by linking them to its effectiveness. We conducted an effective analysis of the correlation and difference between the senescence score and the infiltration level of immune cells, which is critical when predicting the results of immune therapy using risk scores associated with cellular senescence. Notably, both significant differences in infiltration levels of macrophages among different subtypes and a significant correlation between senescence score and macrophages were shown by the results. M0 macrophages are undifferentiated precursor cells that differentiate into either M1 (pro-inflammatory) or M2 (anti-inflammatory) macrophages according to the needs of the TME. Although we cannot accurately evaluate the precise effects of cell senescence signatures on immune cells, we hypothesize that genes in the model may facilitate the polarization of M0 macrophages toward tumor-promoting M2 macrophages during polarization. For example, SOCS3 plays an important role in the polarization of M2 macrophages (10), but the roles of other genes in macrophage polarization remain unclear. Our main goal is to guide immune therapy for patients by reflecting the infiltration of immune cells in complex TMEs and highly heterogeneous backgrounds using senescence scores. Overall, these results suggest that it is possible to evaluate the infiltration of immune cells in the TME or the molecular characteristics of different tumor patients by using senescence scores. Senescence scores can be used as a useful marker to evaluate the infiltration level of immune cells effectively and conveniently in TME.
Although current research provides a relatively clear understanding of the senescence features in ccRCC and enhances our understanding of the role of cellular senescence in ccRCC, some limitations still cannot be ignored. Although this study validated the results using multiple ccRCC patient cohorts, the limited sample size may not fully represent the characteristics of the entire ccRCC patient population. Additionally, the data used in this study may come from limited datasets or databases and, therefore, may not fully cover all relevant information. Finally, this study did not consider potential external factors such as environmental factors and lifestyle, which may impact the conclusions drawn in this study. In summary, this study uncovered abnormal senescence-related pathways in the TME of ccRCC using a combination of single-cell data and bulk RNA-seq. Based on these dysregulated aging signaling pathways, senescencerelated subtypes and risk models were created, which offer new methods to evaluate prognosis, guide clinical drug selection, and assess immune therapy response in ccRCC patients. Targeting essential senescence-related genes may lead to a new understanding of the molecular mechanisms of aging in ccRCC and their clinical applications.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found within the article/Supplementary Materials.

Ethics statement
All the procedures were performed in compliance with the National Institutes of Health Guidelines for the Care and Use of Laboratory Animals and ARRIVE Guidelines.