A circadian rhythm-related gene signature for predicting relapse risk and immunotherapeutic effect in prostate adenocarcinoma

Jin Liu; Zhao Tan; Shijie Yang; Xinda Song; Wenping Li

doi:10.18632/aging.204288

Research Paper Volume 14, Issue 17 pp 7170—7185

A circadian rhythm-related gene signature for predicting relapse risk and immunotherapeutic effect in prostate adenocarcinoma

Jin Liu^1, , Zhao Tan^1, , Shijie Yang^1, , Xinda Song^1, , Wenping Li^1, ,

¹ Department of Urology, The Third Hospital of Hebei Medical University, Shijiazhuang, China

Received: June 9, 2022 Accepted: September 5, 2022 Published: September 13, 2022

https://doi.org/10.18632/aging.204288
How to Cite

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

Abstract

Prostate adenocarcinoma (PRAD) represents the most common male carcinoma in developed countries, its high relapse risk contributes to the second-leading cause of cancer-related deaths. Therefore, it is required to develop an effective signature for predicting the relapse risk of PRAD. To identify a circadian rhythm- (CR-) related predictive signature, we analyzed RNA-seq data of patients with prostate adenocarcinoma (PRAD) from the TCGA and GEO cohort. Seven circadian rhythm- (CR-) related genes (FBXL22, MTA1, TP53, RORC, DRD4, PPARGC1A, ZFHX3) were eventually identified to develop a CR-related signature. AUCs for 3-year overall survival were 0.852, 0.856 and 0.944 in the training set, validation set and an external independent test set (GSE70768), respectively. Kaplan-Meier curve analysis showed that the high-risk group has a reduced relapse-free survival (RFS) than the low-risk group in the training set, validation set, and test set, respectively (P < 0.05). We constructed a prognostic nomogram combining the CR-related signature with T staging to precisely estimate relapse risk of PRAD patients. Finally, we observed that the CR-related gene signature was associated with tumor mutation burden, multiple immune checkpoint molecules and microsatellite instability, and thus could predict response to immune checkpoint inhibitors in PRAD. Conclusively, we developed a circadian rhythm-related gene signature for predicting RFS and immunotherapy efficacy in prostate adenocarcinoma.

Introduction

Prostate adenocarcinoma (PRAD) represents the most common male carcinoma in developed countries [1]. In spite of substantial efforts invested into therapeutic development of PRAD, its high relapse risk contributes to the second-leading cause of cancer-related deaths [2–4]. Therefore, it is clinically necessary to identify a useful signature for PRAD to guide cancer treatment. Presently, the most common clinical indicators of cancer relapse are including prostate-specific antigen (PSA), Gleason score, and tumor stage [5, 6]. However, due to the heterogeneity of PRAD, the predictive ability of these common indicators is not unmet. Thereby, to identify a novel biomarker for predicting relapse risk for PRAD is an urgent and relevant effort.

Circadian rhythms are 24-hour oscillations that affect multiple biological functions in humans [7]. Circadian rhythm disorders are linked to aggressive tumor behaviors and unwanted clinical outcomes. Circadian-related genes have been implicated in the pathogenesis of colorectal cancer [7], prostate cancer [8], and bladder cancer [9]. Meanwhile, emerging evidence points to its involvement in tumor microenvironment [10–12]. While circadian rhythm is a hot topic of cancer research recently, the specific mechanisms of its role in humans are unclear. Besides, it is unreported whether circadian rhythm- (CR-) related gene signature can be a prognostic biomarker for PRAD patients.

To establish a CR-related predicting signature for PRAD patients, we investigated bulk RNA sequencing (RNA-seq) profiles from the Cancer Genomic Atlas (TCGA) and Gene Expression Omnibus (GEO), hoping provide an applicable gene signature for predicting prognosis for PRAD patients.

Materials and Methods

Acquisition of data

