Timing gone awry: distinct tumour suppressive and oncogenic roles of the circadian clock and crosstalk with hypoxia signalling in diverse malignancies

Datasets generated by The Cancer Genome Atlas were downloaded from Broad Institute GDAC Firehose [28]. Genomic, transcriptomic and clinical profiles of 21 cancer types and their non-tumour counterparts were downloaded (Additional file 2).

Firehose Level 4 copy number variation datasets were downloaded. GISTIC gene-level tables provided discrete amplification and deletion indicators [29]. Samples with ‘deep amplifications’ were identified as those with values higher than the maximum copy-ratio for each chromosome arm (+ 2). Samples with ‘deep deletions’ were identified as those with values lower than the minimum copy-ratio for each chromosome arm (− 2). ‘Shallow amplifications’ and ‘shallow deletions’ were identified from samples with GISTIC indicators of + 1 and − 1 respectively.

Recurrently deleted/amplified genes were identified from genes that were deleted/amplified in at least 20% of samples within a cancer type and at least one-third of cancers (> seven cancers). Putative loss-of-function genes (Clock^Loss) were defined as genes that were recurrently deleted and downregulated in tumour versus non-tumour samples. Putative gain-of-function genes (Clock^Gain) were defined as genes that were recurrently amplified and upregulated in tumour versus non-tumour samples. Clock^Loss genes were CLOCK, CRY2, FBXL3, FBXW11, NR1D2, PER1, PER2, PER3, PRKAA2, RORA and RORB. Clock^gain genes were ARNTL2 and NR1D1. For each patient, Clock^Loss and Clock^Gain scores were calculated from the mean log₂ expression values of genes within each set. Molecular assessment of tumour hypoxia was performed using a 52-hypoxia gene signature [30]. Hypoxia scores were determined from the mean log₂ expression of the 52 genes.

We have published detailed methods for the aforementioned analyses [31], hence, the methods will not be repeated here. Briefly, for survival analyses, patients were separated into survival quartiles based on their Clock^Loss and Clock^Gain scores. Cox proportional hazards regression, Kaplan–Meier and receiver operating characteristic analyses were performed using the R survcomp, survival and survminer packages according to previous methods. For analyses in Fig. 5 investigating the crosstalk between hypoxia and the circadian clock, patients were separated into four groups based on their median clock and hypoxia scores. Nonparametric Spearman’s rank-order correlation analyses were performed to determine the relationship between clock and hypoxia scores in Fig. 5. Circular heatmaps in Fig. 6 were generated from Clock^Loss scores and HIF-target genes (CA9, VEGFA and LDHA) log₂ expression values in glioma patients. Clock^Loss scores were ranked from high (purple) to low (yellow) in the heatmap. HIF-target genes were ranked by decreasing order of Clock^Loss scores. Spearman’s correlation analyses were performed on Clock^Loss scores and HIF-target genes expression values. Differential expression analyses between tumour and non-tumour samples and between the 4th and 1st quartile patients determined using Clock^Loss and Clock^Gain scores were performed using the Bayes method and linear model followed by Benjamini–Hochberg procedure for adjusting false discovery rates. Multidimensional scaling analyses in Fig. 2 were performed using the R vegan package (Euclidean distance) and permutational multivariate analysis of variance (PERMANOVA) was used to determine statistical difference between tumour and non-tumour samples.

Differentially expressed genes (DEGs) were fed into GeneCodis [32] and Enrichr [33, 34]. Using GeneCodis, DEGs were mapped against the KEGG and Gene Ontology databases. To determine whether DEGs were enriched for targets of stem cell-related transcription factors, Enrichr was used to map DEGs to ENCODE and ChEA chromatin immunoprecipitation sequencing profiles.

All plots were generated using R pheatmap and ggplot2 packages.

We retrieved 32 genes representing core circadian clock components from Kyoto Encyclopaedia of Genes and Genomes (KEGG) (Additional file 1). To determine the extent of circadian dysregulation across diverse malignancies, we analysed tumour copy number and mRNA differential expression profiles (tumour versus non-tumour) of 18,484 samples across 21 cancer types (Additional file 1) (Fig. 1a). Chromophobe renal cell carcinoma (KICH) exhibited the highest fraction of samples harbouring deleted clock genes (Fig. 1b). This was in contrast with another kidney cancer subtype, papillary renal cell carcinoma (KIRP), which had the lowest frequency of somatic deletions (Fig. 1b). When considering somatic amplifications, we observed that this was the highest in lung squamous cell carcinoma (LUSC) and the lowest in glioma (GBMLGG) (Fig. 1b).

