Construction and experimental validation of a macrophage cell senescence-related gene signature to evaluate the prognosis, immunotherapeutic sensitivity, and chemotherapy response in bladder cancer

Tumor-associated macrophages (TAMs) are pivotal components of tumor microenvironment (TME), and senescent TAMs contribute to the alternation of the profiles of TME. However, the potential biological mechanisms and the prognosis value of senescent macrophages are largely unknown, especially in bladder cancer (BLCA). Based on the single-cell RNA sequencing of a primary BLCA sample, 23 macrophage-related genes were identified. Genomic difference analysis, LASSO, and Cox regression were used to develop the risk model. TCGA-BLCA cohort (n = 406) was utilized as the training cohort, and then, three independent cohorts (n = 90, n = 221, n = 165) from Gene Expression Omnibus, clinical samples from the local hospital (n = 27), and in vitro cell experiments were used for external validation. Aldo-keto reductase family 1 member B (AKR1B1), inhibitor of DNA binding 1 (ID1), and transforming growth factor beta 1 (TGFB1I1) were determined and included in the predictive model. The model serves as a promising tool to evaluate the prognosis in BLCA (pooled hazard ratio = 2.51, 95% confidence interval = [1.43; 4.39]). The model was also effective for the prediction of immunotherapeutic sensitivity and chemotherapy treatment outcomes, which were further confirmed by IMvigor210 cohort (P < 0.01) and GDSC dataset, respectively. Twenty-seven BLCA samples from the local hospital proved that the risk model was associated with the malignant degree (P < 0.05). At last, the human macrophage THP-1 and U937 cells were treated with H2O2 to mimic the senescent process in macrophage, and the expressions of these molecules in the model were detected (all P < 0.05). Overall, a macrophage cell senescence-related gene signature was constructed to predict the prognosis, immunotherapeutic response, and chemotherapy sensitivity in BLCA, which provides novel insights to uncover the underlying mechanisms of macrophage senescence.


Introduction
Bladder cancer (BLCA) is the sixth most common cancer worldwide, with a median age of diagnosis at 73 years old (Flaig et al. 2022). Despite the development of many novel and effective therapeutic regimens, recent years still have seen 573,000 new cases of BLCA and 213,000 deaths caused by it (Sung et al. 2021). Immunotherapy, such as immune checkpoint inhibitors (ICIs), cisplatin-based neoadjuvant chemotherapy, and targeted therapy, has improved the clinical outcomes of numerous bladder cancer (BLCA) patients. These treatment modalities are currently recommended as first-line or secondline treatments for BLCA, based on the latest clinical guidelines (Rouprêt et al. 2023). Unfortunately, many patients do not respond positively to the current therapeutic prescriptions, particularly with regard to immunotherapy and chemotherapy. Therefore, there is an urgent need for more accurate prognostic indicators of immunotherapeutic sensitivity and chemotherapeutic response, which can inform clinical decision-making. Qijun Jiang and Junhao Zhou contributed equally to this work.

