Prognostic values of the core components of the mammalian circadian clock in prostate cancer

Dataset acquisition from the TCGA and GEO database

We identified suitable public datasets of PC patients in the TCGA (https://tcga-data.nci.nih.gov/tcga/) and GEO (https://www.ncbi.nlm.nih.gov/geo/) Databases. We eliminated datasets without intact gene expression data and prognostic information. Both RNA sequencing (RNA-seq) data and complete clinical annotation for each PC patient were downloaded online. In total, GSE70770 with 112 PC patients (Ross-Adams et al., 2015) and TCGA-PC dataset with 493 PC patients (Liu et al., 2018) were eventually gathered in this study for further analysis. Characteristics of the TCGA cohort and the GEO cohort were summarized in Table S1. All participants gave their informed consent for publication.

Construction of prognostic signature based on clock genes

A total of 22 CCMCCs were obtained from previously published reviews (Takahashi, 2017; Partch, Green & Takahashi, 2014). We applied the COX regression analysis to assess the effects of CCMCCs on clinical prognosis. We selected potential prognosis-related clock genes to construct the circadian clock-based risk score and the prognosis prediction model. The definition of the circadian clock-based risk score was as follows:

The circadian clock-based risk score = Σ λi

where i represents the expression of prognosis-related clock genes, λ is the coefficient that extracted from the COX regression analysis. The final formula of the circadian clock-based risk score was as follows: the circadian clock-based risk score = (0.8*expression value of CSNK1D) + (−1.824*expression value of BTRC) + (−1.7645*expression value of CLOCK) + (0.4555*expression value of CSNK1E) − (1.239*expression value of FBXL3) + (−1.56*expression value of PRKAA2) + (1.325*expression value of DBP) + (0.433*expression value of NR1D2) + (1.049*expression value of RORB). The training subset and the internal validation subset was from the TCGA cohort, while GSE70770 was used as an external validation dataset. The receiver operating characteristic (ROC) analysis was applied to validate the performance of the proposed model. By R package termed “survivalROC”, areas under the ROC curve (AUC) were calculated.

Functional enrichment analysis

Differentially expressed genes (DEGs) between the high and low circadian clock-based risk score groups were recognized by the “limma” R package. DEGs with absolute value change of expression more than 2 and p value less than 0.05 were selected for signal pathway analysis. The Gene Ontology (GO, http://www.geneontology.org/) and the Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.kegg.jp/kegg/) enrichment analyses were utilized to reveal unique biological processes and signal pathways between the high and low circadian clock-based risk score groups.

Statistical analysis

The survival analysis of each CCMCC and the circadian clock-based risk score was conducted by the Kaplan–Meier (KM) method and the log-rank test. The optimal cutoff point was evaluated by the R package “survminer”. The differences between the two groups were compared with the t-test. All statistical analyses were two-tailed. And p value less than 0.05 was considered as statistical significance. The visualization of all statistical results was performed by the SPSS 22.0 software and the R 4.0.1 software.

Expression profile of CCMCCs in PC

In the TCGA dataset, the gene expression of 22 CCMCCs and clinical features of 493 PC patients from TCGA database were summarized in Fig. 1A and Table S1. Most subjects were under the age of 65 (327/493, 66.3%). In addition, a majority of patients in the TCGA cohort were at N0 stage (342/493, 69.4%). Pathology T3-4 patients also accounted for a large proportion (300/493, 60.9%). The number of patients who received neoadjuvant therapy was also summarized. We further explored the association between CCMCCs’ expression levels and clinicopathological characteristics. Importantly, the expression level of SKP1 was downregulated in PC patients over 65 years old (Fig. 1B). For lymph node (N) stage, CRY1 (p < 0.006), CSNK1D (p < 0.0014), EBXL21 (p < 0.001), PER1 (p < 0.001), PER2 (p < 0.012), PRKAA2 (p < 0.011), and SKP1 (p < 0.007) were differentially expressed between the N0 and the N1 groups (Fig. 1C). Moreover, for tumor (T) stage, ten CCMCCs differentially expressed between the T2 group and the T3-4 group (all p < 0.05, Fig. 1D). In addition, relative high expression levels of BTRC, CRY2, FBXL21, PER1, PER2, and PRKAA2, were found in the T2N0 group when compared with the T3-4N1 group (both p < 0.05, Fig. S1). In the GEO cohort, we further conducted correlation analysis on Gleason grade, PSA level, and expression levels of 22 CCMCCs (Fig. S2). However, regardless of statistical correlation between Gleason grade and NR1D2 (p = 0.011, r = 0.241), CLOCK (p = 0.045, r = −0.192), as well as SKP1 (p = 0.320, r = −0.205), the correlation was weak (Fig. S2). No statistic relationship was found between PSA level, and expression levels of 22 CCMCCs (Fig. S2).

Relationship between CCMCCs and prognosis

To evaluate the association between CCMCCs expression and clinical outcomes, we performed the survival analysis. In the TCGA cohort, the sample size of prostate cancer patients who owned disease-free survival (DFS), progression-free survival (PFS), and overall survival (OS) data was 333, 493, and 493, respectively (Figs. 2–6). We found that high expression levels of ARNTL (p = 0.037), CLOCK (p = 0.006), PER2 (p = 0.0051), PER3 (p = 0.013), BTRC (p = 0.026), CUL1 (p = 0.0054), FBXL21 (p = 0.00035), FBXL3 (p = 0.012), PRKAA1 (p = 0.013), PRKAA2 (p = 0.015), and RORB (p = 0.0067) were related to longer DFS (Figs. 2 and 3). On the contrary, overexpression of CRY1 (p = 0.0015), CSNK1D (p = 0.00032), CSNK1E (p = 0.0019), DBP (p = 0.00075), NR1D1 (p = 0.04), and SKP1 (p = 0.017) were associated with poor DFS (Figs. 2 and 3). For PFS, KM survival analysis found that PFS was statistically related to the expression levels of 18 CCMCCs (Figs. 4 and 5). In addition, there were high connections between the expression levels of CSNK1D (p = 0.016), FBXL21 (p = 0.039), NFIL3 (p = 0.037), PER1 (p = 0.035), RORB (p = 0.037) and OS in KM curves (Fig. 6).

Considering the influence of TNM stage on prognosis, we applied survival subgroup analysis to the T2N0 cohort and the T3-4N1 cohort (Tables S2–S7). Importantly, at both T2N0 stage and T3-4N1 stage, high BTRC, CLOCK, CRY1, FBXL3, PER3, and RORA expression indicated longer DFS (both p < 0.05, Tables S2–S3). It was also worth to mention that high DBP linked to shorter PFS in T2N0 prostate cancer (p = 0.013, Table S4), while the contrary result was found in T3-4N1 stage patients (p = 0.0016, Table S5). For OS, high NR1D1 expression was only significantly related to better prognosis in T2N0 stage patients (p = 0.041, Table S6), while the negative result was found in T3-4N1 stage patients (p = 0.13, Table S7).

Identification of potential prognostic CCMCCs

In the univariate Cox regression analysis, BTRC (p = 0.020136), CLOCK (p = 0.037606), CSNK1D (p = 0.002343), CSNK1E (p = 0.006002), FBXL3 (p = 0.01807), and PRKAA2 (p = 0.028131) were significantly associated with DFS (Table 1). Additionally, among 22 CCMCCs, six genes, including CLOCK (p = 0.042614), CSNK1E (p = 0.029721), DBP (p = 0.002931), NR1D2 (p = 0.041884), PRKAA2 (p = 0.018837), and RORB (p = 0.00098), showed significant relationship with PFS (Table 2). For OS, CSNK1D was found to be potential prognostic factor (p = 0.000651, Table 3). After combining above findings and removing repetitive gene, 9 CCMCCs were left, including CSNK1D, BTRC, CLOCK, CSNK1E, FBXL3, PRKAA2, DBP, NR1D2, and RORB. We proposed a hypothesis that 9 CCMCCs were key prognostic genes in PC and incorporating them could effectively predict the prognosis of PC patients.

Construction and validation of circadian clock-based risk score

To verify our hypothesis, we separated the TCGA cohort into the test set and the internal validation set in a ratio of 7:3. All enrolled patients in the cohort underwent surgical treatment, thus we mainly investigated the application of the predictive model in DFS prediction. As shown in Figs. 7A–7B, the AUC value of the predictive model in the training cohort and the validation cohort for 3-year DFS was 0.742 and 0.821, respectively. The ROC curves and AUC values (0.728 for the test cohort; 0.753 for the validation cohort) further indicated the satisfactory predictive power of the circadian clock-based signature in DFS prediction (Figs. 7C–7D). PC patients with high circadian clock-based risk score also had lower DFS than the low risk score group (p < 0.0001, Fig. 7E). In an independent cohort with 112 PC patients who prostatectomy, we divided the GSE70770 cohort into high and low risk score groups and conducted the KM survival analysis. Results showed that the high circadian clock-based risk score was correlated with early relapse (vs low circadian clock-based risk score: 13.000 months vs 21.000 months, p = 0.06; Fig. S3).

Moreover, we explored the impact of TNM stage. When compared with clinical prognostic factors, such as T (AUC values range from 0.502 to 0.808) and N stage (AUC values range from 0.500 to 0.515), the predictive performance of the proposed model was superior in DFS (Fig. S4). Similar results were also detected in PFS and OS prediction (Fig. S4). The improvement of utility when combined circadian clock-based risk score model with T and N stage was limited. Then, for better clinical application, we investigated the power of circadian clock-based risk score model in T2N0 disease and T3-4N1 disease (Fig. S5). AUC values of DFS curves indicated that the risk model performed better in patients with T2N0 stage (vs T3-4N1 stage, AUC value: 0.749–0.834 vs 0.515–0.745, Fig. S5). For PFS, the predictive values of risk score model in T2N0 disease and T3-4N1 disease was almost (Fig. S5). In T2N0 stage, high risk score was significantly related to shorter DFS (p = 0.0024) and PFS (p = 0.034), while no statistical significance was found in OS (p = 0.22, Fig. S6). In T3-4N0 stage, high circadian clock-based risk score was significantly related to shorter PFS (p = 0.016), while no statistical significance was found in DFS (p = 0.12) and OS (p = 0.24, Fig. S6).

We also tentatively applied the risk score model to PFS and OS prediction (Figs. S7–S8). The 3-year and 5-year AUC values of PFS curves ranged from 0.607–0.735 (Fig. S7), while the AUC values of OS curves were higher than 0.700 (Fig. S8). The high circadian clock-based risk score indicated poor PFS (p < 0.0001), and higher 5-year death rate (p = 0.007).

Correlation between clinicopathological parameters and circadian clock-based risk score

The relationship between clinicopathological parameters and circadian clock-based risk scores was evaluated in the TCGA cohort. There was no significant difference in the risk score according to age (p = 0.19; Fig. S9A). However, high risk score was significantly linked to high T status (p = 0.00015) and N status (p = 0.00051; Figs. S9B–S9C). When in comparison with T2N0 stage, we found that higher circadian clock-based risk score was found in the T3-4N1 stage (the TCGA cohort, p = 4e−05; the GEO cohort, p = 3.3e−06; Figs. S9D–S9E).

Functional analysis of circadian clock-based risk score

In the TCGA cohort, on the basis of the circadian clock-based risk score, 246 PC patients were assigned to the high risk score group, while 247 PC patients were assigned to the low risk score group. There were a total of 1114 DEGs between the two groups, including 4643 upregulated and 6471 downregulated DEGs (Fig. 8A). According to the GO analysis results, the top 3 enriched biological processes were ribosome biogenesis (GO:0042254), viral gene expression (GO:0019080), and viral transcription (GO:0019083; Fig. 8B). For molecular function, the top 3 enriched GO terms were small GTPase binding (GO:0031267), Ras GTPase binding (GO:0017016), and structural constituent of ribosome (GO:0003735; Fig. 8B). Moreover, Ribosome (hsa03010, p = 5.15E−15) was also one of the top 10 circadian clock-related pathways in PC (Fig. 8C).

Correlation between CCMCCs and several key prognostic genes

By literature consulting, we noticed some key prognostic genes, such as PTEN, TP53, BRCA1, BRCA2, ATM, RB1, PALB2, CHEK2, MLH1, MSH2, MSH6, and PMS2 (Rebello et al., 2021; Muñoz Fontela et al., 2016; Abida et al., 2019). Then, we explored the correlation between CCMCCs and the mentioned genes (Table S8 and Fig. S10). We firstly evaluated the expression levels of CCMCCs between different mutation status of key prognostic genes in PC (Table S8). Importantly, in the TCGA cohort, the mutation rate of PTEN, TP53, BRCA1, BRCA2, ATM, RB1, PALB2, CHEK2, MLH1, MSH2, MSH6, and PMS2, was 21.26%, 15.99%, 2.23%, 5.06%, 6.07%, 9.72%, 1.62%, 1.62%, 0.81%, 0.61%, 0.81%, and 1.01%, respectively. When compared to the PTEN wild-type group, higher expression of NR1D2 (p = 0.0008006), PRKAA1 (p = 0.00005925), and RORB (p = 0.00000054), were found in the PTEN mutation group (Table S8). When compared to the TP53 wild-type group, higher expression of CRY1 (p = 0.003774), CSNK1E (p = 0.0004053), and RORB (p = 0.006969), were found in the TP53 mutation group (Table S8). In PC, ATM mutation was correlated with high expression levels of CLOCK (p = 0.009059). Mutation status of other genes were also related to CCMCCs’ expression (Table S8).

Moreover, we investigated the correlation between the expression levels of CCMCCs and the above important genes. As shown in Fig. S10, a wild correlation was found among these genes. It is worth mentioning that ATM expression had the highly positive relevance to CLOCK expression (r = 0.633, p < 0.001). A highly positive correlation was also existed between RB1 and FBXL3 expression (r = 0.598, p < 0.001, Fig. S10).