At an individual gene level, somatic deletions were observed in 19 circadian clock genes in at least 20% of samples within a cancer type and at least one-third of cancer types (> seven cancers) (Fig. 1a, c). Global somatic losses were observed in core circadian pacemaker genes, PER1, PER2, PER3, ARNTL (BMAL1), CLOCK, NR1D2 (REV-ERBB), RORA and RORB (Fig. 1c). For instance, PER1 was deleted in 16 cancers, while PER3 and CLOCK were lost in 15 and 12 cancers respectively (Fig. 1c). On the other hand, a distinct set of 12 genes exhibited global patterns of somatic gains, which included core clock genes namely ARNTL2 (BMAL2), CRY1, NR1D1 (REV-ERBA) and RORC (Fig. 1d). ARNTL2 was one of the most amplified genes found in 16 cancers, followed by RORC in 14 cancers and NR1D1 and CRY1 in 12 cancers each (Fig. 1d).

Somatic copy number alterations (SCNAs) associated with differential transcript expression may represent putative loss- or gain-of-function events. Somatic losses accompanying transcript downregulation in tumours could indicate a loss-of-function and vice versa, somatic gains linked to transcript upregulation could imply a gain-of-function. Differential expression analyses between tumour and non-tumour samples were performed to identify genes that were significantly altered in tumours (> 1.5-fold-change, P < 0.05) of at least seven cancer types (Fig. 1c, d). Of genes that exhibited SCNAs, 11 and two genes were associated with putative loss-of-function and gain-of-function phenotypes respectively (Fig. 1a, c, d).

Given their global patterns spanning multiple cancer types, we reason that putative loss- and gain-of-function phenotypes would impact patient prognosis. We hypothesize that loss-of-function genes could have tumour suppressive qualities where gene deletions may give rise to cancer. On the contrary, gain-of-function genes are likely to harbour tumour promoting properties. If this is true, patients with low expression of loss-of-function genes (Clock^Loss) would have poorer clinical outcomes. Likewise, high expression of gain-of-function genes (Clock^Gain) would be associated with more advanced disease states and poorer outcomes. To evaluate the effects of Clock^Loss and Clock^Gain on overall survival, each patient was assigned a score based on the average expression values of 11 and two genes respectively: Clock^Loss genes (CLOCK, CRY2, FBXL3, FBXW11, NR1D2, PER1, PER2, PER3, PRKAA2, RORA and RORB); Clock^Gain genes (ARNTL2 and NR1D1). For Kaplan–Meier analyses, patients were separated into survival quartiles based on their Clock^Loss and Clock^Gain scores. Intriguingly, we found that both gene sets conferred prognostic information in seven cancer cohorts, allowing the stratification of patients into risk groups based on circadian dysregulation (Fig. 2a, b). Clock^Loss was prognostic in five cohorts: bladder (P = 0.027), glioma (P < 0.0001), pan-kidney (consisting of chromophobe renal cell, clear cell renal cell and papillary renal cell carcinoma; P = 0.011), clear cell renal cell (P < 0.0001) and stomach (P = 0.0007) (Fig. 2a). Likewise, Clock^Gain was also prognostic in glioma (P < 0.0001), pan-kidney (P = 0.0034) and clear cell renal cell (P = 0.014) cohorts and two additional cohorts; lung (P = 0.046) and pancreas (P = 0.0059) (Fig. 2b).

Interestingly, we observed that patients within the 4th quartile (highest Clock^Loss scores) had the lowest mortality rates in glioma (hazard ratio [HR] = 0.188, P < 0.0001), pan-kidney (HR = 0.520, P = 0.001) and clear cell renal cell (HR = 0.292, P < 0.0001) cohorts (Table 1). This supports our initial hypothesis that downregulation/loss-of-function of putative tumour suppressive clock genes were linked to adverse clinical outcomes while patients with high expression of these genes would perform better. However, this was not the case for bladder (HR = 2.081, P = 0.0093) and stomach (HR = 2.155, P = 0.0054) cancers, where patients with high Clock^Loss scores (4th quartile) had increased death risks, suggesting that the function of Clock^Loss is tumour-type dependent (Table 1). In terms of Clock^Gain, high expression levels were consistently associated with increased mortality rates in all five cohorts, supporting the hypothesis on tumour-promoting effects of these genes: glioma (HR = 3.961, P < 0.0001), pan-kidney (HR = 1.890, P = 0.00066), clear cell renal cell (HR = 1.755, P = 0.0062), lung (HR = 2.023, P = 0.006) and pancreas (HR = 3.034, P = 0.0022) (Table 1).

