Abstract
Immunosenescence refers to the immune system undergoing a series of degenerative changes with advancing age and is tightly associated with the initiation and progression of cancers. However, the immunosenescence-related genes as critical biomarkers for bladder cancer (BLCA) have not been systematically analyzed. We retrieved the immunosenescence-related genes from the public database and verified their association with hallmarks of immunosenescence based on The Cancer Genome Atlas (TCGA) cohort. Through gene pairing, Lasso, and univariate Cox regression, an 8-gene pair model was constructed to evaluate the overall survival of BLCA, which was then validated in the training cohort (), two external validation cohorts (; ), and local samples (). We also downloaded the clinical information and gene expression matrices of other 32 different cancers from TCGA. The established model showed significant predictive value for the prognosis in 15 cancers (). The risk model could also serve as a promising predictor for immunotherapeutic response, which has been verified by the TIDE algorithm (), IMvigor210 dataset (), and other two datasets correlated with immunotherapy (; ). The TCGA dataset, in vitro cell experiments, and pan-cancer analysis displayed that the gene signature was associated with cisplatin sensitivity (). Overall, we proposed a novel immunosenescence-related gene signature to predict prognosis, immunotherapeutic response, and cisplatin sensitivity of BLCA, which were validated in different independent cohorts, local samples, and pan-cancer analyses.
1. Introduction
Bladder cancer (BLCA) is one of the predominant malignancies with high mortality and morbidity worldwide and is characterized by poor prognosis and high recurrence rates [1]. BLCA is mainly divided into nonmuscle invasive bladder cancer and muscle-invasive bladder cancer, which has less than 50% overall survival (OS) rates [2], based on the degree of cancer invasion. Tremendous progress in the treatment of BLCA has been made in recent years. Bacillus Calmette-Guerin (BCG) intravesical perfusion has been utilized for the treatment of BLCA since 1976 and serves as traditional management for the moment [3], demonstrating that BLCA is immunogenic. Recent studies display that BLCA is among the malignancies with the highest tumor mutational burden (TMB), and TMB was reported as the strongest predictor for immunotherapeutic effectiveness [4, 5]. On the whole, BLCA is prone to respond to chemotherapeutic reagents, including the newly proposed immune checkpoint inhibitors (ICIs) [6]. In addition to immunotherapy, neoadjuvant chemotherapy, especially cisplatin-based neoadjuvant chemotherapy, which has been recommended for the first-line treatment before radical cystectomy, dramatically improves the prognosis of BLCA [7]. Despite the advances in medical therapy, a large number of patients remain unresponsive to immunotherapy or chemotherapy and suffer unfavorable prognoses. Hence, seeking more accurate and reliable predictions for the prognosis and sensitivity to immunotherapy and chemosensitivity is a hot issue of much interest in BLCA.
With the rise and popularization of high-throughput sequencing, more and more biomarkers associated with the prognosis or drug response of BLCA have been discovered [8, 9]. The screened biomarkers not only helped to evaluate the clinical outcomes but also provided the clues to study the potential mechanisms of BLCA, such as ferroptosis [10], autophagy [11], and N6-methyladenosine (m6A) modification [12]. Nevertheless, the identification of the novel biomarkers with high accuracy and robustness is still meaningful, which could serve as promising clinical tools to guide the personalized treatment and the cut-in points to expound the correlated biological processes.
Most malignancies exhibit higher incidence rates in the elderly compared with the young, which is probably due in part to immunosenescence. The conception of immunosenescence was proposed by Walford in 1969 and referred to the immune system, including immune organs and immune cells, degenerated with aging [13]. The main hallmarks of immunosenescence include thymic involution, decreased T cell population and diversity, decreased antigen presentation ability of dendritic cells, attenuated phagocytosis of macrophages, and increased myeloid-derived suppressor cells (MDSCs), causing the decline of immune surveillance and cancer cell clearance [14]. The decreased immune surveillance promotes tumor immune escape, which directly leads to the proliferation, metastasis, and drug resistance of cancer cells. Although immunosenescence plays an essential role in cancer initiation and progression, the immunosenescence-related biomarkers have not been systematically analyzed in BLCA.
The present study collected the immunosenescence-related genes from the public database and verified their association with the hallmarks of immunosenescence in The Cancer Genome Atlas-Bladder Cancer (TCGA-BLCA) cohort, which was also set as the training dataset. We adopted a gene pair strategy to construct the risk signature to render the model applicable in different detection platforms [15]. Lasso, univariate Cox, and multivariate Cox regression analyses were used to identify the hub gene pairs associated with the OS in BLCA. The predictive ability of the risk model to immunotherapeutic and chemotherapeutic response was also detected. GSE13507, GSE32894, GSE5287, IMvogor210, GSE35640, GSE78220, pan-cancer analysis, in vitro cell experiments, and local samples were used for external validation.
2. Materials and Methods
2.1. Data Source and Processing
The senescence-related genes were retrieved from the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/), and the immune-related genes were obtained from ImmPort (https://www.immport.org/), as shown in Supplementary Tables 1 and 2, respectively. The gene sets of the main hallmarks of immunosenescence, including senescence-associated secretory phenotype (SASP), mitochondrial biogenesis, glycolysis, and cellular response to reactive oxygen species (ROS), were also downloaded from MSigDB [16]. The transcriptional data with fragments per kilobase per million mapped reads (FPKM) format, corresponding clinical information, and genomic mutation data (varscan software) of the TCGA-BLCA cohort were downloaded from TCGA (https://tcga-data.nci.nih.gov/tcga/, accessed on August 20, 2021). The maftools package in R was used to calculate the TMB of each sample and to visualize the Top 30 genes with the highest mutational rate in the different subgroups. The RNA expression matrices and OS data of the other 32 cancers were also obtained from TCGA. GSE13507 [17], GSE32894 [18], GSE5287 [19], GSE35640 [20], and GSE78220 [21] were directly downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/, accessed on August 21, 2021). The RNA expression data with FPKM format and clinical features of the IMvogor210 cohort were extracted from the IMvigor210CoreBiologies package in R software (version 3.6.3). We utilized R to conduct quality control, and the genes with average and the cases with <30 days’ follow-up were excluded.
2.2. Clinical Sample Collection
The BLCA tissues from 10 patients undergoing partial/radical cystectomy were collected between November 2019 and July 2021 in the Third Affiliated Hospital of Southern Medical University. All the fresh BLCA samples were collected from the center of the tumor during surgery and stored in liquid nitrogen for RNA extraction. The patients undergoing preoperative chemotherapy, immunotherapy, and radiotherapy were not included in this study. The project has been approved by the Ethics Committee of the Third Affiliated Hospital of Southern Medical University (ID: 20211121), and the informed consent files were signed by all patients. The diagnosis of BLCA depends on the histopathological examination, and the detection of the tumor tissues’ tumor-node-metastasis (TNM) stages was based on the eighth TNM staging system defined by the American Joint Commission on Cancer. Among the 10 patients, 1 patient was in TNM stage I, 4 patients were in TNM stage II, 3 patients were in TNM stage III, and 2 patients were in TNM stage IV.
2.3. Single-Sample Gene Set Enrichment Analysis (ssGSEA)
The ssGSEA was performed with the GSVA and the GSEABase packages with method specifications as “ssgsea,” “Gaussian,” and “abs.ranking=TRUE”, and the ssGSEA -score of each sample was calculated to quantify the enrichment level. The collected immunosenescence-related genes were set as the reference gene set.
2.4. Risk Model Construction
To make the risk model applicable in different detection platforms, we adopted the gene pair method to construct the model. We defined a gene combination form, “gene A | gene B,” as a gene pair. If the expression of gene A was higher than that of gene B, the pair would be considered 1; otherwise, it would be regarded as 0. When the amount of the gene pairs with values equal to 1 or 0 accounts for less than 20% of the total pairs, the pairs would be excluded. All the genes were cyclically singly paired, and a 0-or-1 matrix was established. Lasso and univariate Cox regression were performed to identify the significant gene pairs associated with OS of BLCA through the glmnet and the survival packages. In the univariate Cox analysis, was thought to be significant. Subsequently, the gene pairs con-determined by Lasso and univariate Cox regressions were included in the multivariate Cox regression analysis with stepwise by the survminer package, and the risk model was ultimately constructed. The risk score evaluated with the risk model was defined as immunosenescence-related score (IRS). IRS was calculated using the following formula:
2.5. Validation of the Risk Model
The Kaplan-Meier survival analysis with log-rank test was performed to detect the prognosis difference between high- and low-IRS patients through the survival package, and was regarded to be significant. The optimal cut-off value was determined through X-tile software. The time-dependent receiver operating curves (ROCs) were drawn via the survivalROC package to evaluate the predictive ability of IRS to 1-, 3-, and 5-year OS rates and to compare the predictive ability of IRS and the risk clinicopathological traits, including age, gender, grade, stage, and TNM stage, after the transformation of all the continuous variables into binary variables.
2.6. Functional Enrichment Analysis
The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional annotation were conducted via the clusterProfiler package in R with and filtering. The enrichment analyses were visualized with the enrichplot and the ggplot2 packages. The Gene Set Enrichment Analysis (GSEA) was performed in GSEA desktop software (version 4.1.0, https://www.gsea-msigdb.org/gsea/) with 1,000 permutations and default parameters. The gene sets with a and were considered to be significant [22].
2.7. The Analyses of Immune Infiltration and Immunotherapeutic Response
The infiltration level of different immune cells among the samples from the TCGA-BLCA project was evaluated via CIBERSORT, which calculated the infiltration proportion of 22 kinds of immune cells [23]. The infiltration immune cell content difference between the high- and low-IRS groups was measured with the Wilcoxon signed-rank test. We also utilized the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm to predict the response to ICIs of the patients from the TCGA-BLCA cohort [24] and conducted the Chi-squared test to detect the association with the IRS stratification. The IMvigor210 cohort, containing 298 patients with metastatic urothelial cancer (UC) undergoing atezolizumab (an anti-PD1 agent) treatment, were also used for external validation through Kaplan-Meier survival analysis, Wilcoxon signed-rank test, and Chi-squared test. GSE35640, including 56 patients with melanoma, and GSE78220, including 27 patients with melanoma, were selected to verify the predictive ability of IRS to immunotherapy in other malignancies via the Wilcoxon signed-rank test. The filtering threshold was set to .
2.8. The Evaluation of Chemotherapeutic Response
The sensitivity to the common chemotherapeutic reagents, such as cisplatin, doxorubicin, gemcitabine, methotrexate, and vinblastine, of the samples from the TCGA-BLCA project was quantified as half inhibitory concentration (IC50) through the pRRophetic package. Besides, the IRSs among the patients with different clinical statuses, including complete response (CR), partial response (PR), progression disease (PD), and stable disease (SD), were also compared. The Wilcoxon signed-rank test was adopted for difference detection, and was significant.
2.9. Cell Culture and Treatment
Human BLCA cell line T24 was purchased from the Stem Cell Bank, Chinese Academy of Sciences (Shanghai, China), and maintained in McCoy’s 5A Medium (Gibco, USA), which has been added with 100 U/ml penicillin, 100 μg/ml streptomycin, and 10% fetal bovine serum (Gibco, USA), in a humidified atmosphere with 5% CO2 at 37°C. The cells were treated with 20 μM cisplatin (Sigma-Aldrich, USA) for 24 hours before total RNA extraction [25].
2.10. Real-Time Quantitative PCR (RT-qPCR)
The total RNA of the human and cell samples was isolated with the TRIzol/chloroform method (TRIzol, Invitrogen, USA) after the homogenization of human samples, following the manufacturer’s instruction. The cDNA was synthesized with the PrimeScript RT reagent Kit (Takara, Japan), and the RT-qPCR was performed with the SYBR Premix Ex Taq II reagent (Takara, Japan) based on the Applied Biosystems 7300 Real-Time PCR System (Applied Biosystems, USA). GAPDH was chosen as the interreference to normalize the data, and the 2-ΔΔCt method was used to calculate the mRNA expression value. The statistical analysis was performed with Student’s -test. The primers utilized in this study are shown in Supplementary Table 3.
2.11. Statistical Analysis
The statistical analyses of the present study were conducted in R software (version 3.6.3). The data were presented as the or (%). The Wilcoxon signed-rank test was used to compare the difference in different groups if not otherwise specifically stated. The Pearson Chi-squared test was performed to analyze the association of categorical data, which was visualized via the ggplot2 package. The ROCs as the performance evaluation method for binary classification ability were drawn via the pROC package. .
3. Results
3.1. Identification of the Immunosenescence-Related Genes
The workflow of the present study is shown in Figure 1. First, 105 common genes con-determined by the 590 senescence-related genes from MSigDB and 1793 immune-related genes from ImmPort were identified as the potential immunosenescence-related genes (Figure 2(a)). The ssGSEA was implemented to calculate the enrichment level of the retrieved immunosenescence-related gene set among the 396 BLCA samples from the TCGA-BLCA project (Figure 2(b)). It was found that the BLCA samples with high ssGSEA -score exhibited worse OS in the TCGA-BLCA cohort (, Figure 2(c)), implying that the collected gene set was significantly associated with the prognosis. The prognosis value of the gene set was also validated in the patients with BLCA from the GSE5287 cohort (, Figure 2(d)) and GSE32894 cohort (, Figure 2(e)). GSEA indicated that SASP (Figure 2(f)), cellular response to ROS (Figure 2(g)), and glycolysis (Figure 2(h)) were positively associated with the ssGSEA -score, while mitochondrial biogenesis (Figure 2(i)) was negatively associated with the ssGSEA -score. The increase of SASP, cellular response to ROS, and glycolysis and the decrease of mitochondrial biogenesis have been reported as the main hallmarks of immunosenescence [16]. Besides, the association between the ssGSEA -score and the known immunosenescence marker genes was also detected (Figure 2(j)). Overall, the immunosenescence-related genes collected from the public databases could reflect the immunosenescence process to some extent and were chosen for further analysis.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
3.2. Risk Model Construction
After the cyclical pairing, a total of 1533 gene pairs were established. Lasso regression identified that 40 of 1533 gene pairs were significantly associated with the OS (Figures 3(a) and 3(b)). At the same time, 31 gene pairs were determined by univariate Cox regression with filterings. 17 gene pairs were con-determined by these dimension-reduction methods (Figure 3(c)), 8 of which were ultimately included in the risk model via multivariate Cox analysis with stepwise (Figure 3(d)). The IRS was calculated as follows: , where represented a gene pair. According to the optimal cut-off value detected by X-tile, which was equal to 2.90, all the patients in the training cohort and the external validation cohorts were divided into the high- and low-IRS subgroups. The Sankey plot displayed the association among IRSs, ssGSEA -scores, and the survival statuses in the TCGA-BLCA project (Figure 3(e)). Figures 3(f) and 3(g) indicated the GO and KEGG functional annotation of the 15 genes comprising the risk signature, respectively, where some important pathways both associated with immunosenescence and tumorigenesis, such as chemical carcinogenesis-reactive oxygen species and PD-L1 expression and PD-1 checkpoint pathway in cancer, were significantly enriched.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
3.3. IRS Is a Robust Biomarker for Prognosis Prediction
We validated the robustness of the risk signature in the training dataset, 2 external validation datasets from the public database, local BLCA samples, and pan-cancer analysis. The baseline information of the training dataset (TCGA-BLCA) and the 2 external validation dataset (GSE13507 and GSE32894) is displayed in Table 1. Kaplan-Meier survival analyses showed that the BLCA patients with high IRS suffered a lower OS rate in the TCGA-BLCA cohort (, Figure 4(a)), GSE13507 cohort (, Figure 4(b)), and GSE32894 cohort (, Figure 4(c)). The time-dependent ROCs for the evaluation of 1-, 3-, and 5-year OS in TCGA-BLCA (Figure 4(d)), GSE13507 (Figure 4(e)), and GSE32894 (Figure 4(f)) verified the predictive value. Besides, with the increase of IRS, more deaths were observed, as shown in the scatter plots (Figure 4(g)–4(i)). The patients with high IRS in the TCGA-BLCA (Figure 4(j)), GSE13507 (Figure 4(k)), and GSE32894 (Figure 4(l)) were more likely to exhibit immunosenescence statuses via the GSEA. The expression level of the 15 genes in the risk model in BLCA samples and adjacent normal samples is shown in Supplementary Figure 1, and their predictive performance for the OS of BLCA in the TCGA-BLCA cohort (Supplementary Figure 2), GSE13507 cohort (Supplementary Figure 3), and GSE32894 cohort (Supplementary Figure 4) was measured via Kaplan-Meier survival analyses.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
Next, the tumor samples of the 10 patients with BLCA from the local hospital were collected for experimental validation. The 10 patients were divided into the TNM stage I-II and TNM stage III-IV groups, and the gene expression value was measured via RT-qPCR. The expression difference of the 15 genes in the risk signature between TNM stage I-II and TNM stage III-IV is shown in Supplementary Figure 5a, and the statistical method was Welch’s -test. According to the detected gene expression value, the IRS of each patient was calculated (Supplementary Figure 5b). The patients in TNM stage III-IV have higher IRSs compared with those in TNM stage I-II through Welch’s -test (, Supplementary Figure 5c), which partly proved the reliability of the risk model.
To confirm whether the established signature could serve as a predictor for the prognosis in pan-cancer, we downloaded the datasets of other 32 cancers from TCGA, and the detailed information is shown in Table 2. As shown in Figure 5, the IRS was a significant biomarker for the prognosis in 15 different cancers, implying that IRS was a powerful predictor in pan-cancer.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
3.4. The Clinical Association of IRS
The Chi-squared test indicated that IRS was significantly associated with age (), gender (), tumor grade (), TNM stage (), and pathological T stage () based on the TCGA-BLCA project, as shown in Figure 6(a). After transforming all the factors into binary variables, we compared the predictive ability of the clinicopathological parameters and IRS to OS using ROC analyses. The optimal cut-off of age was detected by the X-tile software, and the patients over 64 years old suffered the most significantly worse OS rates among the subjects in the TCGA-BLCA cohort (). The areas under curves (AUCs) of IRS were all higher than those of the risk clinical features in the 1-year (, Figure 6(b)), 2-year (, Figure 6(c)), 3-year (, Figure 6(d)), 4-year (, Figure 6(e)), and 5-year (, Figure 6(f)) OS evaluation. IRS was also an independent risk factor for OS among the BLCA patients from the TCGA-BLCA cohort through the univariate (hazard ratio, HR; confidence interval, CI; , -7.17, ) and the multivariate (, 95% -6.61, ) analyses (Table 3).
(a)
(b)
(c)
(d)
(e)
(f)
3.5. IRS Is a Promising Tool to Evaluate the Immunotherapeutic Response
Compared with the low-IRS group, the high-IRS group exhibited a higher infiltration level of resting memory CD4 T cells (), M0 macrophages (), M2 macrophages (), activated mast cells (), and neutrophils () and the lower infiltration proportion of naïve B cells (), memory B cells (), CD8 T cells (), activated memory CD4 T cells (), follicular helper T cells (), and activated NK cells (), indicating the tremendous changes of tumor immune microenvironment between high- and low-IRS patients (Figure 7(a)). Besides, the low-IRS patients harbored higher TMB (, Figure 7(b)), suggesting that the low-IRS subjects were more likely to benefit from immunotherapy, and the TIDE algorithm verified the assumption (, Figure 7(c)). Next, we validated the predictive ability of IRS in different independent cohorts receiving immunotherapy. It was found that the patients with metastatic UC showing high IRS exhibited poorer survival rates in the IMvigor210 cohort (, Figure 7(d). Compared with the subjects with no response to atezolizumab, the patients with high sensitivity to atezolizumab had lower IRSs via the Wilcoxon signed-rank test (, Figure 7(e)) and Chi-squared test (, Figure 7(f)). In addition to the BLCA patients, 2 melanoma cohorts who had received immunotherapy, GSE35640 and GSE78220, were also used for validation. Figures 7(g) and 7(h) displayed that the patients with resistance to immunotherapy had higher IRSs, redemonstrating the potential of IRS to evaluate the immunotherapy sensitivity.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
The expression difference of the 15 genes in the risk signature between the responsive subjects and nonresponsive subjects from the IMvigor210 cohort is shown in Supplementary Figure 6a. Supplementary Figure 6b displays the predictive value of each variable to the OS of the IMvigor210 cohort.
3.6. IRS Can Predict Cisplatin Sensitivity
A total of 109 patients receiving chemotherapy were extracted from the TCGA-BLCA project, and IRS could still discriminate against the high-risk patients (, Figure 8(a)). Figure 8(b) indicates that IRS was significantly associated with the IC50 of cisplatin () and methotrexate (). However, the Kaplan-Meier survival plots showed that IRS was not a significant prognosis predictor among the patients receiving methotrexate treatment (, Figure 8(c)), but IRS could significantly distinguish the high-risk subjects receiving cisplatin treatment (, Figure 8(d)). The IRS level of the patients with CR, PD, PR, and SD to cisplatin treatment is shown in Figure 8(e), and a significant difference between the CR and PR groups was observed (). However, it should be stated that most cases in the TCGA-BLCA cohort received chemotherapeutic agent combination therapy. Therefore, a series of in vitro cell experiments to reconfirm the association was conducted, and the T24 cells with cisplatin treatment were used for external validation. As shown in Figure 8(f), the level of ADIPOR2 (), CCN2 (), CTSS (), GBP2 (), IRF1 (), MAP2K1 (), TFRC (), and THBS1 () was obviously increased after cisplatin treatment, while the level of EGFR (), ELAVL1 (), MAPK1 (), NOX4 (), and SRC () was markedly decreased, suggesting that most of the genes in the risk signature have a tight relationship with the pharmacology of cisplatin.
(a)
(b)
(c)
(d)
(e)
(f)
Pan-cancer analysis was also used for validation, and the other 8 types of cancer, which contained >20 patients receiving cisplatin treatment in the TCGA project, were selected (Table 2). IRS showed powerful predictive capability for the prognosis among the cases receiving cisplatin treatment with mesothelioma (MESO, , Figure 9(a)), head and neck squamous cell carcinoma (HNSC, , Figure 9(b)), ovarian serous cystadenocarcinoma (OV, , Figure 9(c)), and stomach adenocarcinoma (STAD, , Figure 9(d)). Figures 9(e)–9(h) displayed the predictive ability of IRS in cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC, ), lung squamous cell carcinoma (LUSC, ), lung adenocarcinoma (LUAD, ), and testicular germ cell tumors (TGCT, ), respectively, and the heterogeneity and limited sample size might account for the nonsignificance. Overall, IRS was a reliable predictor for cisplatin sensitivity, which has been validated in different cohorts and cell experiments.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
3.7. The Mutational Landscape in High- and Low-IRS Patients
Since the IRS was tightly associated with TMB, we then explored the mutational divergence between high- and low-IRS patients based on the TCGA-BLCA project. The Top 30 genes with the highest mutational rate in the high- and low-IRS subgroups are shown in Figures 10(a) and 10(b), respectively. We found that IRS could predict the mutational statuses of CNTN2 (, Figure 10(c)), FAM129A (, Figure 10(d)), FGFR3 (, Figure 10(e)), KIAA1257 (, Figure 10(f)), NCF1 (, Figure 10(g)), SIGMAR1 (, Figure 10(h)), and STAG2 (, Figure 10(i)) with high efficacy, implying that the mutation of these genes was correlated with IRS.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
4. Discussion
With the population aging, the prevalence and mortality of BLCA increase gradually in many regions. Immunosenescence is accepted as one of the leading causes of tumorigenesis, which directly promotes the evasion of immune surveillance and the immunosuppressive atmosphere in tumor microenvironment [26, 27]. With the advance of immunosenescence, the infiltration level of dendritic cells, T cells, M1 macrophages, and N1 neutrophils decreases, while the infiltration abundance of the immunosuppressive cells, such as Treg cells, M2 macrophages, and N2 neutrophils, increases (Figure 11). The depletion of functional T cells, the attenuation of antigen-presenting ability, and the accumulation of immunosuppressive cells lead to unfavorable drug treatment effectiveness and poor prognosis [28]. Hence, the identification of the immunosenescence-related biomarkers is important and meaningful, but no immunosenescence-related gene model has been reported in BLCA to date.
In this study, we retrieved the immunosenescence-related genes from the MSigDB and ImmPort and verified the association of the collected genes with the known hallmarks of immunosenescence via ssGSEA and GSEA. A sum of 105 immunosenescence-related genes was screened. Subsequently, we adopted a gene pair method to construct the risk signature, and the 105 genes were cyclically singly paired. Lasso, univariate Cox, and multivariate Cox regressions identified an 8-gene pair model to evaluate the OS based on the TCGA-BLCA project. Here, we defined the risk score calculated by the model as IRS. GSE13507, which included 165 BLCA samples, GSE32894, which included 224 BLCA samples, 10 BLCA samples collected from the local hospital, and pan-cancer analysis were used for external validation to confirm the prognostic value of IRS. The Wilcoxon signed-rank test displayed that IRS was significantly associated with TMB and immune infiltration level in the tumor microenvironment, which was calculated via CIBERSORT. Therefore, we then explored whether IRS could serve as a predictor for the immunotherapeutic response. The validation in the TCGA-BLCA, IMvigor210 cohort, GSE35640 cohort, and GSE78220 cohort proved the assumption. IRS was also a marker for cisplatin treatment, which has been validated in the cohort receiving cisplatin treatment from the TCGA-BLCA project, T24 cells treated with cisplatin, and pan-cancer analyses. The mutational landscape in the high- and low-IRS groups was also detected.
The more accurate and reliable prediction for the prognosis and treatment response is always one of the hottest issues in BLCA. Many predictive gene models have been proposed, providing useful clinical tools for clinicians and BLCA patients [29, 30]. However, the number of the risk model used for the evaluation of prognosis, immunotherapeutic response, and chemosensitivity at the same time is limited to date. A predictive model applied to multiple different scenarios at the same time is more practical and more portable. Besides, we adopted the gene pair strategy to construct the signature, which means that the risk model does not depend on the definite gene expression value and could be applied for different detection platforms, including RT-qPCR, RNA-seq, microarray, and NanoString. Overall, the proposed model in this study was a promising and practical tool for BLCA and can be used to guide personalized treatment.
Some genes were first reported as biomarkers for BLCA. For instance, GBP2 was significantly upregulated in the BLCA samples with low malignancy (Supplementary Figure 5), and the high expression GBP2 heralded a favorable prognosis (Supplementary Figures 4 and 6). Godoy et al. reported that GBP2 was significantly associated with the favorable prognosis in 766 patients with breast cancer, implying that GBP2 was a tumor suppressor [31]. In general, the established signature helped to identify novel biomarkers in BLCA.
However, some shortcomings should not be neglected. First, this study was designed as a retrospective and multicenter research due to the limited financial support, and a prospective, randomized, and double-blind clinical trait would be more beneficial to verify the clinical usefulness. Second, experimental validation, such as transfection assays and tumor-forming experiments in athymic mice, should be conducted to clarify the biological functions of the screened novel biomarkers.
5. Conclusions
To sum up, a novel immunosenescence-related gene signature was developed to estimate the prognosis, immunotherapeutic response, and cisplatin sensitivity in BLCA, which has been validated in different independent cohorts, local BLCA samples, in vitro cell experiments, and pan-cancer analyses.
Data Availability
The data used to support the findings of this study are available from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/), GEO (https://www.ncbi.nlm.nih.gov/geo/), MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/), and ImmPort (https://www.immport.org/), and the code would be supplied from the corresponding author upon request.
Conflicts of Interest
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Authors’ Contributions
Ranran Zhou and Jingjing Liang contributed equally to this work.
Acknowledgments
Special thanks are due for the contribution of TCGA, GEO, MSigDB, and ImmPort databases. This study is supported by the National Natural Science Foundation of China (NO. 81772257), Youth Cultivation Program of Southern Medical University (NO. PY2018N076), and Medical Scientific Research Foundation of Guangdong Province (NO. A2019557).
Supplementary Materials
Supplementary 1. Supplementary Table 1: the collected senescence-related genes.
Supplementary 2. Supplementary Table 2: the collected immune-related genes.
Supplementary 3. Supplementary Table 3: the primer sequence utilized in this study.
Supplementary 4. Supplementary Figure 1: the expression difference of the 15 genes in the risk signature between adjacent normal and BLCA tissues in the TCGA-BLCA cohort (a), GSE13507 cohort (b), and GSE32894 cohort (c).
Supplementary 5. Supplementary Figure 2: the predictive value of the 15 genes to OS in the TCGA-BLCA cohort.
Supplementary 6. Supplementary Figure 3: the predictive value of the 15 genes to OS in the GSE13507 cohort.
Supplementary 7. Supplementary Figure 4: the predictive value of the 15 genes to OS in the GSE32894 cohort.
Supplementary 8. Supplementary Figure 5: validation of the prognostic value of IRS in the local cohort. (a) The expression level of the 15 genes in the risk signature in the subjects with BLCA in stage I-II and the subjects with BLCA in stage III-IV. (b) The calculated IRS of each patient. (c) The patients in stage III-IV exhibited higher IRS levels compared with those in stage I-II.
Supplementary 9. Supplementary Figure 6: the expression divergence of the 15 genes in the signature between the patients with high immunotherapeutic sensitivity and the patients with low immunotherapeutic sensitivity (a) and the prognostic value of each gene in the IMvigor210 cohort (b).