Data source and processing
The 1044 macrophage-related cell markers were extracted from GSE145137 within BLCA tumor, whereas 279 senescence-related genes were retrieved from the Rui Sun et al. (Sun et al. 2022) and the CellAge public database (http:// genom ics. senes cence. info/ cells) (Supplementary Table 1). The gene expression data with fragments per kilobase per million mapped reads (FPKM) format and clinicopathological features of 408 samples were accessed from the TCGA-BLCA cohort (Supplementary Table 2). The external validation cohorts, which included GSE31684, GSE32894 and GSE13507, were obtained on (https:// www. ncbi. nlm. nih. gov/ geo/) database. The RNA expression data with the FPKM format and therapy information of the IMvigor210 cohort were extracted from the IMvigor210CoreBiologies package in R software (version 4.4.1). The genomic difference inclusion criteria were |logFC| > 1 and adjusted P < 0.05 calculated by EdgeR package of R (Song et al. 2021). The minimum overall survival of clinical data required 30 days. To discover the developmental trajectories of macrophage and cell-cell communications with other tumor cells, we integrated two bladder cancer SC-RNA sequencing datasets, GSE145137 and GSE135337.

Clinical sample collection
The BLCA tissues from 27 BLCA patients were collected from the local hospital to further verify our findings. All tumor specimens were collected in the center of the tumor and well preserved during the surgery, and all the samples were divided into two parts: one was embedded with paraffin for immunohistochemistry, and other one was stored in liquid nitrogen for RNA extraction. The study received approval from the Ethics Committee of the Third Affiliated Hospital of Southern Medical University (ID: 2023ER024), and all patients provided their signed informed consent files. Among these samples, 12 samples were at TNM stage I-II, while 15 patients were at TNM stage III-IV, where the TNM stages were classified by the eighth TNM staging system defined by the American Joint Commission on Cancer.

NMF cluster and validation
The macrophage cell senescence subtypes of BCa were identified through consensus non-negative matrix factorization (NMF) clustering by using the R package, with the "brunet" standard and 10 iterations. The optimal k value was determined based on cophenetic dispersion and silhouette. The range of K value was 2-10 and the minimum member of each subtype was 10 (Jiang et al. 2021). The principal component analysis (PCA) algorithm verified the clustering results of each signature, while its OS was evaluated by Kaplan-Meier survival curve for each cluster, divided by NMF.

Expression of differential genes and construction of the risk model
Lasso and univariate Cox regression were inducted to identify the number of genes in our model associated with OS of each sample through the glmnet and the survival packages. In the univariate Cox analysis, P value < 0.01 of survival analysis was used to determine statistically significant variables. Subsequently, the multivariate Cox regression analysis compressed the variables further and prevented our data from overfitting. The results determined by Lasso and univariate Cox regressions analysis were included in the multivariate Cox regression analysis with stepwise by the survminer package, with risk model was ultimately constructed. The evaluation index of the risk model was defined based on the macrophage cell senescence-related scores (MSR). This model was constructed by significant markers that were calculated by lasso regression. Coef was their influence value in Fig. 1 The work flowsheet of the present study the model, and gene expression represents gene expression collected from training set (has been converted to TPM format). MSR score was computed using the following formula:

Validation of the biomarkers risk model
Three external validation cohorts: GSE13507, GSE32894, andGSE13507 was employed to verify the robustness of the risk model. The Kaplan-Meier survival analysis with logrank test was performed to detect the prognosis difference between high-and low-ILMRS samples defined by training set through survival package, and P < 0.05 was regarded to be significant. The forest map was drawn to comprehensively evaluate four K-M survival curves, with hazard ratio and its 95% confidence interval, through the desktop software Review Manager (version 5.4). Meanwhile, the risk score distribution was also plotted by the time-dependent receiver operating curves (ROCs) for 1-, 3-, and 5-year OS rates of MSR via the survivalROC package.

Functional enrichment analysis and scRNA-seq data analysis
The GO functional annotation composed by Biological Process, Cellular Component, and Molecular Function 3 categories were conducted through the clusterProfiler R package with P < 0.05 and FDR < 0.05 filtering. To screen macrophage signature markers with BLCA, the expression matrix of 2058 barcodes isolated with primary BCa patients was downloaded from GEO database. The inclusion barcodes using the following criteria: (1) more than 200 genes and less than 7800, (2) less than 3% of mitochondria gene proportion, and (3) less than 1% of Erythrocyte gene proportion. 12 clusters were then classified by UMAP and t-SNE methods subject to top 2000 variable feature genes across all cell samples. The cell categories annotation according to CellMarker website (http:// biocc. hrbmu. edu. cn/ CellM arker/) and retrieved through Singler R package, whereas Coef i * (gene expression) i pseudotime analysis was performed to the monocle package of R.

Immunotherapeutic response and immune cell infiltration
The immune cell infiltration states computed with the datasets (QUANTISEQ, MCPcounter, XCELL, TIMER, EPIC, CIBERSORTx, and CIBERSORT) and the infiltration estimated with R package immunedeconv (https:// grst. github. io/ immun edeco nv). The content of immune cell infiltration between the high-and low-MSR groups was measured by the Wilcoxon signed-rank test. We employed IMvigor210 cohort including 348 malignant metastatic urothelial cancer samples who were treated with atezolizumab and immune checkpoint blocking (ICB) therapy to validate the predictive ability of MSR via the K-M survival curve. After excluding the samples missing immunotherapeutic response, the respective boxplots were drawn by the Wilcoxon signed-rank test, with P < 0.05 filtering to verify whether the model was related to the tumor immunotherapeutic response.

Tumor mutation and cell-cell communication analysis
Utilized by the maftool package in R software, tumor gene mutation frequency of each cluster characterized by NMF from TCGA was analyzed and visualized. The top 15 genes of mutation were plotted at heatmap. Receptor-ligand pairs interaction with cellular communication, after divided by MSR score into high group and low group, were conducted with the cellchat LR database: CellChatDB (Jin et al. 2021). Calculating the communication presumptions of corresponding LR pair interactions, the signaling pathway within intercellular communication was assessed.

Immunohistochemistry
The tumor tissues from the patients with BCa were fixed with 4% paraformaldehyde for > 24 h and then dehydrated through an alcohol gradient method. Tissues were paraffin embedded, sectioned, dewaxed, and unmasked in boiling 10 mM citrate buffer pH 6.0. Sections were washed three times in PBS, blocked in 5% normal goat serum for 1 h, and incubated with primary antibody: rabbit polyclonal anti-TGFB1I1 (proteintech-12869-1-AP), anti-ID1 (Affinity Biosciences Cat# DF2932), and rabbit polyclonal anti-AKR1B1 (proteintech-10624-1-AP) overnight at 4 °C. Both of them were at dilution of 1:500. Sections were then washed and incubated for 30 min with biotinylated secondary antibody (Goat Anti-Rabbit IgG(H+L) antibody, Biyot-orb688927). Nuclei were counterstained with hematoxylin. The sections were washed twice in acidified water and then mounted Fig. 2 Investigation of three BLCA-MSR subtypes. A The consensus non-negative matrix factorization (NMF) clustering. B Different factors distributions at rank = 2-10, the cophenetic and dispersion reflect the stability of the cluster established from NMF, and the silhouette coefficient refers to the description of method. C K-M survival curve for OS among the NMF subgroups. D-F Macrophage infiltration difference calculated by XCELL, MCPCOUNTER, CIBERSORTx, and EPIC algorithms, results showed by boxplot (C1:47, C2:109, C3:243; ns no significance, *P < 0.05, **P < 0.01, ***P < 0.001). G TIDE method performed immunotherapeutic efficacy intergroup ◂ for microscopy. Five visual fields were randomly selected for each section, and their integrated optical density (IOD) values were evaluated by the software Image-Pro Plus 6.0, which was filtered by Welch's two-sample T-test, P < 0.05.

Cell culture and cell viability
The THP-1 and U937macrophage cell line were purchased from Shanghai Institutes for Biological Sciences (Shanghai, China). THP-1 cells maintained between 5 × 105 and 7 × 105 cells/mL were cultivated in RPMI 1640 medium supplemented with 10% FBS and 100 U/mL each of streptomycin and penicillin. U937 cells were cultured in RPMI 1640 medium with 10% FBS and 1% pen/strep, depriving FBS to 1% 12 h before treatments. Cell viability was discovered using CCK8 assay. The cell lines were starved in 2% FBS with 96-well plates. After H 2 O 2 treatment, 20 μL (5.0 mg/ mL) CCK8 solution was added to plates and incubated for 2 h at 37 °C in a 5% CO 2 atmosphere, while the absorbance was measured at 450 nm by enzyme-labeling instrument.

Drug treatment and real-time PCR
The macrophage cell lines were cultured in 24-well plates. After incubation with H 2 O 2 (170 μg/mL) for 48 h, the oxidative-senescence of THP-1 and U397 was established. Then, the medium in each well was changed every 36 h and washed away with PBS (He et al. 2019). The total RNA of the samples was collected employing Trizol (ThermoFisher Scientific, Germany). Subsequently, PrimeScript RT reagent kit (TaKaRa, Osaka, Japan) and SYBR Premix ExTaq kit (Takara, China) were used to synthesize and amplify the cDNA. The CFX 192 Connect Real-Time PCR system (Bio-Rad, U.S.A.) was utilized to perform qRT-PCR procedure, and the mRNA expression of target genes was normalized with the levels of GAPDH expression. The sequence of the primers used for qPCR analysis were as followed: GAPDH-R (human) (5'-AAG TGG TCG TTG AGG GCA ATG-3').

ELISA Kit of cytokines secretion
THP-1 and U937 cells were stimulated with H 2 O 2 , as previously described, and the resulting supernatants were collected after a 48-h treatment. We quantified the release of the senescence-associated secretory phenotype (SASP) IL-8, CCL5, and SERPINE1 cytokines (Elder and Emmerson 2020;Rossi et al. 2023) in response to H2O2 stimulation using commercially available ELISA kits (IL-8 Cat#BY-EM230561, CCL5 Cat#BY-EM220320, and SERPINE1 Cat#BY-EM220493) obtained from BOYAN Biotech in China. The absorbance of the plates was measured at 450 nm, and the analysis was performed using InfiniteM200 (Tecan) according to the manufacturer's protocol. The levels of cytokines released by THP-1 and U937 cells were determined by calculating the plate standard curve.

Statistical analyses
In this research, statistical analyses were performed using the R software tool (version 4.1.1). For continuous variables, the Wilcoxon signed-rank test was employed to compare two groups, while the Kruskal-Wallis H-test was employed for comparison when the group number > 2. The correlation of categorical variables was estimated by the Pearson chi-square test. The visualization of overall statistical analyses results was performed through the ggplot2 R package. The ROC curve was drawn by using the pROC package. All analyses were a two-sided comparison, with an alpha level of 0.05. Benjamini Hochberg's method was used for multiple hypothesis tests to control the false-discovery rate (FDR).

Three senescent-macrophage subtypes identified by non-negative matrix factorization
The flowsheet of the present study is displayed in Fig. 1. This study acquired somatic mutation and transcriptome profiles of 406 BLCA patients from TCGA after excluding samples without clinicopathological or follow-up information. The Fig. 3 Investigation of different CNV and SNP profiles in the MSR subtypes. A Heatmap showed the top 15 frequently mutated genes in 3 subtypes, respectively. The key at the bottom of the plot represents each mutation type includes amplification, missense mutations, and deep deletions. The title of the waterfall plot showing the mutation patient number. B The co-occurrence or mutually exclusive landscape for each subtype (ns no significance, *P < 0.05, **P < 0.01). C The CNV amplification or deletion aberrations of each sample in MSR subgroups. The results were standardized by Log10 and conducted by the Wilcoxon test (C1:47, C2:109, C3:243; ns no significance, *P < 0.05, **P < 0.01, ***P < 0.001). D Circos performed CNV frequency regions among different clusters ◂ researchers identified 279 senescence-related genes from CellAge and 1044 macrophage cell markers (Supplementary  Table 1, Figure S1), and extracted 23 overlapping genes for further analysis. Using the NMF algorithm with the "brunet" standard, the selected MSR Bca genes were clustered, revealing three distinct modification pattern clusters: pattern cluster 1 included 243 cases, pattern cluster 2 included 109 cases, and pattern cluster 3 included 47 cases. The accuracy and robustness of the clustering results was determined by three main factors cophenetic, dispersion, and silhouette coefficient ( Fig. 2A-B). The Kaplan-Meier survival analysis showed that cluster 3 had a significantly better survival rate compared to cluster 1 and cluster 2 among the macrophagesenescence subtypes found in the TCGA dataset (Fig. 2C). Similarly, as demonstrated in GSE19423 array (Fig. S2), C3 showed a statistically significant survival difference compared with C1 while no significant difference was observed between C3 and C2. To analyzed the heterogeneity of immune infiltration in the three subtypes, we use XCELL, MCPCOUNTER, CIBERSORT, and EPIC algorithms to quantify the immune microenvironment of BLCA samples. The different immune infiltration pattern was exhibited in patients with 3 MSR subgroups, C3 had a higher abundance of activated Tregs cells, T cells follicular helper, and dendritic cells, but a lower abundance of T cells CD4 naive and B cells naïve compared with C2 ( Figure S3). In particular, among monocyte macrophage infiltration, the results showed that the abundance of infiltrating macrophages was higher in the C1 than in the C2 and C3 (Fig. 2D). We observed that C1 had a higher proportion of M0 stage and M1 stage cells, whereas a slightly lower infiltration abundance of M2 stage cells, indicating that C1 could be identified as proinflammatory phenotype (Fig. 2E-F). Simultaneously, we utilized TIDE algorithm to predict the potential immunotherapeutic response. The result implied C3 suffered worse immunotherapeutic response rate than C1, which showed C1 subgroups had poor prognosis and were positively correlated with higher proportion of immunotherapy responders (Fig. 2G). Furthermore, we employed differential gene analysis to identify enriched immune HALLMARKER pathways among the three subtypes and explored the functional signatures intergroup via GSVA ( Figure S2). The immunerelated pathways, such as IL2-STAT5 and complement system, were predictably enriched in C1. Taken together, MSR subtypes displayed distinct characteristics where C1 had a higher senescent-macrophage infiltration ratio and lower overall survival, whereas C3 demonstrated negative immunotherapeutic response due to the highly expressed immunerepressive genes.

Construction of MSR-related signature and functional analyses
Initially, each gene expression from the training set underwent screening through univariate Cox regression with a P value threshold of less than 0.05. Subsequently, 23 genes with statistically significant correlations to overall survival (OS) were identified. Using univariate and multivariate Cox algorithms, we developed a primary prognostic risk model through LASSO regression with an optimal lambda value of 0.02 and 10-fold cross-validation (Fig. 4A-B). The threemarker signatures model genes were jointly determined using these dimension-reduction techniques ( Figure S4). The eventual formula as follows: MSR Score = 0.0068*(TGFB1I1)-0.00023*(ID1) + 0.00092*(AKR1B1). The notable KEGG pathways (Fig. 4C) and GO terms (Fig. 4E) were illustrated by the GSVA method according to low-risk and high-risk TCGA-BLCA samples. The results implied that high-MSR score group was mainly associated with senescence-related GO terms, including regulatory senescence, cellular senescence, and inhibitory activity. Furthermore, WNT signaling Fig. 4 The construction and functional analyses of the risk model. A-B Lasso regression identified 3 MSR genes associated with overall survival. C GSVA results of KEGG pathways and D HALLMARKS gene sets, E GO terms based on the differentially expressed genes within the TCGA Array. F Differential immune-related genes expression and immunomodulators pattern between MSR values in BLCA patients. FDR false-discovery rate, LogFC log2(fold change) ◂ and tyrosine metabolism KEGG pathways were enriched in high-risk samples, while P53 signaling was enriched in lowrisk samples. In any case, the hallmark gene sets derived from the MsigDB website were identified using singlesample Gene Set Enrichment Analysis (GSEA). The results demonstrated that high-risk samples showed enrichment for hedgehog signaling, complement pathway, and inflammatory response (Fig. 4D). These enrichment figures indicated that macrophage-senescence-related pathways were active in the high-risk group. Furthermore, we conducted an evaluation of immune-related gene expression levels across two risk groups (Fig. 4F). Our analyses revealed that low-MSR values may indicate stronger immune activities. The high-risk group demonstrated downregulation of positive immunomodulatory factors, which were highly correlated with impaired anti-tumor immunity.

Validation of the risk model
The research validated the robustness of the risk signature by training on a specific set and assessing on three external datasets from GEO (see Supplementary Table 2). The Kaplan-Meier survival plots, utilizing the best cutoff of MSR determined by X-tile, effectively stratified BCa patients into high and low stages based on overall survival in the TCGA-BCa training cohort (P = 0.001, n = 406), the GSE13507 cohort (P = 0.024, n = 165), the GSE31684 cohort (P = 0.005, n =90), and the GSE32894 cohort (P < 0.001, n = 221). Our findings indicated that BLCA patients with high MSR had a lower overall survival rate across all four datasets (Fig. 5A-B). Principal component analysis satisfactorily classified all samples into two groups based on MSR (Fig. 5C). Additionally, we conducted a meta-analysis of hazard ratio values for MSR across all four datasets, which comprehensively evaluated the potential utility of the risk model in clinical practice (Fig. 5D).

Real-world immunotherapeutic responses and potential drug sensitivity
To explore the efficacy of predicting PD1/PD-L1 treatment responses, we employed IMvigor210 cohort as validation set (Mariathasan et al. 2018). All patients were identified into low-and high-MSR score groups based on the optimal threshold. The Kaplan-Meier survival plot indicated that patients with higher levels of MSR scores had a worse prognosis compared to those with lower levels (P-value < 0.05; Fig. 6A). The study then investigated the association between low-or high-MSR scores and tumor cells (TC) level, tumor-infiltrating immune cells (IC) level, and immunologically excluded phenotype (Fig. 6B). Patients with higher MSR scores showed a higher correlation with an immunologically excluded phenotype in BCa patients, as well as fewer IC and TC levels compared to those with lower MSR scores. Additionally, patients who responded favorably to immunotherapy had higher MSR scores, particularly for significant signatures such as AKR1B1 (P-value < 0.05) and TGFB1I1 (P-value < 0.05; Fig. 6C).
Subsequently, we evaluated potentially effective antitumor drugs by analyzing their half maximal inhibitory concentration (IC50) in GDSC databases to uncover the relationship between the MSR score and drug sensitivity. Out of the total 672 possible drugs, we identified the top 15 that were most negatively correlated with the MSR score. Then, we summarized the specific relationship between those LN IC50 and risk score using Fig. 6D. Based on the GDSC database, the experimental treatment linsitinib (He et al. 2022) (P < 0.001), an inhibitor of INSR and IGF1R, and RSK inhibitor SL0101 (Chelakkot et al. 2020) (P < 0.001) exhibited lower AUC at the high-MSR group than low-MSR group (Fig. 6F). Generally, five of chemotherapeutic agent cisplatin, vinblastine, methotrexate, gemcitabine, and doxorubicin were utilized as neoadjuvant or adjuvant treatment. We also analyzed the potential efficacy among the MSR subgroups and each of IC50 value according to CCLE database. The Wilcox test indicated that IC50 of cisplatin, methotrexate, and vinblastine ( Fig. 6E) had significant negative correlation with low MSR, which were demonstrated their effectiveness for high-MSR samples.

SC-RNA analysis in tumor environment
The macrophage barcodes were extracted from an integrated single-cell RNA dataset among various cells in the TME and divided into two subgroups based on the middle MSR value (Fig. 7A). The part of differentially expressed genes (DEGs) between high-MSR and low-MSR subgroups were depicted as Fig. 7B. It was obviously implied that high score group had high expression of TGFB1I1, CEBPB (Montes et al. 2021), and low expression of NUPR1 (Grasso et al. 2015), most of which were senescence-related factors. The GSEA analysis of Reactome annotation showed ROS activity pathway, glycolysis, and senescence-associated secretory phenotype were involved in high-MSR cells compared to low (Fig. 7D). This finding suggests that tumor-associated Fig. 5 Prediction of biomarkers model's prognostic value. A Kaplan-Meier survival method indicated that the patients with higher MSR risk scores exhibited worse prognosis in the TCGA training cohort, GSE13507 cohort, GSE32894 cohort, and GSE31684 cohort. B Distribution of MSR score according to the survival time and status in 4 cohorts. C Scatter diagram for PCA analysis based on the MSR in training and external cohorts. D The meta forest-analysis revealed that robustness of MSR model validated by TCGA-BLCA, GSE13507, GSE31684, GSE32894 (due to the high heterogeneity between results, a random effect model is chosen as the final criterion for judgment) ◂ macrophages (TAMs) could be effectively classified into distinct cell groups based on their MSR scores. Pseudotime analysis was conducted on the subpopulations of high-MSR and low-MSR cell types. Consistent with previous studies, low-MSR macrophages were predominantly located at the root, while high-MSR macrophages were located on the branches (Fig. 7C).
Tumor-associated macrophages with low or high MSR keep critical role in cooperation with other TME cells to promote tumor growth or immune suppression. We employed CellChat package to delineate the cellular interactions and receptor-ligand pairs (LR pairs) in tumor tissue ( Figure S5). This model quantified the regulatory pathways and expression levels of LR pairs between TME cells and senescent macrophages in different bladder cancers. From the overall pairs among cell communications, we found that CXCL receptors and the SPP1 pathway showed significant heterogeneity in cellular interaction for high-and low-MSR cells (Fig. 7E). Consistent with previous studies, cellular senescence induces the activity of SASP complement IL-1α, IL-6, Mmp-3, Mmp-9, Cxcl-1, and Cxcl-10 pathways , which might affect the proliferation of tumor epithelial cells via cell-cell communication. The specific regulatory genes of signaling between TME cells and each macrophage group in bladder tumors are exhibited in Fig. 7F. Taken together, these results demonstrate that the cellular communication network and functional differences between TME cells and MSR subgroup macrophages could affect the tumor immune microenvironment and proliferation.

Validation in vivo and vitro
In order to further confirm the expression of three senescence-related signatures in a model, we employed tumor tissue specimens from 27 patients with bladder cancer (BCa) (12 patients at stage I-II and 15 at stage III-IV) to verify differential expression of these signatures via immunohistochemical staining. The results showed that transforming growth factor beta 1 (TGFB1I1) and aldo-keto reductase family 1 member B (AKR1B1) had higher protein expression of stage III-IV patients than stage I-II patients. In contrast, inhibitor of DNA binding 1 (ID1) had lower expression in advanced BCa patients (Fig. 8A). Integrated optical density (IOD) represented the density of target protein-positive cells; the boxplots were evaluated by Welch's T-test (Fig. 8B). To determine whether TGFB1I1, AKR1B1, and ID1 expression is involved in the senescent-macrophage process, we treated THP-1 cells and U937 macrophage cell lines with the senescence inducer H 2 O 2 for 24 h (Fig. 8C). The CCK8 assay results showed that cell viability of both cell lines decreased continuously after 3 to 15 mM of culture, whereas no significant difference with each other (Fig. 8D). And the IC50 values of H 2 O 2 with THP-1 and U937 were 8.6 mM and 7.8mM (Stepanić et al. 2019), respectively. Meanwhile, we observed alterations in pro-inflammatory cytokines and SASP factors in senmacrophages subjected to hydrogen peroxide stimulation (Rossi et al. 2023). Measurement of the levels of IL-8, CCL5, and SERPINE1 by multiplex ELISA exhibited a significant increase in the treatment group, thus confirming the establishment of senescence phenotype in THP-1 and U937 cells following H 2 O 2 exposure ( Figure S6). Furthermore, RT-PCR analysis showed that the relative mRNA level across control group and H 2 O 2 -treated cell lines below its IC50 value. The result supported, as expected, senescent cell lines had higher AKR1B1 and TGFB1I1 level than control group, while ID1 expression was inversely lower than the control group (Fig. 8E). These outcomes indicated that novel macrophage-senescence-related signatures would be effective in vivo and vitro.

Discussion
Immunosenescence plays a crucial role in promoting the development of tumors and age-related diseases, directly reducing the proportion of immune cells and impairing the immune response capacity within the tumor microenvironment. Additionally, cytotoxic T lymphocytes and dendritic cells infiltration levels, which have an anti-tumor effect, are also compromised due to cellular senescence (Feehan et al. 2021;Liu et al. 2023). Researchers have found that intercellular signals from different macrophage phenotypes can adjust T lymphocyte response activity . For instance, the M1 phenotype macrophage, which is proinflammatory, mostly polarizes to pro-tumorigenic M2 phenotype or immunosuppressive M2-like macrophages elevate in the aging-tumor microenvironment. Anti-tumor strategies for elder cancer patients have focused on therapeutic agents targeting the regulation of TAM polarization . Hence, it is essential to construct macrophage-senescence-related biomarkers in the TME, while the relative signature model still remains unknown in BLCA.  6 IMvigor210 dataset for immunotherapeutic responses and GDSC drug sensitivity. A Overall survival difference of low-and high-MSR subgroups in immunotherapeutic response dataset. B The IC level, TC level, and immune clinical features of the IMvigor210 data tested by chi-square, IC0 depicts < 1%, IC1 depicts ≥ 1% but < 5%, IC2+ depicts ≥ 5% in PD-L1 positive patients, TC0 depicts < 1%, TC1 depicts ≥ 1% but < 5%, TC2 depicts ≥ 5% in specimens of PD-L1-positive patients. C The Kruskal-Wallis test displayed the relationship of the risk score between SD/PD group and CR/PR group. D The LN IC50 of 15 potential effective drug extracted from GDSC database were negatively correlated with MSR, tested by the Wilcox method. E The violin plots revealed IC50 difference of cisplatin, methotroxate, and vinblastine in the high/low-MSR group. F AUC values of linsitinib and SL0101 in the high/low-MSR group. Label: ns no significance, *P < 0.05, **P < 0.01, ***P < 0.001 ◂ Page 15 of 18 228 We identified 23 common macrophage-senescencerelated genes in BLCA patients using SC-RNA cell markers from the CELLAGE and TCGA databases. The genes were clustered using the NMF method and their clinical prognosis and immune infiltration were analyzed. The three MSR subtypes displayed significantly different functional fields. C1 showed a higher level of M1 macrophage infiltration and a better immunotherapeutic response compared to the other subtypes. The results of SNV landscape implied that C1 had HMCN1 mutation feature while C3 had RYR2 and MUC16 mutation. It was reported that age-related macular degeneration (AMD) had been proved to be positively correlated with HMCN1 mutation in the elder (Nudelman et al. 2023), and the patients with RYR2 mutation had superior anti-tumor immunity and better prognosis in esophageal adenocarcinoma . Previous findings suggest that they can suppress the harmful release of mitochondrial reactive oxygen species and inhibit oxidative damage, thus delaying the progression of macrophage senescence (Hamilton et al. 2020). Du et al. (2019) reported that upregulation of TGFB1I1 plays a crucial role in constructing a tumor-promoting stroma and leads to a senescence-like phenotype. However, depletion of BAG3 induces G1 arrest and cellular senescence through the accumulation of p27 (Chillappagari et al. 2022).
Multivariate Cox regressions and Lasso iteration were assembled to identify a 3-gene risk model for assessing OS and RNA average expression originating from the TCGA cohort. We applied the MSR score along with datasets GSE13507, GSE32684, and GSE31684 external Fig. 7 Subpopulations of MSR validation and cellular interactions in TME. A The t-SNE plot for high-and low-MSR cells extracted from integrated data GSE145137 and GSE135337. B DEGs of TAMs between opposite subgroups exhibited as scatterplot. C The pseudotime analysis of differentiation trajectories with subpopulations of MSR score. D The GSEA results of senescence-related Reactome enrichment in different MSR subgroups. E Summary of CellChatDB ligand-receptor pairs among high-, low-MSR cells, and TME cells in bladder tumor tissues. F The specific regulatory genes in CXCL signaling among TAMs and all tumor cells ◂ Fig. 8 Immunohistochemical staining and RT-PCR for MSR signature validation. A specimens from local patients at TNM stages, the protein expression of AKR1B1, ID1, and TGFB1I1 was compared with the stage I-II and stage III-IV. B The boxplot summarized by IOD (ns no significance, *P < 0.05, **P < 0.01, ***P < 0.001). D Revealed by CCK-8 assay, THP-1, and U937 cell viability when treated with H 2 O 2 for 3 to 15 (ns no significance, *P < 0.05, **P < 0.01, ***P < 0.001, mM = μmol/L). E The qPCR analysis of 3 MSR signature markers between control group and senescence inducer-treated group (ns no significance, *P < 0.05, **P < 0.01, ***P < 0.001, points represent the means, the error bars lines represent standard deviations) data, consisting of 165 primary BCa samples, 221 bladder urothelial carcinoma samples, and 90 cases of malignant tumors, respectively. We explored the signature's association with the immunotherapy response by referencing the IMvigor210 immunotherapeutic-related cohort. Moreover, prediction analysis of potential target drugs revealed that high MSR had a negative correlation with some small molecule drugs, indicating a possible subsequent treatment regimen. Using the Wilcoxon signed-rank test with a P-value threshold of < 0.05, we discovered that the MSR signature significantly influenced the immune infiltration level in the tumor microenvironment (TME). We also conducted GSEA analysis to classify macrophage-like cells into high-and low-MSR cell groups and studied intercellular communication across multiple types of tumor cells to gain insight into malignant cell-cell communication heterogeneity. Evidence strongly suggests that CXCL-ASKR LR pairs play a pivotal role in monocyte recruitment signaling intercellularly and immunosuppressive M2 phenotype differentiation with the senescence environment (Orecchioni et al. 2019). Our findings reveal molecular patterns at the level of LR pairs that are characteristic of the senescent immune cell phenotype. These results present a potential therapeutic strategy for targeting selective macrophage senescence. Additionally, we compared the pathological differential expression of the selected model genes using immunohistochemical sections extracted from 27 local BLCA patients. In vitro, we cultured THP-1 and U937 senescent-macrophage cell lines with oxidant H 2 O 2 and included RT-PCR analysis (Zhou et al. 2022). As expected, the results of experimental validation were consistent with the prediction of model parameters. Altogether, the MSR score was an efficient model with superior traits of survival assessment.
We appointed some putative biomarkers, AKR1B1, ID1, and TGFB1I1, as a risk signature through the present study. According to Tanawattanasuntorn et al. (2022), AKR1B1 plays a role in increasing cellular oxidative stress in breast and ovarian cancer cells, and its upregulation is strongly linked to tumor migration aggressiveness, suggesting its involvement in the senescence TME process. TGFB1I1 and ID1 are inhibitors of transforming growth factor and DNA binding, respectively. The former serves as a tumor suppressor gene in colorectal cancers (CRCs) and pancreatic cancers (Pcas), and its depletion can increase tumor cell growth and reduce tumor cell apoptosis (Luo et al. 2022). ID1 exhibits dual roles in anti-tumor effects. Although ID1 has been shown to increase EPC angiogenesis within TME through PI3K/AKT in ovarian cancer, other studies have demonstrated its function in inhibiting proliferation, which is regulated by PGC1α, in various tumor types. Therefore, ID1 is considered a novel prognostic marker and an ideal chemotherapeutic target (Fei et al. 2023;Oh et al. 2021). In general, the specific regulatory pathway of these biomarkers warrants further clarification based on the established signature, which contributed to the identification of new biomarkers in BLCA.
The current study is subject to certain limitations. Firstly, the sequencing depth and coverage of the two single-cell RNA datasets retrieved were inadequate. This restricted further exploration of additional molecular factors and tumor-associated macrophages within the senescence tumor microenvironment due to limited cell barcodes. Additionally, as noted in the immunotherapy data cited in this paper, some patients who received immune checkpoint blockade (ICB) may have also undergone chemotherapy or radiotherapy prior to treatment, thereby compromising the homogeneity of the samples.

Conclusions
In conclusion, our study has established the senescencerelated macrophage signature within the tumor microenvironment, as demonstrated by the results of multiple independent cohorts, including SC-RNA analysis and bulk-seq data. Utilizing local clinical samples and in vitro cell analysis, we have developed a reliable MSR risk score model that enables assessment of clinical overall survival, immunotherapeutic response, and immune microenvironment infiltration in patients with BLCA