Given the prognostic significance of Clock^Loss and Clock^Gain, we predict that their expression profiles would differ between tumour and non-tumour samples in these cancers. Indeed, as confirmed by multidimensional scaling analyses using Clock^Loss, there was a clear separation between tumour and non-tumour samples in all five cohorts, suggesting that circadian dysregulation is a hallmark of cancerous cells (Fig. 2c). Since Clock^Gain only involved two genes, we employed the Mann–Whitney–Wilcoxon test to compare the distribution of Clock^Gain scores in tumour and non-tumour samples. Clock^Gain was significantly upregulated in pan-kidney (P < 0.00001), clear cell renal cell (P < 0.00001) and lung cancer cohorts (P < 0.00001) (Fig. 2d). Glioma and pancreatic cancer cohorts had limited number of non-tumour samples, five and four samples respectively. Because of this, we did not observe any significant upregulation of Clock^Gain in these tumours (Fig. 2d).

Multivariate Cox proportional hazards regression was used to determine whether Clock^Loss and Clock^Gain were independent of other clinicopathological features. Despite accounting for tumour, node and metastasis (TNM) staging, both gene sets remained independent predictors of overall survival; Clock^Loss: bladder (HR = 1.776, P = 0.043), pan-kidney (HR = 0.569, P = 0.0055), clear cell renal cell (HR = 0.433, P = 0.00085), stomach (HR = 2.070, P = 0.0084) and Clock^Gain: pan-kidney (HR = 1.941, P = 0.00054), clear cell renal cell (HR = 1.856, P = 0.0032), lung (HR = 1.832, P = 0.018) and pancreas (HR = 2.890, P = 0.0042) (Table 1). The glioma cohort consisted of low- and high-grade subtypes. Clock^Gain remained a prognostic factor in histological subtypes of astrocytoma (HR = 3.048, P = 0.0018) and oligodendroglioma (HR = 2.764, P = 0.018) (Table 1). Since both gene sets were independent of tumour stage, we evaluated their ability to improve the resolution of TNM staging. Kaplan–Meier analyses revealed that further delineation of risk groups within similarly staged tumours is afforded by both gene sets; Clock^Loss: bladder (P < 0.0001), pan-kidney (P < 0.0001), clear cell renal cell (P < 0.0001), stomach (P = 0.024) (Fig. 3a) and Clock^Gain: pan-kidney (P < 0.0001), clear cell renal cell (P < 0.0001), lung (P < 0.0001), pancreas (P = 0.048), astrocytoma (P < 0.0001) and oligodendroglioma (P = 0.023) (Fig. 3b).

Receiver operating characteristic (ROC) analyses were employed to determine the predictive performance of Clock^Loss and Clock^Gain in comparison to TNM staging. Circadian gene sets were superior to TNM staging in predicting 5-year overall survival; Clock^Loss: bladder (area under the curve [AUC] = 0.639 vs. 0.626), stomach (AUC = 0.697 vs. 0.561) and Clock^Gain: pancreas (AUC = 0.757 vs. 0.593) (Fig. 3c, d). Importantly, when clock gene sets and TNM staging were considered as a combined model, predictive performance was greater than clock genes or TNM when measured separately; Clock^Loss: bladder (AUC = 0.688), pan-kidney (AUC = 0.794), clear cell renal cell (AUC = 0.801), stomach (AUC = 0.703) and Clock^Gain: pan-kidney (AUC = 0.791), clear cell renal cell (AUC = 0.775), lung (AUC = 0.677) and pancreas (AUC = 0.759) (Fig. 3c, d). Within the glioma cohort, AUCs for Clock^Loss and Clock^Gain were 0.844 and 0.727 respectively (Fig. 3c, d). Clock^Gain was a prognostic indicator in glioma subtypes and ROC analyses confirmed its predictive performance in astrocytoma (AUC = 0.727) and oligodendroglioma (AUC = 0.670) (Fig. 3d).