Gene expression profiling and clinical information for PRAD patients from TCGA were obtained from UCSC Xena on Sep 9, 2021 (https://xena.ucsc.edu/). Microarray RNA-seq data and survival information of 29 PRAD patients was retrieved from GSE70768 in GEO. GSE70768 was sequenced using the platform of GPL10558 (Illumina HumanHT-12 V4.0 expression beadchip). RNA-seq data from TCGA and GSE70768 were normalized in the form of TPM values and then log2(x + 1) transformed.

Somatic nucleotide variation (SNV) data of PRAD patients were obtained from the TCGA database.

Estimation of enrichment scores for individual patients

To quantify the expression levels of CR gene set in individual patients, we estimated enrichment score (ES) of the CR-related gene set for individual PRAD patients using single-sample gene set enrichment analysis (ssGSEA) [13]. ssGSEA is a mathematic methodology to estimate relative expression levels of a given gene set using RNA-seq data. The parameters used in this study were as follows: min.sz = 1, max.sz = Inf, tau = 0.25.

Construction of the gene signature

338 PRAD patients from TCGA were randomly divided into the training set (n = 236) and the validation set (n = 102). The CR-associated genes were screened for eligible genes for establishing the predictive signature using univariate Cox regression firstly, and then further analyzed using least absolute shrinkage and selection operator (LASSO) regression. Then, the eligible genes in LASSO were utilized to construct a gene signature based on eligible genes’ expression levels and their corresponding coefficients in LASSO, using the following formula: FBXL22×(-0.3746637) + MTA1×(0.9070002) + TP53×(-0.2111043) + RORC×(-0.4931651) + DRD4×(0.4730237) + PPARGC1A×(-0.2315328) + ZFHX3×(-0.5221173).

Assessment of the predictive performance of the gene signature

The predictive ability of the gene signature was assessed mainly by two analyses including receiver operating characteristic (ROC) curve and Kaplan-Meier (KM) curve. Area under the curve (AUC) and log-rank test in were performed in the training set, the validation set, and an independent test set.

Functional enrichment analysis

Functional enrichment analysis was carried on using the R package “clusterProfiler” (version: 3.18.1) [14]. This toolkit can determine whether canonical biological processes and signaling pathways are significantly enriched in a given patient cohort based on gene expression profiles. The information of canonical biological functions and signaling pathways are available in Bioconductor annotation data GO.db and KEGG.db.

Construction of a nomogram

To better predict the prognosis of PRAD patients, we conjointly analyzed the gene signature and several common clinical characteristics using multivariate Cox regression analysis and established a predicting nomogram using R package ‘rms’.

Statistics

Statistical analysis was carried on using R software (Version 4.0.1). Independent sample t test or Wilcoxon signed rank test was utilized according to the homogeneity of variance and normal distribution of data. Spearman’s correlation coefficient was performed to investigate the relationship between two continuous variables. Statistical significance was considered when P value is less than 0.05.

Availability of data and materials

All data generated and described in this article are available from the corresponding web servers, and are freely available to any scientist wishing to use them for noncommercial purposes, without breaching participant confidentiality. Further information is available from the corresponding author on reasonable request.

Results

Establishment of a circadian rhythm- (CR-) related gene signature

Circadian rhythm (CR) is reported to be implicated in cancer [7], whereas it remains unclear whether it has an effect on prostate adenocarcinoma (PRAD). To investigate its association with PRAD, we compared the expression levels of circadian rhythm signaling pathway between PRAD and normal tissues using gene expression profiles of PRAD patients from TCGA. The expression levels of circadian rhythm signaling pathway was quantified as enrichment score (ES) using ssGSEA algorithm based on RNA-seq data of 551 samples from the TCGA cohort of PRAD patients, and the findings showed that ES was significantly increased in normal tissue than in tumor tissue (Figure 1A). Moreover, we divided PRAD patients into the low-ES and high-ES groups according to a median value of ES, and performed a survival analysis. High-ES patients had an improved relapse-free survival (RFS) than low-risk patients (Figure 1B). Gene set enrichment analysis (GSEA) also supported that CR signaling pathway was significantly enriched in normal tissue compared with tumor tissue (Figure 1C). These findings indicated that circadian rhythm had a relationship with PRAD.

Figure 1. Establishment of a circadian rhythm- (CR-) related gene signature in prostate adenocarcinoma (PRAD). (A) Enrichment score (ES) of CR-related gene set was significantly enriched in normal tissue than tumor tissue. (B) High-ES patients had an improved relapse-free survival (RFS) than their counterparts. (C) CR-related gene set was positively enriched in normal tissue compared with tumor tissue. (D) 13 CR-related genes were qualified in univariate Cox regression analysis (FBXL22, FBXL6, MTA1, NR2F6, TP53, BTRC, GHRL, RORC, CIPC, DRD4, MTOR, PPARGC1A, ZFHX3; P<0.05). (E, F) 13 qualified genes were further filtered using LASSO regression analysis to eliminate multicollinearity and seven eligible genes (FBXL22, MTA1, TP53, RORC, DRD4, PPARGC1A, ZFHX3) were eventually acquired for establishment of a CR-related gene signature for predicting RFS in PRAD patients.

Following the finding of the relationship between CR and PRAD, we wondered if CR could predict prognosis for PRAD patients and then planned to develop a circadian rhythm- (CR-) related gene signature to predict survival of PRAD patients. To identify eligible CR-related genes in PRAD, we first performed univariate Cox regression analysis for 117 CR-related genes, and obtained 13 genes (FBXL22, FBXL6, MTA1, NR2F6, TP53, BTRC, GHRL, RORC, CIPC, DRD4, MTOR, PPARGC1A, ZFHX3; P<0.05; Figure 1D). Then, 13 eligible genes were further filtered using LASSO regression analysis to eliminate multicollinearity and seven eligible genes (FBXL22, MTA1, TP53, RORC, DRD4, PPARGC1A, ZFHX3) were eventually acquired for establishment of a CR-related gene signature for predicting RFS in PRAD patients (Figure 1E, 1F). The CR-related gene signature was quantified based on the mRNA expression levels and the corresponding coefficients of seven CR-related genes, using the following formula: FBXL22×(-0.3746637)+MTA1×(0.9070002)+TP53×(-0.2111043)+ RORC×(-0.4931651)+DRD4×(0.4730237)+ PPARGC1A×(-0.2315328)+ ZFHX3×(-0.5221173). The coefficients represented the influence of genes on relapse risk; the positive represented a risk factor for relapse, while the negative represented a protective factor for relapse.

Assessment of predicting performance of the gene signature

The predicting capacity of the CR-related gene signature was validated using ROC curve and KM curve in the training set (n = 236) and the validation set (n = 102). Area under the curve (AUC) for predicting 3-year overall survival were 0.852 and 0.856 in the training set and validation set, respectively (Figures 2A–2B). Then, we quantified the risk of relapse for PRAD patients as risk score using the above mentioned formula, and divided patients into the low- and high-risk patients according to the median value of risk score. Consistent with the above results, survival analysis also revealed that low-risk group had an improved RFS than high-risk group both in the training set and in the validation set (log-rank test, P <0.001).

Evaluation of the performance of the gene signature. (A) AUC for predicting three-year RFS is 0.852 in the training set (n=236). (B) AUC for predicting three-year RFS is 0.856 in the validation set (n=102). (C) Low-risk patients had an improved RFS than high-risk patients in the training set. (D) Low-risk patients had an improved RFS than high-risk patients in the validation set.

Figure 2. Evaluation of the performance of the gene signature. (A) AUC for predicting three-year RFS is 0.852 in the training set (n=236). (B) AUC for predicting three-year RFS is 0.856 in the validation set (n=102). (C) Low-risk patients had an improved RFS than high-risk patients in the training set. (D) Low-risk patients had an improved RFS than high-risk patients in the validation set.

Discriminative ability of the CR-related signature

To investigate the discriminative ability of the CR-related signature, we performed principal component analysis (PCA) using seven selected CR-related genes for PRAD patients, and observed a distinction between the low- and high-risk patients (Figure 3A, 3C). Then, we explored the association among gene signature, relapse status and seven gene expression levels. The results demonstrated that risk score was linked to relapse status and the expression of seven genes (Figure 3B, 3D).

The discriminative power of the CR-related signature. (A) Principle component analysis (PCA) demonstrated that low-risk group was apparently distinct from high-risk group in Dim 1 in the training set. (B) As risk score increased, chance of tumor recurrence increased in the training set. (C) Principle component analysis (PCA) showed that the low-risk group was apparently distinct from the high-risk group in Dim 1 in the validation set. (D) As risk score increased, chance of tumor recurrence increased in the validation set.

Figure 3. The discriminative power of the CR-related signature. (A) Principle component analysis (PCA) demonstrated that low-risk group was apparently distinct from high-risk group in Dim 1 in the training set. (B) As risk score increased, chance of tumor recurrence increased in the training set. (C) Principle component analysis (PCA) showed that the low-risk group was apparently distinct from the high-risk group in Dim 1 in the validation set. (D) As risk score increased, chance of tumor recurrence increased in the validation set.

Comparison of predicting performance of the CR-related gene signature with other indicators

To further evaluate the predicting performance of the CR-related gene, we compared the predicting capability of the CR-related gene with other indicators of RFS, including clinical characteristics and other four reported gene signatures. The results showed that CR-related gene signature showed an improved predictive performance than other clinical indicators (PSA, Gleason, clinical T stage, clinical N stage), with an AUC of 0.831 (Figure 4A). Meanwhile, CR-related gene signature (Signature 1) also showed an improved predictive performance than other four reported gene signatures: the immune-related gene signature [15] (Signature 2), the metabolism-related gene signature [16] (Signature 3), the PPP1R12A-related gene signature [16] (Signature 4), the TMB-related gene signature [17] (Signature 5), with an AUC of 0.856 VS 0.763, 0.689, 0.626 and 0.488 (Figure 4B).

Comparison of the CR-related gene signature with other indicators for RFS in PRAD. (A) CR-related gene signature showed an improved predictive performance than other clinical indicators, with an AUC of 0.831. (B) CR-related gene signature (Signature 1) showed an improved predictive performance than other four reported gene signatures. (C) CR-related gene signature also showed an ideal predictive performance in an independent cohort of PRAD (GSE70768), with an AUC of 0.944. (D) Low-risk group had an improved RFS than high-risk group in an independent cohort of PRAD (GSE70768) (with a cutoff point of median value; log-rank test, P = 0.003).

Figure 4. Comparison of the CR-related gene signature with other indicators for RFS in PRAD. (A) CR-related gene signature showed an improved predictive performance than other clinical indicators, with an AUC of 0.831. (B) CR-related gene signature (Signature 1) showed an improved predictive performance than other four reported gene signatures. (C) CR-related gene signature also showed an ideal predictive performance in an independent cohort of PRAD (GSE70768), with an AUC of 0.944. (D) Low-risk group had an improved RFS than high-risk group in an independent cohort of PRAD (GSE70768) (with a cutoff point of median value; log-rank test, P = 0.003).

Furthermore, we assessed the predicting ability of CR-related gene signature in an independent cohort of PRAD patients (GSE70768) using ROC curve. Surprisingly, the CR-related gene signature manifested an impressive predictive ability, with AUC of 0.944 (Figure 4C). Consistently, we also quantified the risk of RFS for patients using the above mentioned formula and divided patients into the low- and high-risk group based on a cutoff point of median value. Low-risk group had a significantly improved RFS than high-risk group (log-rank test, P = 0.003; Figure 4D). These results further underlined the predicting capacity of the CR-related gene signature for predicting RFS in PRAD patients.

Functional enrichment analysis

To interrogate the biological functions that were linked to the CR-related gene signature, we carried out functional enrichment analysis for genes that correlated with the CR -related gene signature. We calculated the correlation coefficients between the gene signature and all genes. 695 genes were selected as the signature-related genes (P < 0.01, R > 0.4) [18]. Afterwards, these 695 genes were investigated using R package “clusterProfiler”. The results demonstrated that these genes were mainly enriched in cytoplasmic translation, oxidative phosphorylation, mitochondrial respiratory chain complex assembly, oxidoreduction-driven active transmembrane transporter activity and electron transfer activity (Figure 5A–5D). Gene set enrichment analysis (GSEA) displayed that several enriched biological processes were associated with cancer invasion and metastasis, including AMPK signaling pathway, dopaminergic synapse, central carbon metabolism in cancer, NF-kappa B signaling pathway and others (Figure 5E, 5F).

Figure 5. Functional enrichment analysis for CR-related gene signature. (A) Enriched biological processes included cytoplasmic translation, oxidative phosphorylation, mitochondrial respiratory chain complex assembly. (B) Enriched cell components included mitochondrial protein-containing complex, ribosome, inner mitochondrial membrane protein complex. (C) Enriched molecular functions included structural constituent of ribosome, oxidoreduction-driven active transmembrane transporter activity and electron transfer activity. (D) Enriched KEGG pathways included ribosome, diabetic cardiomyopathy, chemical carcinogenesis (E, F) Gene set enrichment analysis (GSEA) showed the top 10 KEGG signaling pathway including AMPK signaling pathway, dopaminergic synapse, central carbon metabolism in cancer, NF-kappa B signaling pathway and others.

Establishment of a nomogram

For better prediction of the survival of the PRAD patients, we endeavored to combine the signature with clinical features to establish a prognostic nomogram. We first analyzed the ability of age, clinical T stage, clinical N stage, Gleason score, prostate specific antigen (PSA) and the CR-related gene signature (risk score) to predict RFS using multivariate Cox regression. The results demonstrated that the gene signature and clinical T stage remained to be a valid predictor (Figure 6A, P≤0.05).

Clinical application of CR-related gene signature for PRAD patients. (A) Multivariate Cox regression was used to investigate the ability of age, clinical T stage, clinical N stage, Gleason score, prostate specific antigen (PSA) and the CR-related gene signature (risk score) to predict RFS. Risk score remained to be a valid predictor (PB) A nomogram combing risk score and clinical T stage was constructed to predict 1-, 3- and 5-year RFS for individual PRAD patient. (C) ROC analysis for nomogram showed an impressive predictive performance, with AUC of 0.944. (D) Calibration curve showed agreement between actual and predicted RFS, indicating an ideal predictive capability.

Figure 6. Clinical application of CR-related gene signature for PRAD patients. (A) Multivariate Cox regression was used to investigate the ability of age, clinical T stage, clinical N stage, Gleason score, prostate specific antigen (PSA) and the CR-related gene signature (risk score) to predict RFS. Risk score remained to be a valid predictor (P<0.05). (B) A nomogram combing risk score and clinical T stage was constructed to predict 1-, 3- and 5-year RFS for individual PRAD patient. (C) ROC analysis for nomogram showed an impressive predictive performance, with AUC of 0.944. (D) Calibration curve showed agreement between actual and predicted RFS, indicating an ideal predictive capability.

Then, we established a prognostic nomogram consisting of the gene signature and pathological T (Figure 6B). To assess prognostic ability of the nomogram, we performed a ROC curve and calibration curve. ROC curve showed an amazing predictive capability, with 5-year AUC of 0.944; consistently, C-index of the nomogram was 0.75±0.05, indicating the discriminatory ability of the nomogram; calibration curve also displayed a consistency between predicted probability and observed probability, suggesting the effectiveness and robustness of the established nomogram.

Exploration of the association of CR-related gene signature with tumor immune microenvironment

Since CR-related gene signature could predict RFS in PRAD, we next wondered to explore its relationship with tumor immune microenvironment. Tumor immune microenvironment has been reported to be implicated in identification of pivotal molecules, drug response, and novel therapeutic methods. To interrogate the role of CR-related gene signature in PRAD, we computed the proportion of immune cells in cancerous tissue using CYBERSORT and investigated the association between risk score and immune cells. The results showed that T cell follicular helper was significantly elevated, while B cell naive, B cell plasma, T cell CD4 memory resting and monocyte were reduced in high-risk group than low- risk group (P < 0.05; Figure 7A). Further correlation analysis supported that T cell follicular helper was positively correlated with gene signature (Figure 7E), while B cell naive, B cell plasma, T cell CD4 memory resting and monocyte were negatively correlated with gene signature (Figure 7B–7D, 7F).

Association of the CR-related gene signature with tumor infiltrating immune cells in PRAD. (A) Comparison of tumor infiltrating immune cells between low- and high-risk groups demonstrated that there existed a significant difference in the abundance of B cell naive, B cell plasm, T cell CD4 memory, T cell follicular helper, and monocyte (P B–D) Risk score was inversely correlated with B cell naive, B cell plasm, T cell CD4 memory, T cell follicular helper, and monocyte (R P E) Risk score was positively correlated with T cell CD4 memory (R P F) Risk score was inversely correlated with monocyte (R P

Figure 7. Association of the CR-related gene signature with tumor infiltrating immune cells in PRAD. (A) Comparison of tumor infiltrating immune cells between low- and high-risk groups demonstrated that there existed a significant difference in the abundance of B cell naive, B cell plasm, T cell CD4 memory, T cell follicular helper, and monocyte (P < 0.05). (B–D) Risk score was inversely correlated with B cell naive, B cell plasm, T cell CD4 memory, T cell follicular helper, and monocyte (R < 0, P < 0.05). (E) Risk score was positively correlated with T cell CD4 memory (R < 0, P < 0.05). (F) Risk score was inversely correlated with monocyte (R < 0, P < 0.05).

The profiling of somatic nucleotide variation for PRAD patients between the low- and high-risk patients

To investigate the association of CR-related gene signature with the mutation levels in PRAD patients, we profiled the mutation landscape of the low-risk and high-risk groups, respectively, by analyzing somatic nucleotide variation data of PRAD patients from the TCGA cohort using the maftool R package. To better compare the potential distinction of the mutation landscape between different risk-score groups, we here defined the patients with the top 25% of risk score as the high-risk group and defined the patients with bottom 25% of risk score as the low-risk group. The waterfall plot demonstrated that the high-risk group had a higher nucleotide variation rate than the low-risk group (50.6% vs 35.8%, Figure 8A, 8B). Consistent with the findings of the waterfall plot, the bar plot and box plot analyses also showed that the risk scores were critically increased in the high-mutation group than in the low-mutation group (Figure 8C, 8D). A high somatic nucleotide variation rate has been reported to be associated with tumor mutation burden (TMB) that is reflective of effective response to immune checkpoint inhibitors (ICIs) [19, 20]; Thus, we asked whether the risk score had a relationship with TMB. We first quantified tumor mutation burden (TMB) using the maftool R package based on somatic nucleotide variation data of PRAD patients, and found that TMB was significantly correlated with risk score (P = 0.02, Figure 8E). Collectively, these results showed that the risk score could reflect the genomic instability and the TMB level, indicating its potential capability to predict response to ICIs in PRAD patients.

The profiling of somatic nucleotide variation for PRAD patients of different risk groups. (A, B) The waterfall plot of nucleotide variation rate in high-risk group and low-risk group. (C, D) The bar plot and box plot analyses of risk scores in the high-mutation group and low-mutation group. (E) The correlation of tumor mutation burden and risk score. *P P P

Figure 8. The profiling of somatic nucleotide variation for PRAD patients of different risk groups. (A, B) The waterfall plot of nucleotide variation rate in high-risk group and low-risk group. (C, D) The bar plot and box plot analyses of risk scores in the high-mutation group and low-mutation group. (E) The correlation of tumor mutation burden and risk score. *P <0.05, ** P <0.01, ***P<0.001.

The effects of CR-related gene signature on response to ICIs

Elevated TMB, high-microsatellite instability (MSI) and immune checkpoint molecules have been demonstrated to be predictive of response to immune checkpoint inhibitors (ICIs) in cancer patients [21]. The previous results had suggested that CR-related risk score had a significantly impact on RFS, genomic stability and TMB; thereby we asked if it was also associated with MSI status and immune checkpoint molecules in PRAD. To enhance the reliability of results, we quantified MSI levels using three different strategies, including ssGSEA, PreMSIm and UCSCXenaShiny. MSI levels was significantly elevated in high-risk group than in low-risk group (P <0.001; Figure 9A), and positively correlated with risk score (P < 0.05, R = 0.37; Figure 9B). Consistently, risk score was also significantly higher in MSI-H group than MSI-L/MSS group using PreMSIm method (Figure 9C; P <0.001). Meanwhile, risk score was significantly higher in MSI-H group than MSI-L/MSS group using UCSCXenaShiny method (P <0.05; Figure 9D).

The effects of CR-related risk score on response to ICIs. (A) MSI levels was elevated in the high-risk patients than in low-risk patients (P B) MSI levels were positively correlated with risk score (P R = 0.37). (C) Risk score was significantly increased in MSI-H group than MSI-L/MSS group (data of MSI generated from PreMSIm R package; P D) Risk score was significantly higher in MSI-H group than MSI-L/MSS group (data of MSI generated from UCSCXenaShiny R package; P E) Comparison of multiple immune checkpoint mRNA (PDCD1, LAG3, CD40, CTLA4, PDCD1LG2) in the high-risk patients and the low-risk patients. *P P P

Figure 9. The effects of CR-related risk score on response to ICIs. (A) MSI levels was elevated in the high-risk patients than in low-risk patients (P <0.001). (B) MSI levels were positively correlated with risk score (P < 0.05, R = 0.37). (C) Risk score was significantly increased in MSI-H group than MSI-L/MSS group (data of MSI generated from PreMSIm R package; P < 0.001). (D) Risk score was significantly higher in MSI-H group than MSI-L/MSS group (data of MSI generated from UCSCXenaShiny R package; P < 0.05). (E) Comparison of multiple immune checkpoint mRNA (PDCD1, LAG3, CD40, CTLA4, PDCD1LG2) in the high-risk patients and the low-risk patients. *P <0.05, ** P <0.01, ***P<0.001.

Consistently, the expression levels of LAG3, PDCD1 and CTLA4 were also critically upregulated in the high-risk patients than in the low-risk patients (Figure 9A), underlining the ability of CR-related gene signature to predict drug response (Figure 9E). Collectively, these findings further underlined the ability of CR-related gene signature to predict immunotherapeutic response in PRAD.

Discussion

We established a seven-gene CR-related signature for predicting relapse risk in PRAD patients, developed a nomogram consisting of the CR-related gene signature and T staging. Several enriched biological processes were related to AMPK signaling pathway, dopaminergic synapse, central carbon metabolism in cancer, and NF-kappa B signaling pathway. Moreover, we observed that the CR-related gene signature was associated with tumor mutation burden, multiple immune checkpoint molecules and microsatellite instability, and thus could predict response to immune checkpoint inhibitors in PRAD. Overall, this study provided a novel predictor of relapse risk and immunotherapeutic effects of patients with PRAD, which would help to develop precision medicine.

One contribution of this research is the development of a CR-related gene signature for predicting clinical outcomes of PRAD patients. The predicting ability of the signature has been validated in an independent cohort (GSE70768) and compared with common clinical index and some reported gene signatures to prove the superiority of its predicting ability.

Another important finding that several important biological processes were identified to be involved in CR and high relapse risk, including AMPK signaling pathway, dopaminergic synapse, central carbon metabolism in cancer, and NF-kappa B signaling pathway. Consistent with our findings, AMPK signaling pathway has been reported to involve in PRAD [22–24], and targeting AMPK signaling pathway with CO or metformin can suppress prostate cancer cell growth [25, 26]. Remarkably, it is unreported that whether dopaminergic synapse is related to PRAD. Here, we observed the association between dopaminergic synapse and PRAD, which requires further investigation. Central carbon metabolism has a key role in metabolic programming [27], and regulation of carbon metabolism is essential for the treatment of prostate cancer. Current studies have pointed out that NF-kappa B signaling pathway is involved in castration-resistant prostate cancer [28, 29] and associated with prostate cancer cell EMT and bone metastasis [30, 31]. Here, we found CR-related gene signature was associated with these oncogenic signaling pathways, suggesting the role of circadian rhythm disruption in the development of cancer.

Moreover, the seven genes consisting of the CR-related gene signature could be potential biomarkers for the development and relapse risk of prostate cancer. FBXL22 has been identified to correlated with ER- breast cancer in an exome-wide analysis [32]; nevertheless, its specific effects on cancer, including prostate cancer, is unclear. Here, the present study indicated its role in predicting recurrence risk of prostate cancer. MTA1 has been reported to drive malignant progression and bone metastasis in prostate cancer in several studies [33, 34], whereas the protective effect of TP53 is well-known [35, 36]. RORC encoding ROR-γ could induce androgen expression in prostate cancer [37] and is a potential therapeutic target in castration-resistant prostate cancer [38]. PPARGC1A and ZFHX3 are also reported to be associated with prostate cancer [39, 40], which could partly explain their ability to predict relapse risk.

This research has several practical meanings for clinical management of PRAD patients. Firstly, we developed a new gene signature and a novel predicting nomogram that would benefit therapeutic strategies for PRAD. In addition, we identified multiple pivotal genes and signaling pathways that would help development of new therapeutic strategies for PRAD. Needless to say, this study was mainly performed using bioinformatics, and thus laboratory experiments are needed for further research.

The present study has several limitations that need attention. Firstly, the pivotal gene and signaling pathways in this study were identified using bioinformatics analysis, laboratory experiments are needed to investigate their biological rationales. In addition, the risk score was identified to predict response to immunotherapy of PRAD, warranting further clinical investigation.

Conclusively, we successfully developed a new CR-related signature and a nomogram that could predict relapse risk of PRAD patients, which could provide insights into further researches on PRAD.

Author Contributions

Jin Liu, Shijie Yang and Zhao Tan collected the data. Jin Liu, Shijie Yang, Xinda Song and Wenping Li analyzed the data and wrote the manuscript. All authors reviewed and approved the final manuscript.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Funding

The study was supported by 2020 Government-Funded Clinical Medical Talent Training Program by Hebei Provincial Finance Department.

References

1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. CA Cancer J Clin. 2021; 71:7–33. https://doi.org/10.3322/caac.21654 [PubMed]
2. Ward PS, Thompson CB. Metabolic reprogramming: a cancer hallmark even warburg did not anticipate. Cancer Cell. 2012; 21:297–308. https://doi.org/10.1016/j.ccr.2012.02.014 [PubMed]
3. Ferlay J, Colombet M, Soerjomataram I, Dyba T, Randi G, Bettio M, Gavin A, Visser O, Bray F. Cancer incidence and mortality patterns in Europe: Estimates for 40 countries and 25 major cancers in 2018. Eur J Cancer. 2018; 103:356–87. https://doi.org/10.1016/j.ejca.2018.07.005 [PubMed]
4. Matsumoto K, Niwa N, Hattori S, Takeda T, Morita S, Kosaka T, Mizuno R, Shinojima T, Kikuchi E, Asanuma H, Oya M. Establishment of the optimal follow-up schedule after radical prostatectomy. Urol Oncol. 2018; 36:341.e9–14. https://doi.org/10.1016/j.urolonc.2018.04.003 [PubMed]
5. Numan A, Singh S, Zhan Y, Li L, Khalid M, Rilla K, Ranjan S, Cinti S. Advanced nanoengineered-customized point-of-care tools for prostate-specific antigen. Mikrochim Acta. 2021; 189:27. https://doi.org/10.1007/s00604-021-05127-y [PubMed]
6. Epstein JI, Zelefsky MJ, Sjoberg DD, Nelson JB, Egevad L, Magi-Galluzzi C, Vickers AJ, Parwani AV, Reuter VE, Fine SW, Eastham JA, Wiklund P, Han M, et al. A Contemporary Prostate Cancer Grading System: A Validated Alternative to the Gleason Score. Eur Urol. 2016; 69:428–35. https://doi.org/10.1016/j.eururo.2015.06.046 [PubMed]
7. Papagiannakopoulos T, Bauer MR, Davidson SM, Heimann M, Subbaraj L, Bhutkar A, Bartlebaugh J, Vander Heiden MG, Jacks T. Circadian Rhythm Disruption Promotes Lung Tumorigenesis. Cell Metab. 2016; 24:324–31. https://doi.org/10.1016/j.cmet.2016.07.001 [PubMed]
8. Wendeu-Foyet MG, Menegaux F. Circadian Disruption and Prostate Cancer Risk: An Updated Review of Epidemiological Evidences. Cancer Epidemiol Biomarkers Prev. 2017; 26:985–91. https://doi.org/10.1158/1055-9965.EPI-16-1030 [PubMed]
9. Zhou R, Chen X, Liang J, Chen Q, Tian H, Yang C, Liu C. A circadian rhythm-related gene signature associated with tumor immunity, cisplatin efficacy, and prognosis in bladder cancer. Aging (Albany NY). 2021; 13:25153–79. https://doi.org/10.18632/aging.203733 [PubMed]
10. Li M, Chen Z, Jiang T, Yang X, Du Y, Liang J, Wang L, Xi J, Lin M, Feng M. Circadian rhythm-associated clinical relevance and Tumor Microenvironment of Non-small Cell Lung Cancer. J Cancer. 2021; 12:2582–97. https://doi.org/10.7150/jca.52454 [PubMed]
11. Aiello I, Fedele MLM, Román F, Marpegan L, Caldart C, Chiesa JJ, Golombek DA, Finkielstein CV, Paladino N. Circadian disruption promotes tumor-immune microenvironment remodeling favoring tumor cell proliferation. Sci Adv. 2020; 6:eaaz4530. https://doi.org/10.1126/sciadv.aaz4530 [PubMed]
12. Xuan W, Khan F, James CD, Heimberger AB, Lesniak MS, Chen P. Circadian regulation of cancer cell and tumor microenvironment crosstalk. Trends Cell Biol. 2021; 31:940–50. https://doi.org/10.1016/j.tcb.2021.06.008 [PubMed]
13. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14:7. https://doi.org/10.1186/1471-2105-14-7 [PubMed]
14. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–7. https://doi.org/10.1089/omi.2011.0118 [PubMed]
15. Fu M, Wang Q, Wang H, Dai Y, Wang J, Kang W, Cui Z, Jin X. Immune-Related Genes Are Prognostic Markers for Prostate Cancer Recurrence. Front Genet. 2021; 12:639642. https://doi.org/10.3389/fgene.2021.639642 [PubMed]
16. Zhang Y, Zhang R, Liang F, Zhang L, Liang X. Identification of Metabolism-Associated Prostate Cancer Subtypes and Construction of a Prognostic Risk Model. Front Oncol. 2020; 10:598801. https://doi.org/10.3389/fonc.2020.598801 [PubMed]
17. Wang L, Yao Y, Xu C, Wang X, Wu D, Hong Z. Exploration of the Tumor Mutational Burden as a Prognostic Biomarker and Related Hub Gene Identification in Prostate Cancer. Technol Cancer Res Treat. 2021; 20:15330338211052154. https://doi.org/10.1177/15330338211052154 [PubMed]
18. Dong C, Dang D, Zhao X, Wang Y, Wang Z, Zhang C. Integrative Characterization of the Role of IL27 In Melanoma Using Bioinformatics Analysis. Front Immunol. 2021; 12:713001. https://doi.org/10.3389/fimmu.2021.713001 [PubMed]
19. Goodman AM, Kato S, Bazhenova L, Patel SP, Frampton GM, Miller V, Stephens PJ, Daniels GA, Kurzrock R. Tumor Mutational Burden as an Independent Predictor of Response to Immunotherapy in Diverse Cancers. Mol Cancer Ther. 2017; 16:2598–608. https://doi.org/10.1158/1535-7163.MCT-17-0386 [PubMed]
20. Carbone DP, Reck M, Paz-Ares L, Creelan B, Horn L, Steins M, Felip E, van den Heuvel MM, Ciuleanu TE, Badin F, Ready N, Hiltermann TJN, Nair S, et al, and CheckMate 026 Investigators. First-Line Nivolumab in Stage IV or Recurrent Non-Small-Cell Lung Cancer. N Engl J Med. 2017; 376:2415–26. https://doi.org/10.1056/NEJMoa1613493 [PubMed]
21. Samstein RM, Lee CH, Shoushtari AN, Hellmann MD, Shen R, Janjigian YY, Barron DA, Zehir A, Jordan EJ, Omuro A, Kaley TJ, Kendall SM, Motzer RJ, et al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet. 2019; 51:202–6. https://doi.org/10.1038/s41588-018-0312-8 [PubMed]
22. Green AS, Chapuis N, Lacombe C, Mayeux P, Bouscary D, Tamburini J. LKB1/AMPK/mTOR signaling pathway in hematological malignancies: from metabolism to cancer cell biology. Cell Cycle. 2011; 10:2115–20. https://doi.org/10.4161/cc.10.13.16244 [PubMed]
23. Penfold L, Woods A, Muckett P, Nikitin AY, Kent TR, Zhang S, Graham R, Pollard A, Carling D. CAMKK2 Promotes Prostate Cancer Independently of AMPK via Increased Lipogenesis. Cancer Res. 2018; 78:6747–61. https://doi.org/10.1158/0008-5472.CAN-18-0585 [PubMed]
24. Yuan H, Han Y, Wang X, Li N, Liu Q, Yin Y, Wang H, Pan L, Li L, Song K, Qiu T, Pan Q, Chen Q, et al. SETD2 Restricts Prostate Cancer Metastasis by Integrating EZH2 and AMPK Signaling Pathways. Cancer Cell. 2020; 38:350–65.e7. https://doi.org/10.1016/j.ccell.2020.05.022 [PubMed]
25. Yan Y, Du C, Li G, Chen L, Yan Y, Chen G, Hu W, Chang L. CO suppresses prostate cancer cell growth by directly targeting LKB1/AMPK/mTOR pathway in vitro and in vivo. Urol Oncol. 2018; 36:312.e1–8. https://doi.org/10.1016/j.urolonc.2018.02.013 [PubMed]
26. Chen C, Wang H, Geng X, Zhang D, Zhu Z, Zhang G, Hou J. Metformin exerts anti-AR-negative prostate cancer activity via AMPK/autophagy signaling pathway. Cancer Cell Int. 2021; 21:404. https://doi.org/10.1186/s12935-021-02043-2 [PubMed]
27. Pällmann N, Deng K, Livgård M, Tesikova M, Jin Y, Frengen NS, Kahraman N, Mokhlis HM, Ozpolat B, Kildal W, Danielsen HE, Fazli L, Rennie PS, et al. Stress-Mediated Reprogramming of Prostate Cancer One-Carbon Cycle Drives Disease Progression. Cancer Res. 2021; 81:4066–78. https://doi.org/10.1158/0008-5472.CAN-20-3956 [PubMed]
28. Shang Z, Yu J, Sun L, Tian J, Zhu S, Zhang B, Dong Q, Jiang N, Flores-Morales A, Chang C, Niu Y. LncRNA PCAT1 activates AKT and NF-κB signaling in castration-resistant prostate cancer by regulating the PHLPP/FKBP51/IKKα complex. Nucleic Acids Res. 2019; 47:4211–25. https://doi.org/10.1093/nar/gkz108 [PubMed]
29. Thomas-Jardin SE, Dahl H, Nawas AF, Bautista M, Delk NA. NF-κB signaling promotes castration-resistant prostate cancer initiation and progression. Pharmacol Ther. 2020; 211:107538. https://doi.org/10.1016/j.pharmthera.2020.107538 [PubMed]
30. Ren D, Yang Q, Dai Y, Guo W, Du H, Song L, Peng X. Oncogenic miR-210-3p promotes prostate cancer cell EMT and bone metastasis via NF-κB signaling pathway. Mol Cancer. 2017; 16:117. https://doi.org/10.1186/s12943-017-0688-6 [PubMed]
31. Chen L, De Menna M, Groenewoud A, Thalmann GN, Kruithof-de Julio M, Snaar-Jagalska BE. A NF-ĸB-Activin A signaling axis enhances prostate cancer metastasis. Oncogene. 2020; 39:1634–51. https://doi.org/10.1038/s41388-019-1103-0 [PubMed]
32. Haddad SA, Ruiz-Narváez EA, Haiman CA, Sucheston-Campbell LE, Bensen JT, Zhu Q, Liu S, Yao S, Bandera EV, Rosenberg L, Olshan AF, Ambrosone CB, Palmer JR, Lunetta KL. An exome-wide analysis of low frequency and rare variants in relation to risk of breast cancer in African American Women: the AMBER Consortium. Carcinogenesis. 2016; 37:870–7. https://doi.org/10.1093/carcin/bgw067 [PubMed]
33. Kumar A, Dholakia K, Sikorska G, Martinez LA, Levenson AS. MTA1-Dependent Anticancer Activity of Gnetin C in Prostate Cancer. Nutrients. 2019; 11:2096. https://doi.org/10.3390/nu11092096 [PubMed]
34. Kumar A, Dhar S, Campanelli G, Butt NA, Schallheim JM, Gomez CR, Levenson AS. MTA1 drives malignant progression and bone metastasis in prostate cancer. Mol Oncol. 2018; 12:1596–607. https://doi.org/10.1002/1878-0261.12360 [PubMed]
35. Mu P, Zhang Z, Benelli M, Karthaus WR, Hoover E, Chen CC, Wongvipat J, Ku SY, Gao D, Cao Z, Shah N, Adams EJ, Abida W, et al. SOX2 promotes lineage plasticity and antiandrogen resistance in TP53- and RB1-deficient prostate cancer. Science. 2017; 355:84–8. https://doi.org/10.1126/science.aah4307 [PubMed]
36. Nyquist MD, Corella A, Coleman I, De Sarkar N, Kaipainen A, Ha G, Gulati R, Ang L, Chatterjee P, Lucas J, Pritchard C, Risbridger G, Isaacs J, et al. Combined TP53 and RB1 Loss Promotes Prostate Cancer Resistance to a Spectrum of Therapeutics and Confers Vulnerability to Replication Stress. Cell Rep. 2020; 31:107669. https://doi.org/10.1016/j.celrep.2020.107669 [PubMed]
37. Stone L. Prostate cancer: ROR-γ drives androgen receptor expression. Nat Rev Urol. 2016; 13:237. https://doi.org/10.1038/nrurol.2016.71 [PubMed]
38. Wang J, Zou JX, Xue X, Cai D, Zhang Y, Duan Z, Xiang Q, Yang JC, Louie MC, Borowsky AD, Gao AC, Evans CP, Lam KS, et al. ROR-γ drives androgen receptor expression and represents a therapeutic target in castration-resistant prostate cancer. Nat Med. 2016; 22:488–96. https://doi.org/10.1038/nm.4070 [PubMed]
39. Canto P, Granados JB, Feria-Bernal G, Coral-Vázquez RM, García-García E, Tejeda ME, Tapia A, Rojano-Mejía D, Méndez JP. PPARGC1A and ADIPOQ polymorphisms are associated with aggressive prostate cancer in Mexican-Mestizo men with overweight or obesity. Cancer Biomark. 2017; 19:297–303. https://doi.org/10.3233/CBM-160467 [PubMed]
40. Hu Q, Zhang B, Chen R, Fu C, A J, Fu X, Li J, Fu L, Zhang Z, Dong JT. ZFHX3 is indispensable for ERβ to inhibit cell proliferation via MYC downregulation in prostate cancer cells. Oncogenesis. 2019; 8:28. https://doi.org/10.1038/s41389-019-0138-y [PubMed]

Research Paper Volume 14, Issue 17 pp 7170—7185