To further investigate the underpinning biological consequences of circadian clock dysregulation and determine how they link to unfavourable patient outcomes, we performed differential expression analyses on all transcripts to determine genes that were altered as a result of circadian perturbation by comparing patients from the 4th survival quartile to those from the 1st quartile. Interestingly, patients stratified using Clock^Loss had significantly higher number of differentially expressed genes (DEGs; − 1.5 > log₂ fold change > 1.5; P < 0.01) (Fisher’s exact test, P = 2.2e⁻¹⁶) compared to Clock^Gain (Fig. 4a) (Additional file 3). Many DEGs were found to be in common between cancer types, more so for Clock^Loss, suggesting the existence of conserved signalling cascades associated with circadian dysregulation in driving disease pathogenesis (Fig. 4a; Additional file 3). Gene Ontology (GO) and KEGG functional enrichment analyses revealed that circadian reprogramming of tumours resulted in the activation of a myriad of oncogenic pathways. Signalling pathways associated with cancer stem cell function (MAPK, Wnt, JAK-STAT [35], TGF-β and PPAR [36]), metabolism, cell adhesion, cell proliferation, cell death, transmembrane transport and extracellular matrix organisation were among some of the most altered biological processes that likely underpin tumour aggression and decreased survival in these patients (Fig. 4b, c). Interestingly, with the exception of MAPK signalling, pathways associated with cancer stem cell function were only enriched in patients stratified using Clock^Loss, suggesting that Clock^Loss genes are keenly linked to stem cell homeostasis (Fig. 4c). To further understand how DEGs were regulated, we analysed transcription factor (TF) binding using Enrichr and observed that DEGs were enriched for targets of TFs associated with self-renewal and cancer stem cell function (SUZ12, SOX2, REST, EZH2, SMAD4 and NANOG) (Fig. 4d). These TFs have previously been shown to promote metastasis, disease aggression and cancer stem cell maintenance [37–39].

Circadian oscillations of physiological processes such as metabolism, temperature and cortisol levels are affected by oxygen levels [40–42]. Hypoxic responses are gated by the circadian clock and at the genomic level, BMAL1 and HIF-1A synergistically interact to coregulate downstream genes [27]. Crosstalk between hypoxia and the clock has profound implications on cancer pathophysiology [43, 44]. We reason that tumour hypoxia could synergise with the circadian clock to impact disease progression. To determine the extent of the hypoxia-clock crosstalk, a hypoxia gene signature consisting of 52 genes was employed to calculate hypoxia scores in each patient [30]. Remarkably, we observed significant negative correlations between hypoxia and Clock^Loss scores in glioma (rho = − 0.57, P < 0.0001), pan-kidney (rho = − 0.13, P < 0.0001) and clear cell renal cell (rho = − 0.34, P < 0.0001) cohorts (Fig. 5a). Patients were separated into four groups based on their hypoxia and Clock^Loss scores for survival analyses. Kaplan–Meier analyses employing the joint hypoxia-Clock^Loss model revealed that patients with more hypoxic tumours who concurrently had lower levels of Clock^Loss genes performed the worst: glioma (HR = 7.218, P < 0.0001), pan-kidney (HR = 2.512, P < 0.0001) and clear cell renal cell (HR = 1.893, P = 0.0054) (Fig. 5b) (Table 2). This observation is consistent with the predicted tumour suppressive roles of Clock^Loss.

The trend is flipped when considering Clock^Gain; significant positive correlations between hypoxia and Clock^Gain scores were observed in glioma (rho = 0.35, P < 0.0001), pan-kidney (rho = 0.17, P < 0.0001), clear cell renal cell (rho = 0.13, P = 0.0027), lung (rho = 0.50, P < 0.0001) and pancreas (rho = 0.66, P < 0.0001) (Fig. 5c). Consistently, patients with more hypoxic tumours and high Clock^Gain scores had the poorest survival outcomes: glioma (HR = 9.210, P < 0.0001), astrocytoma (HR = 5.684, P < 0.0001), oligodendroglioma (HR = 4.085, P = 0.0022), pan-kidney (HR = 3.079, P < 0.0001), clear cell renal cell (HR = 1.877, P = 0.0041), lung (HR = 2.037, P = 0.0012) and pancreas (HR = 1.976, P = 0.012) (Fig. 5d) (Table 2). Loss of tumour suppression or increase in tumour promoting properties resulted from circadian dysregulation is further exacerbated by hypoxia. The clock-hypoxia model may be used for delineation of patients into additional risk groups to support adjuvant therapy with hypoxia-reducing drugs in combination with mainstream chemotherapy and radiotherapy.

Not applicable.

The authors declare that they have no competing interests.

The datasets supporting the conclusions of this article are included within the article and its additional files.

Not applicable.

Not applicable.

None.

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

The datasets supporting the conclusions of this article are included within the article and its additional files.