Transcription of Clock Genes in Medulloblastoma

1. Introduction

The 2017 Nobel Prize in Physiology and Medicine was awarded to Hall, Rosbash and Young for their work on clock genes and proteins and led to the transcription translation feedback model as essential to account for the circadian oscillation of clock components [1]. Circadian rhythms in transcription of clock genes occurs in the suprachiasmatic nucleus (SCN) of the hypothalamus as well as in peripheral tissues [2,3]. Most peripheral tissues express clock genes [4]. Expression of clock genes in some peripheral tissues may be synchronized by the SCN [5], but may be independent of the SCN in others. Seven of the core clock genes have been reported to undergo a circadian oscillation in all major tissues tested including liver, kidney, lung, adipose tissue, muscle, heart, and brain [6]. Disruption of circadian rhythms can result in a variety of diseases. There are reports that disruption of circadian rhythms may be involved in the development of several cancers including lung, breast, liver, and prostate cancer [7,8,9,10,11,12]. Disruption of the circadian clock is regarded as a risk factor for cancer [13,14], and circadian clock proteins have been proposed as therapeutic targets in cancer [9,15].

The core clock genes include genes for transcription factors CLOCK, BMAL1 (aka ARNTL), BMAL2 (aka ARNTL2), NPAS2, the transcription repressor genes PER1, PER2, PER3 (PERIOD 1, 2, 3), and the cryptochrome genes CRY1 and CRY2. [5,16]. The BMAL1/CLOCK complex stimulates the transcription of the period and cryptochrome genes, by binding to E-Box elements in their promoters. The PER and CRY proteins in turn bind to CLOCK/BMAL1 forming a PER/CRY transcription/translation feedback loop (Figure 1) which drives the circadian clock [5].

In addition to the primary PER/CRY feedback loop, a secondary transcriptional/translational feedback loop regulating the core clock genes, the REV-ERB and ROR (Related Orphan Receptor) feedback loop, has been described and is regarded as part of the transcriptional structure of the circadian clock [5]. REV-ERB-alpha and REV-ERB-beta proteins are encoded for by the genes NR1D1 and NR1D2. NR1D1 is transcriptional repressor and nuclear heme receptor that has been reported to coordinate circadian rhythms of activity [17,18].

Finally, a third feedback loop was described by Cox and Takahashi [19] as a DBP/NFIL3 (D-Box binding transcription factor/Nuclear factor IL3 binding protein) loop. In the model of Cox and Takahashi DBP and NFIL3 dimerized proteins bind to D-box elements of the PERIOD, RORA/B and NR1D1/2 genes to regulate transcription [19].

Several components of the ubiquitin proteasome system have been found to regulate the expression of circadian clock genes [5], and are included in our analysis of MB clock gene expression. The best documented are the F-box ubiquitin ligase adaptors FBXL3 [20,21], FBXL21 [22,23,24], [25,26,27], FBXW7 [16,28], and the deubiquitinase USP2 [6,29,30,31,32,33]. FBXL21 is highly expressed in the suprachiasmatic nucleus, the hypothalamic nucleus which is the major pacemaker in mammals whereas FBXL3 and USP2 are expressed in many tissues [22]. A mutation in FBXL3 referred to as OVERTIME led to the discovery that the FBXL3 protein is required for the reactivation of CLOCK and BMAL1 by facilitating the degradation of CRY1 and CRY2 proteins and resulting in the loss of negative feedback [34]. Another mutation of FBXL3, AFTERHOURS, slows the degradation of CRY protein [35]. The FBXL21 protein binds to the CRY proteins and antagonizes the action of FBXL3 on the degradation of these proteins [23].

Additional clock related genes include TIMELESS, and CIPC (Clock Interacting Pacemaker), and the clock regulator kinases CSNK1δ and CSNK1ε.

Current availability of public datasets on genome wide gene transcription make it possible to investigate the expression of circadian clock genes in various cancers and in their subgroups. Medulloblastoma (MB) datasets provide an interesting source in which to study the interaction of circadian rhythm genes and cancer since Group 3 MB tumors were initially delineated, in part, by over-expression of photoreceptor and phototransduction genes [36,37,38,39,40].

In the current study we investigate the expression of clock genes, in public datasets of medulloblastoma. These include the Cavalli dataset [40] and the Swartling dataset, a reanalysis of 23 MB (including the Cavalli dataset) and non-tumor datasets, for batch variation [41]. The data enabled us to address the question of whether circadian clock genes were expressed in medulloblastoma tissue samples and if so, whether transcription of these genes was differentially expressed among the 4 consensus subgroups of medulloblastoma, Group 3, Group 4, SHH and WNT [36]. The Cavalli data allowed for further study of clock genes by the twelve subtypes within the four subgroups in this dataset; it also facilitated the examination of gene expression based on chromosome arm copy number variation. The Swartling data, moreover, enabled us to determine differences in expression of clock genes between MB tumor tissue and non-tumor (NT) cerebellum.

The Cavalli dataset also had available survival data, which allowed us to determine whether the expression of any of the core clock genes in MB tissues were associated with patient survival. The data facilitated the identification of potential therapeutic targets among the clock genes that were associated with survival.

Reports relating clock genes and melatonin synthesis have been reported [42,43]. Therefore we also examined the medulloblastoma datasets to determine whether expression of the gene for the rate-limiting step in melatonin synthesis, AANAT, was related to MB subgroups.

Several clock genes, as well as AANAT, are located on chromosome 17, the chromosome associated with copy number gains of 17q in MB and the chromosomal abnormality isochromosome 17q. The genes located on chromosome 17q include the clock gene NR1D1, the clock related gene CSNK1D as well as the gene encoding AANAT. The availability of copy number variation data in supplemental datasets of Cavalli et al. [40] facilitated the determination of whether the expression of clock genes in MB were related to copy number gain of chromosome 17q.

Finally, based on the evidence that several clock proteins act as transcription factors, or transcription regulators, we examined the most significant KEGG pathways associated with genes whose transcription was correlated with that of each of the clock genes.

2. Methods

Data Sources

Gene expression was mined from publicly available datasets of gene expression in MB, primarily the Cavalli dataset (n = 763) and the Swartling dataset (n = 1350 MB samples and 291 NT brain samples), a reanalysis of several datasets (including the Cavalli data) for batch variation. These datasets were made available through the R2 Genomics Analysis and Visualization Platform (R2 Genomics) (https://r2.amc.nl). Each of the datasets provided differential expression of genes for the four consensus subgroups of MB, the WNT, SHH, Group 3, and Group 4. The Cavalli dataset (n = 763) provided further information on gene expression in subtypes of each group, for a total of 12 subtypes, and also provided survival data. The Swartling dataset (n = 1641), a batch-corrected meta-analysis, in addition to providing gene expression of the four subgroups of MB, also included gene expression data of non-tumor tissue derived from 291 normal cerebellum tissue samples. Subsequently we refer to this group as the NT group.

Additional MB datasets found the R2 Genomics Platform included the Pfister dataset, the Northcott MAGIC MB dataset, and the Gilbertson dataset. Additional NT datasets included the NT brain dataset of Berchtold, and the NT cerebellum dataset of Roth.

Data Analysis

The R2 Genomics platform was also used for basic statistical analysis, including analysis of variance (anova), of differentially expressed clock genes and, in the Cavalli dataset, for Kaplan-Meier analysis of survival data related to these genes. In addition, Cox proportional hazards regression analysis was performed through the R2 Genomics platform for the clock genes listed in Table 1. The Morpheus program (Broad Institute) was used to produce a heatmap and cluster analysis of clock gene expression in the Cavalli dataset (Figure 2).

Chromosomal copy number variation data, found in the supplemental tables of Cavalli et al. [40], facilitated the examination of clock gene expression as it relates to copy number gain of Chromosome 17q.

The KEGG pathways associated with groups of genes correlated with each of the clock genes were determined using the gene set analysis application in the R2 Genomics platform.

3. Results

3.1. Clock Gene Expression in Medulloblastoma Subgroups

Our list of clock genes included the core clock genes as well as genes that encode proteins that regulate clock gene expression. Differential expression of various clock genes between MB subgroups was highly significant by analysis of variance (anova). Table 1 shows the significance of MB group differences in the Cavalli dataset and the relative significance in the Swartling meta-analysis, which included a NT group in addition to the 4 MB subgroups.

Our gene expression results are discussed in the order of clock proteins presented in the 2021 review of Cox and Takahashi [19]: genes of the core clock proteins CLOCK and BMAL1, genes of the PERIOD/CRYPTOCHROME feedback loop (Figure 1), genes of the REV-ERB and ROR feedback loop, and genes of a third feedback loop described as a DBP/NFIL3 loop. This will be followed by presentation of gene expression associated with the ubiquitin proteasome system, which has been shown to have a key role in circadian regulation of clock proteins; these genes include FBXL3, FBXL21, FBXW7, BTRCP, and USP2. Our focus is on the results with the greatest statistical significance, and on results that are confirmed across datasets. Differential expression of clock genes by MB subgroup in the Cavalli dataset are illustrated in the heatmap and cluster analysis of Figure 2.

Table 2 shows the statistics for the Kaplan-Meier survival curve scans (chi squared and p values) up to 144 months, as well as the Hazards Ratios determined using the Cox proportional HR analysis. The most significant Kaplan-Meier scans (high vs low gene expression) were found for the following clock related genes: USP2, CRY1, and CLOCK (Table 2).

3.1.1. CLOCK and BMAL Expression in MB Subgroups

Differential expression of CLOCK among the four major medulloblastoma (MB) subgroups was observed across multiple datasets. Analysis of the Cavalli dataset revealed significant differences in CLOCK expression (F = 77.07, p = 1.59 x 10⁻⁴³), which were validated in the Pfister (F = 21.78, p = 2.21 x 10⁻¹²), Northcott (F = 52.76, p = 2.44 x 10⁻²⁰), and recalculated Swartling datasets (F = 155.46, p = 6.27 x 10⁻⁸⁶). The Swartling dataset, which included a normal tissue (NT) group, allowed for subgroup comparisons against NT (Figure 3). CLOCK expression was lowest in Group 3 MB, with a highly significant reduction compared to NT (p = 1.25 x 10⁻³⁹). In contrast, no significant difference in CLOCK expression was found between Group 4 and NT.

The CLOCK protein, together with its binding partner BMAL1, acts as a transcription factor regulating genes such as PERIOD, CRYPTOCHROME, and NR1D1/NR1D2. Low CLOCK expression, as observed in Group 3 MB, may disrupt cellular responses to DNA damage, as CLOCK has been implicated in this process. Notably, low CLOCK expression was associated with worse overall survival (Figure 3, Table 2). These findings suggest that the role of reduced CLOCK protein in the development and proliferation of Group 3 MB warrants further investigation.

Expression of ARNTL (coding for BMAL1) was also differentially expressed across MB subgroups in the Cavalli dataset (F = 59.29, p = 1.92 x 10⁻³⁴) and validated in other datasets, including the Swartling dataset (F = 137.94, p = 1.88 x 10⁻⁷⁷). While ARNTL expression was depressed in the SHH subgroup compared to NT (p = 2.8 x 10⁻⁸⁷), it was significantly elevated in Group 3 MB (Figure 3). A strong statistical correlation between ARNTL and RORA expression (Figure 8B) supports a potential functional relationship, such as RORA promoting BMAL1 transcription.

The CLOCK homolog, NPAS2, whose protein can dimerize with BMAL1 as an alternative to CLOCK, exhibited highly significant differential expression (F = 283.34, p = 1.59 x 10⁻¹⁸³). Elevated NPAS2 expression was observed in Group 4 MB compared to other groups and NT, suggesting it may serve as a marker for Group 4 MB. NPAS2 has been reported as a tumor suppressor involved in DNA damage repair in breast cancer cells; however, its role in MB remains unexplored.

Expression of ARNTL2 (coding for BMAL2) was lowest in the WNT subgroup (p = 1.82 x 10⁻³⁷), with reduced ARNTL2 expression associated with better survival outcomes (Chi = 12.91, p = 3.27 x 10⁻⁴) (Table 2). Interestingly, ARNTL2 was specifically elevated in Group 4α, one of the 12 Cavalli subtypes.

3.1.2. BHLHE40 and Survival in WNT MBs

BHLHE40 (aka DEC1), a transcriptional repressor that inhibits the transcriptional activity of the CLOCK/BMAL complex [44], was differentially expressed at a highly significant level (F = 209.02, p = 7.20 x 10⁻¹⁴⁴). Expression levels of BHLHE40 were elevated in the WNT subgroup compared to other subgroups and normal cerebellar tissue (NT) (Figure 3).

Our analysis revealed that this upregulation of BHLHE40 in the WNT subgroup aligns with observations in other cancers, as reported by Kiss et al., who listed various malignancies with increased BHLHE40 expression [45]. This finding adds the WNT subtype of medulloblastoma (MB) to the list of cancers where BHLHE40 is upregulated.

Kaplan-Meier survival analysis (Table 2) demonstrated that elevated BHLHE40 expression is associated with improved survival outcomes in WNT MB patients, further highlighting its potential role as a prognostic marker.

3.1.3. Cryptochrome (CRY) and Period (PER) Gene Expression

Cryptochrome proteins are required to maintain circadian rhythmicity in mammals Cryptochrome proteins are blue light photoreceptors in plants and in insects including Drosophila, but knockout experiments in mice lead to the conclusion that they do not function as direct photoreceptors in mammals [46]. Differential expression of CRY1 (F = 160.65, p = 1.29 x 10^-80) andCRY2(F = 35.25; 2.48 x 10^-21) was noted among the four major subgroups of MB in the Cavalli dataset and in the Swartling dataset. Compared to the NT group in the Swartling dataset CRY1 expression (Figure 4A) was elevated in Group 3 MB by t-test (p = 5.73 x 10^-48), but depressed in the SHH (p = 4.81 x 10^-47) and WNT (p = 7.49 x 10^-17) groups. CRY2 (Figure 4B) expression was depressed in Group 3 (p = 5.76 x 10^-18), SHH (p = 1.41 x 10^-46), and WNT (p = 3.06 x 10^-30) groups compared to the NT group.

Among the 12 MB subtypes of the 4 subgroups, highest expression of CRY1 was found in the Cavalli Group 3 alpha subtype (p = 2.77 x 10^-56 vs NT). Individuals with high levels of CRY1 were mostly less than 10 years of age (Figure 4C). High expression was associated with poor survival (Table 2).

PERIOD proteins adjust circadian rhythms to changes in photoperiod [47]. They function as transcriptional repressors in the transcription/translation feedback loop (Figure 1). Differential expression of the Period genes, PER1 and PER2, was observed between MB subgroups in the Cavalli dataset and in the Swartling dataset (Figure 5). Compared to the NT group in the Swartling dataset the major statistical differences were a depression of PER1 expression in the WNT subgroup (p = 2.34 x 10^-20) and in Group 4 (p = 5.00 x 10^-20). Decreased expression of PER2 expression was observed in the WNTgroup (p = 1.72 x 10^-32) and in Group 3 (p =5.29 x 10^-31) groups compared to the NT group.

3.1.4. TIMELESS and MB

The TIMELESS gene codes for a protein, TIM, that inhibits activation of the PER1 gene and regulates minor changes in circadian rhythms by interacting with PER2 and CRY2 proteins [48,49]. It also contributes to the DNA damage response [50]. In the Swartling dataset TIMELESS expression was substantially elevated in the MB groups G3, G4, and WNT compared to the non-tumor group (Figure 6). In the Cavalli dataset, high expression was associated with worse survival based on the Kaplan-Meier (Table 2).

3.1.5. CSNK1D and Chromosome 17q

CSNK1D, a gene located on chromosome 17q, encodes a casein kinase which phosphorylates PER1 and PER2 proteins. Phosphorylation together with ubiquitination enables the degradation of these proteins by the proteasome [51] (Figure 1). Figure 7 shows an increase in expression of CSNK1D in Groups 3 and 4 and major decrease in expression of CSNK1D in the SHH subgroup compared to the other subgroups, including the NT group. Our analysis found that high expression of this gene was statistically associated with copy number gain of chromosome 17q; its expression was elevated in individuals with 17q copy number gain.

3.1.6. REV-ERB and RORA Expression

The ROR (RORA, RORB, RORC) and NR1D (NR1D1 and NR1D2) genes encode proteins that bind to ROR-binding elements on the BMAL1 (ARNTL) gene. RORA and REV-ERBα compete for these binding sites, regulating BMAL1 transcription [19,52].

The most notable finding was the significant depression of RORA expression (Figure 8a) in the SHH group compared to the normal tissue (NT) group (t = 17.53, p = 3.00 x 10⁻⁵⁷). RORA expression was strongly correlated with BMAL1 (ARNTL) expression in the Cavalli dataset. Among 18,267 gene combinations analyzed, RORA had the most significant correlation with BMAL1 expression (r = 0.70, p = 3.19 x 10⁻¹¹¹), supporting the role of RORA in promoting BMAL1 transcription (Figure 8B) [53].

RORA also functions as a component of the second major feedback pathway regulating clock genes (Figure 8C). Its association with the phototransduction pathway was particularly striking. Correlation analysis revealed that the KEGG pathway most statistically associated with RORA expression was the phototransduction pathway, including genes encoding proteins of the phototransduction cascade [54] (Table 4) .

Seven genes of this pathway—CNGB1, GNAT1, GRK1, GUCY2D, PDE6B, PDE6G, and RCVRN—were highly correlated with RORA (r > 0.50). These genes were markedly overexpressed in the MB Group 3 alpha subtype in the Cavalli dataset. Furthermore, AANAT, the gene encoding the rate-limiting enzyme for melatonin synthesis, was also overexpressed in Group 3 MBs (Figure 17). The KEGG pathway statistically associated with AANAT and USP2 expression was represented by the same 4 phototransduction-related genes (GNAT1, GRK1, GUCY2D, and RCVRN).

These findings are consistent with prior research by Hooper et al. [55], who compared Group 3 MB expression profiles to normal molecular developmental events. They found that Group 3 MB expression profiles resembled those of rod cell precursor cells at around 15 weeks of human retinal development. Hooper et al. described a neoplastic switch from granule cell precursors to rod cell precursors in a Group 3 MB cohort, suggesting a molecular relationship between Group 3 MBs and pineal tumors.

The refined classification of Group 3 MBs into subtypes—Group 3 alpha, beta, and gamma—in the Cavalli dataset highlights increased expression of phototransduction genes in the Group 3 alpha subtype. This suggests that Group 3 alpha MBs may originate from photoreceptor precursors. Further investigation is needed to test this hypothesis and better understand the developmental origins of Group 3 alpha MBs.

Figure 8. A. Differential expression of RORA and RORB in the Swartling dataset. RORA, F = 85.84, p = 5.18 x 10-66; RORB, F = 19.47, p = 1.18 x 10-15. Red – Group 3, Green – Group 4, Grey – NA (normal) in Swartling dataset is the non-tumor group, Blue – SHH, Yellow – WNT. B. Correlation of RORA and ARNTL (BMAL1) expression in Cavalli dataset (r = 0.70, p = 1.19 x 10^-111). MB subgroups identified by color. C. RORA and REV-ERBα (NR1D1) proteins share binding sites on the BMAL1 gene. RRE – Retinoic acid response element.

Figure 8. A. Differential expression of RORA and RORB in the Swartling dataset. RORA, F = 85.84, p = 5.18 x 10-66; RORB, F = 19.47, p = 1.18 x 10-15. Red – Group 3, Green – Group 4, Grey – NA (normal) in Swartling dataset is the non-tumor group, Blue – SHH, Yellow – WNT. B. Correlation of RORA and ARNTL (BMAL1) expression in Cavalli dataset (r = 0.70, p = 1.19 x 10^-111). MB subgroups identified by color. C. RORA and REV-ERBα (NR1D1) proteins share binding sites on the BMAL1 gene. RRE – Retinoic acid response element.

3.1.7. THRA/NR1D1, NR1D2, and Their Implications in Medulloblastoma (MB)

NR1D1and NR1D2 encode the proteins REV-ERBα and REV-ERBβ, respectively, which play key roles in circadian rhythm regulation [56]. Differential expression of the THRA/NR1D1 locus (F = 423.43, p = 1.43 x 10⁻¹⁸¹) was observed across MB subgroups, with the highest expression in Group 4 MB (Figure 9). This locus includes two genes, THRA (Thyroid Hormone Receptor A) and NR1D1, which are transcribed from opposite strands of DNA on chromosome 17q21.1 [57]. This arrangement complicates the interpretation of gene expression data, as datasets like Cavalli aggregate their expression into a single locus labeled THRA/NR1D1. Meanwhile, the Swartling dataset lists only THRA, further limiting the clarity of NR1D1-specific data.

The high expression of THRA/NR1D1 in Group 4 MB aligns with the frequent occurrence of isochromosome 17q in this subgroup, particularly in the Group 4 beta subtype, where over 95% of cases exhibit this chromosomal abnormality .[40]. Statistical analyses revealed that THRA/NR1D1 expression in Group 4 MB is strongly associated with copy number gain of chromosome 17q (Table 3). Kaplan-Meier analysis showed a significant correlation between high THRA/NR1D1 expression and poor survival (p = 1.84 x 10⁻³) (Table 2), suggesting its potential as a prognostic marker.

NR1D1 (REV-ERBα) is a transcriptional regulator that inhibits BMAL1 transcription [19,58] and regulates numerous circadian-expressed genes across multiple chromosomes (Figure 8C). The CLOCK/BMAL1 protein complex activates NR1D1 transcription [59], further linking this gene to circadian rhythm pathways. High levels of NR1D1 protein (REV-ERBα) in Group 4 MB could lead to the suppression of BMAL1 activity, disrupting downstream transcriptional networks.

For NR1D2, the most significant difference by subgroup was observed as a depression in expression in the SHH subgroup compared to normal tissue (p = 3.86 x 10⁻²⁹) (Figure 8).

The high level of differential expression of the THRA/NR1D1 reporter (Hugene 11t platform probe 8007008) in Group 4 MB suggests it could serve as a biomarker for this subgroup (Table 1, Figure 9). This elevated expression likely reflects the impact of isochromosome 17q on transcription levels.

Down-regulation of REV-ERBα (NR1D1) has been observed in various cancers, as reviewed by Gomatou et al.[60]. However, there is limited literature specifically linking NR1D1 or its encoded protein REV-ERBα to MB. Our findings that high THRA/NR1D1 expression is associated with poor survival in Group 4 MB underscore the need for follow-up studies. Exploring the role of REV-ERBα as a therapeutic target may yield new insights into treatment strategies for Group 4 MB.

3.1.8. DBP (D-Box of Albumin Promoter) and NFIL3 Gene Expression

Differential expression of DBP (D-box binding protein) (F = 55.93, p = 1.14 x 10^-32) was highly significant in the Swartling dataset, with the lowest expression values observed in the SHH group (Figure 10A). Compared to the NT group in the Swartling dataset, the depression of DBP in the SHH group was particularly pronounced (t = 15.65, p = 1.54 x 10^-47). Similarly, in the Cavalli dataset, NFIL3 expression was lowest in the SHH group compared to the other three subgroups (F = 106.98, p = 9.10 x 10^-58).

DBP gene transcription is a known target of the CLOCK/BMAL1 dimer complex. Both DBP and NFIL3 proteins bind to D-box elements of promoters of selected clock genes, playing opposing roles in their regulation. Studies by Ueda et al. [61] and Keniry et al.[62], as reviewed by Cox and Takahashi [19], highlight these opposing effects on D-box elements of clock gene promoters (Figure 10B).

DBP and NFIL3 encode proteins that are components of the third clock gene feedback pathway (Figure 10) described by Cox and Takahashi [19]. Correlation analysis of DBP (Table 4) reveals that the most significant biological pathway associated with DBP gene correlates is the GABAergic pathway, which has been linked to Group 3 MB in the 2012 consensus report [36]. Elevated DBP expression is illustrated in Group 3 MB (Figure 10), suggesting that DBP may contribute to defining this subgroup.

Highest transcription of NFIL3 was observed in Groups 3 and 4 MB in the Cavalli dataset. High expression of NFIL3 was associated with poor survival (Table 2). Zeng et al. have reviewed the abnormal expression of NFIL3 in various cancers and proposed NFIL3 as a therapeutic target [63]. However, to date, neither the NFIL3 gene nor its protein has been specifically studied as a therapeutic target in MB.

3.1.9. Clock Interacting Pacemaker (CIPC) (aka KIAA1737)

CIPC expression is included in the present study due to the report that the CIPC protein binds directly to the CLOCK/BMAL1 complex [64] and due to its extremely high statistical differential expression by MB subgroups (Swartling, F =379.47, p = 1.67 x 10^-150, Cavalli; F = 587.87, p = 2.87 x 10^-240). Compared to the NT group in the Swartling dataset CIPC was elevated in Group 4 (p = 1.29 x 10^-48), but depressed in the SHH (p = 1.28 x 10^-83) and WNT (p = 1.46 x 10^-46) groups (Figure 11A). Among the clock genes, a high correlation of CIPC and FBXL21 was noted (r= 0.79, p = 7.75 x 10^-164). Figure 11B shows the correlation of CIPC and FBXL21, with subgroup expression identified by color. High expression of both genes was associated with copy number gain of chromosome 17q (Table 3).

3.2. Ubiquitin Proteasome Pathway Regulation of Clock Genes in MB

FBXL3

The differential expression of FBXL3 was significant in the Swartling dataset (F = 55.21, p = 1.31 x 10^-43 (Figure 12A). By t test the most significantly different expression of FBXL3 from the normal group was depression in the WNT group.

FBXL15

Figure 12B shows that expression of FBXL15 was depressed in the SHH MB subgroup compared to the other 3 groups as well as to the NT group.

FBXL21

Expression values of FBXL21 were not available in the Swartling dataset (due to its classification as a pseudogene in humans), preventing a comparison of expression values to the NT group. However, using the Megasample application of the R2 Genomics program, we were able to compare expression of FBXL21 in medulloblastoma (Pfister and Gilbertson datasets) to normal brain (Berchtold) and normal cerebellum (Roth dataset) using the same gene chip (Affytmetrix u133p2) and the same reporter (1555412_at). This analysis showed that FBXL21 expression was elevated in MB compared to that of the NT group (Figure 13A).

Differential expression of FBXL21 transcription was also highly significant in the Cavalli dataset (F = 289.07, p = 4.14 x 10^-125), with the mimimal expression in the SHH MB group (Figure 13B); expression of FBXL21 in Group 3 was 25-fold greater than in the SHH group, while expression in Group 4 was approximately 35-fold greater than in the SHH group. This over-expression is primarily in children less than 15 years of age (Figure 13C). The findings of elevated expression of FBXL21 in groups 3 and 4 were confirmed in the Northcott MAGIC MB dataset, in the Pfister MB dataset, and in the Gilbertson dataset, all available in the R2 Genomics platform.

FBXW7

FBXW7 expresssion was depressed in all 4 subgroups compared to the NT group in the Swartling dataset (Figure 14A). However FBXW7 expression was specifically elevated in one of the SHH MB subtypes in the Cavalli dataset, the SHH gamma subtype (Figure 14B), a subtype whose subjects are infants of less than 3 years of age. The SHH gamma subtype was the only subtype in which expression of FBXW7 was significantly elevated compared to all other subtype groups, including NT controls. Kaplan-Meier analysis showed that high expression of FBXW7 was protective (Table 2).

MYCBP2

MYCBP2 is a gene which encodes a ubiquitin ligase (aka PAM, protein associated with MYC) reportedly necessary for normal CNS development [65] and required for the degradation of the NR1D1 encoded protein REV-ERBα [66]. Expression of MYCBP2 was highest in Cavalli Group 4 and lowest in the WNT group. Differential expression of this gene was among the highest statistical significance. (Figure 15; Table 1).

USP2 – A Gene Encoding a Deubiquitinase for Clock Genes

USP2 is an enzyme which deubiquitinates various proteins, including BMAL1, CRY1 and PER1 [29,30,33], and regulates the distribution of PER1 between nucleus and cytoplasm [33]. Elevated expression of the USP2 gene was associated with Group 3 MB (Figure 16). High expression of USP2 was associated with poor survival (Chi squared = 21.73, p = 3.14 x 10^-06). By Cox proportional hazard ratio (HR) analysis of all the clock genes in Table 2, USP2 had the most significant HR (p = 3.5 x 10^-5).

3.3. AANAT Is Over-Expressed in MB Group 3

Since melatonin has been associated with circadian rhythms and clock gene expression we examined the expression of the gene encoding the rate-limiting enzyme of melatonin synthesis, AANAT (Aralkylamine N-Acetyltransferase) [67] in MB subgroups. While AANAT data were not available in the Swartling dataset, the Cavalli data show that AANAT expression was signficantly greater in Group 3 MBs than in the other groups (Fig, 17A). More specifically the elevation in AANAT expression was observed primarily in a subtype of Group 3 MBs, the Group 3 alpha subtype. AANAT expression was elevated primarily in young children (Figure 17B).

Figure 17. A. AANAT expression high in MB Group 3. F= 127.15.09, p = 9.93 x 10^-67 B. Age related expression of AANAT.

Figure 17. A. AANAT expression high in MB Group 3. F= 127.15.09, p = 9.93 x 10^-67 B. Age related expression of AANAT.

3.4. Clock Genes and Survival in the Cavalli Dataset

Several clock genes were significantly associated with survival in the Cavalli dataset. Table 2 shows the Chi squared and p values for Kaplan-Meier analysis of survival curves (high vs low scan) of clock related genes, in order of significance. Also included in Table 2 are the significant hazard ratios, as determined by the Cox proportional hazard analysis.

The three most significant Kaplan-Meier scans, were of CRY1, USP2, and CLOCK expression. Kaplan-Meier analysis showed that high expression of CRY1 and of USP2 were associated with poor survival (CRY1, Chi squared = 22.33, p = 2.29 x 10^-06, USP2, Chi squared =21.73, p = 3.14 x 10^-06) (Table 2, Figure 18); high expression of CLOCK, on the other hand, was associated with better survival (Chi squared = 20.21, p = 6.96 x 10^-06). Expression of CRY1, USP2, and CLOCK were also associated with significant Hazard ratios (Table 2).

Expression of ARNTL2 (BMAL2) and CSNK1D, was also associated with significant Kaplan-Meier scans as well as with signficant Hazard ratios (Table 2). High expression of both of these genes was associated with worse survival (and low expression with survival protection).

Expression of TIMELESS, FBXL21, and BTRC showed significantly different (high vs low expression) Kaplan-Meier scans, but no significant Hazard ratios. At a lower level of significance, high expression of RORB and of NFIL3 were associated with worse survival (Table 2). High expression of CIPC was also associated with lower survival, albeit at a low level of significance (p = 0.046).

3.5. Clock Genes Correlates: Pathway Analysis in the Cavalli Dataset

Table 4 shows the correlates for each of the clock genes and the most significant KEGG pathways associated with them. The number of correlates of clock genes (at r > 0.50) varied from zero to over 1000. The highest number of significant correlates were for the clock related gene loci CIPC (1063 genes), THRA/NR1D1 (1041 genes), FBXL21 (605 genes) and MYCBP2 (183 genes).

KEGG analysis identified the most over-represented pathway for correlates of each of these 4 genes as the ribosome pathway (CIPC, p = 1.34 x 10^-73; THRA/NR1D1, p = 1.44 x 10^-80; FBXL21, p = 2.32 x 10^-74; MYCBP2, p = 1.34 x 10^-48). The expression of these four genes was negatively correlated (r > 0.50) with numerous genes for ribosome subunits (CIPC, 57 genes; THRA/NR1D1, 59 genes; FBXL21, 44 genes, MYCBP2, 29 genes). It should be noted that the most over-represented KEGG pathway for all survival related genes (p < 0.001 by Kaplan Meier analysis) in the Cavalli dataset was the ribosome biogenesis pathway.

The KEGG pathway of phototransduction was the most significant pathway associated with correlates of RORA (with 7 contributing genes), correlates of USP2 (4 contributing genes), as well as for AANAT (4 contributing genes). For each of these 3 genes the correlates included the phototransduction genes GNAT1, GRK1, GUCY2D, RCVRN, genes whose expression was selectively elevated in the Cavalli MB subtype, Group 3 apha. Numerous phototransduction genes were over-expressed in the Group 3 alpha subtype in the Cavalli dataset; they include the rhodopsin gene (RHO) and the gene encoding phosducin (PDC). A detailed analysis of enhanced expression of genes coding for proteins of the phototransduction cascade in MB is beyond the scope of the present manuscript.

In addition for each of the three genes, RORA, USP2, and AANAT the second most significant pathway associated with their correlates was the GABAergic synapse, including the four GABAergic related genes, CACNA1F, GNB3, GNB5, and TRAK2.

3.6. Clock Gene Expression Is Related to Copy Number Gain of Chromosome 17q

A major statistical factor determining expression of clock genes was copy number gain of chromosome 17q (Table 3). This would be expected for genes located on chromosome 17q, such as the THRA/NR1D1 locus, which showed the most significant effect of 17q copy number gain (p = 1.57 x 10^-84) vs normal 17q copy number. The increase in expression was substantial (Table 3). The data however show that increased expression associated with copy number gain of chromosome 17q is also found for clock genes not on chromosome 17q, including CIPC (chr 14q), FBXL21 (chr 5q), and MYCBP2 (chr 13). Table 3 shows the statistical effect of 17q copy number gain on clock gene expression in order of significance by t-test.

PER1 is located on chromosome 17p (17p13.1), on the short arm. In isochromosome 17 a copy of chromosome 17p is lost concurrently with a gain of 17q [68]. This would account for the reduction in PER1 expression in inidividuals with copy number gain of 17 in Table 3.

4. Discussion

It has been estimated that the circadian clock regulates up to 50% of transcription in eukaryotes [69]. Evidence for a role of circadian clock genes in cancer is gradually emerging [70,71]. Here we present gene expression profiling of clock genes in the four major subgroups of medulloblastoma; furthermore the availability of the batch adjusted meta-analysis dataset of Weishaupt et al. [41] enabled a comparison of gene expression to that of NT brain tissue. We show that the prognostic significance of selected clock genes depends on the molecular subgroups of MB defined (Group3, Group 4, SHH, and WNT) and to some extent, on the 12 molecular subtypes reported by Cavalli et al. [40]. While we observed significant MB subgroup differences in expression for most the of clock genes, we focus the discussion on those that are associated with transcription of other genes, those that are most significantly associated with survival and those that have potential therapeutic relevance.

The transcription translation feedback model provides a functional paradigm that explains circadian variations in clock genes [1]. A number of the clock genes are active as transcription factors or transcription repressors. Querying Cytoscape for molecular functions of these genes we found the most significant GO (Gene ontology) molecular term associated with the genes of Table 1 was E-box binding. The reviews of Cox and Takahashi [5,19] illustrate the role of E-box binding in clock controlled genes. The E-box regulators BMAL1, CLOCK, CRY1 and PER1, as proteins (Figure 1), play a role in several cancers [72,73,74] and MB should be added to the list of those cancers. E-box binding transcription factors (EBTFs) in cancer have been reviewed by Pan, Watt and Kay [75].

Cox and Takahashi point out that that the characterization of the human clock mechanisms are not nearly complete [19]. Our analysis of clock gene expression, showing MB group specific variations in clock genes suggest that transcription translation feedback loops in clock genes are present in MB tissue and that they operate differentially in the tissues of the 4 consensus MB subtypes. We also identify the clock genes whose expression is related to copy number gain of chromosome 17q, a common genetic aberration in Groups 3 and 4 MB. Furthermore, our correlation analysis shows the major biological pathways associated with expression of clock genes The present results suggest that the MB subgroups provide novel in vitro models for further study of regulation of clock genes and proteins, and furthermore show that expression of some of the clock genes are strongly associated with patient survival in MB.

CRY1 and USP2 and Survival in Group 3 MB

The gene expression data (Figure 4A) highlight CRY1 as a defining factor in Group 3 medulloblastoma (MB), particularly in pediatric cases (Figure 4C). High CRY1 expression, noted in Group 3—the MB group with the worst survival—correlates with poor prognosis (Table 2). CRY1, a circadian rhythm regulator, is also implicated in DNA repair and tumorigenesis in cancers such as prostate cancer, making it a potential therapeutic target in Group 3 MB. High CRY1 levels likely disrupt BMAL1/CLOCK-induced transcription, further contributing to the malignancy.

CRY1 stability is regulated by ubiquitination and deubiquitination (Figure 1). USP2, a deubiquitinase, stabilizes CRY1 by protecting it from proteasomal degradation, enhancing its negative feedback role in the circadian loop. USP2 is overexpressed in Group 3 MB, particularly in young children, and is linked to survival outcomes (Table 2). Pathway analysis (Table 4) indicates that USP2 correlates significantly with phototransduction and visual perception pathways, suggesting its circadian function is dysregulated in MB.

As an oncogene, USP2 has been suggested as a therapeutic target in cancers. USP2 inhibitors have been reported [76,77] and could be tested for their effects on MB tumorigenesis and tumor cell proliferation. Additionally, the SCF-FBXL3 ubiquitin ligase complex plays a role in CRY1 degradation, with mutations in FBXL3 leading to CRY1 accumulation. Targeting the CRY1-USP2 feedback loop (Figure 1) offers a promising avenue for therapeutic intervention in Group 3 MB.

PERIOD Genes and CSNK1D in MB: Role of Phosphorylation

The PER1 and PER2 genes encode transcriptional repressors that inhibit CLOCK/BMAL (Figure 1). The USP2 protein binds and deubiquitinates PER1, regulating its nuclear localization [31,78]. Ubiquitination of period genes is mediated by the BTRC ligase [19]. CSNK1D, which encodes CK1δ, phosphorylates PER1 and PER2 proteins prior to their degradation [51]. Over-expression of CSNK1D is observed in various cancers [79], and in Group 3 and 4 MB (Figure 7), with high expression linked to poor survival (Table 2). Inhibition of CK1δ using PF-670462 has shown potential in reducing PER1, PER2, and other clock genes in rats [80] and is suggested as a cancer therapeutic [79].

The TIMELESS Connection to MB

The TIMELESS gene, encoding the TIM protein, interacts with PER2 and CRY2 [48,49]. Elevated TIMELESS expression is found in Group 3, Group 4, and WNT MB (Figure 6), and high expression correlates with poor survival (Table 2). Over-expression of TIMELESS is linked to poor prognosis in various cancers [49] and has been suggested as a therapeutic target. TIMELESS also contributes to DNA stability during replication [50], with pathway analysis (Table 4) highlighting the Fanconi anemia pathway, associated with medulloblastoma risk [81].

Clock Genes, Ribosome Connection, and Copy Number Gain of Chromosome 17q

Our analysis (Table 4) shows that the expression of four clock gene loci—THRA/NR1D1, CIPC, FBXL21, and MYCBP2—highly correlates with the ribosome biogenesis pathway. High expression of these genes inversely correlates with ribosome subunit gene transcription, with THRA/NR1D1 showing the strongest correlation (59 ribosome subunit genes). These four genes are most highly expressed in Group 4 MB (Table 1), suggesting dysregulation of clock genes affects ribosome biogenesis.

THRA/NR1D1, located on 17q, is elevated in individuals with 17q copy number gain, as shown in Table 3, along with CIPC, FBXL21, and MYCBP2. We propose the hypothesis that the isochromosome 17q aberration in MB leads to elevated clock gene expression and reduced ribosomal protein gene transcription.

NR1D1 (REV-ERBα) and CIPC influence CLOCK/BMAL1-induced transcription (Figure 1 and Figure 8), with CIPC inhibiting this pathway [64]. Jouffe et al. [82] and Pelletier et al. [83] highlighted the coordination between clock genes and ribosomal biogenesis. Our findings suggest that clock genes are linked to subgroup-specific ribosome expression in MB, supporting the view that MB cell lines could provide insights into ribosome biogenesis. Our Kaplan-Meier survival analysis indicates that the ribosome biogenesis pathway is the most significant KEGG pathway associated with all survival related genes in the Cavalli dataset.

Clock Genes, the Casein Kinase Connection, and Copy Number Gain of Chromosome 17q

CSNK1D1 is one of two clock gene loci located on chromosome 17q. Expression of both THRA/NR1D1 and CSNK1D1 are associated, at a very high statistical level, with 17q copy number gain Table 3). CSNK1D1 encodes a kinase protein that can phosphorylate a number of proteins including the PERIOD proteins and as such has been described as a major non-transcriptional regulator of circadian rhythms [84]. We suggest further study of CK1δ inhibitors (particularly those that cross the blood brain barrier) as potential therapeutic agents in individuals with copy number gain of chromosome 17q.

AANAT Another Chromosome 17q Gene Associated with Phototransduction

Since melatonin has been associated with circadian rhythms and reported to interact with BMAL1 [67], we examined the expression of AANAT, the gene encoding the rate limiting enzyme in synthesis of melatonin. AANAT is another gene located on chromosome 17q. However, no significant increase in expression was found in individuals with copy number gain of 17q (Table 3).

The expression of AANAT was elevated specifically in Group 3 (Figure 17), more specifically in Group 3 alpha. Correlation analysis showed that the biological pathways most significantly correlated with AANAT expression were the phototransduction pathway (Table 4). Our analysis relates AANAT expression to numerous phototransduction genes overexpressed in the Group 3 subgroup of the Cavalli dataset. One interpretation, that has not been disproved, is that the cell type of origin of Group 3 alpha is a photoreceptor precursor cell expressing phototransducton genes and AANAT. This developmental pattern of expression also appears to be found in the neonatal rat pineal AANAT [85].

Ubiquitin Proteasome Pathway Regulation of Clock Genes in MB

Figure 1 shows the role of ubiquitin ligases in the degradation of clock proteins, including F-box motif proteins that serve as adaptors for E3 ubiquitin ligase complexes. Several genes coding for ubiquitin ligases or adaptors are differentially expressed in MB at high significance.

FBXL15 has the lowest expression in the SHH MB subgroup (Figure 12B). FBXL15, a substrate adaptor, plays a role in resetting the circadian clock in Drosophila by degrading the TIM protein and also contributes to the degradation of SMURF1 in the SHH pathway [86,87]. FBXL21 shows one of the most significant differential expressions (Table 1). The FBXL21 protein binds to CRY proteins and opposes FBXL3's action on CRY degradation [23]. It is implicated in the central circadian clock in mammals [22], although its classification as a pseudogene in humans remains debated.

FBXW7 encodes a protein that regulates the degradation of REV-ERBα [28] and MYC [88,89]. Yeh et al. [90] discuss FBXW7 as a tumor suppressor, with its under-expression found in MB (Figure 14) for most subtypes.

MYCBP2, a gene linked to the oncogene MYC, contributes to the degradation of NMNAT2 [91] and REV-ERBalpha [66].

Summary of Pathway Analysis of Clock Gene Correlates

Table 4 summarizes the results of the KEGG pathway analysis of genes whose expression was most significantly correlated with that of clock genes. Table 4 lists the clock genes, the number of genes correlated with their expression, and the biological pathways over-represented in these lists of genes; the clock genes are listed in order of the p value for the most significant KEGG pathway. It is not a comprehensive study of all biological pathways associated with clock genes.

Table 4. Pathway analysis of clock genes.

CLOCK genes	Number of correlates r > 0.50	Most significant KEGG pathways over-represented	# of genes representing pathway	p value for pathway
THRA/NR1D1	1041	ribosome	59	1.44 x 10-80
CIPC (KIAA1737)	1063	ribosome	57	1.34 x 10-73
FBXL21	605	ribosome	44	2.32 x 10-74
MYCBP2	423	ribosome	29	1.34 x 10-48
UTS2/PER3	18	circadian rhythm	3	2.62 x 10-43
RORA	184	phototransduction	7	2.30 x 10-33
NPAS2	85	synaptic vesicle cycle	6	2.09 x 10-21
CRY2	52	GABAergic synapse	5	4.38 x 10-19
TIMELESS	31	Fanconi anemia	3	1.07 x 10-18
BTRC	34	Oocyte meiosis	4	5.61 x 10-12
USP2	161	phototransduction	4	1.74 x 10-11
AANAT	157	phototransduction	4	1.74 x 10-11
BHLHE40	148	WNT signaling	10	2.47 x 10-09
DBP	62	GABAergic synapse	4	3.98 x 10-09
FBXW7	88	morphine addiction	5	3.98 x 10-08
NFIL3	82	Hippo signaling	5	3.28 x 10-06
CSNK1D	105	prolactin signaling	4	3.45 x 10-06
CRY1	139	peroxisome	4	1.12 x 10-03
PER2	22	circadian rhythm	1	2.13 x 10-03
FBXL3	37	circadian rhythm	1	4.90 x 10-02
ARNTL (BMAL1)	8	*
ARNTL2 (BMAL2)	0	*
PER1	5	*
RORB	0	*
LINGO4/RORC	0	*
FBXL15	0	*

*no KEGG pathway detected at cutoff of r > 0..

Table 4. Pathway analysis of clock genes.

CLOCK genes	Number of correlates r > 0.50	Most significant KEGG pathways over-represented	# of genes representing pathway	p value for pathway
THRA/NR1D1	1041	ribosome	59	1.44 x 10-80
CIPC (KIAA1737)	1063	ribosome	57	1.34 x 10-73
FBXL21	605	ribosome	44	2.32 x 10-74
MYCBP2	423	ribosome	29	1.34 x 10-48
UTS2/PER3	18	circadian rhythm	3	2.62 x 10-43
RORA	184	phototransduction	7	2.30 x 10-33
NPAS2	85	synaptic vesicle cycle	6	2.09 x 10-21
CRY2	52	GABAergic synapse	5	4.38 x 10-19
TIMELESS	31	Fanconi anemia	3	1.07 x 10-18
BTRC	34	Oocyte meiosis	4	5.61 x 10-12
USP2	161	phototransduction	4	1.74 x 10-11
AANAT	157	phototransduction	4	1.74 x 10-11
BHLHE40	148	WNT signaling	10	2.47 x 10-09
DBP	62	GABAergic synapse	4	3.98 x 10-09
FBXW7	88	morphine addiction	5	3.98 x 10-08
NFIL3	82	Hippo signaling	5	3.28 x 10-06
CSNK1D	105	prolactin signaling	4	3.45 x 10-06
CRY1	139	peroxisome	4	1.12 x 10-03
PER2	22	circadian rhythm	1	2.13 x 10-03
FBXL3	37	circadian rhythm	1	4.90 x 10-02
ARNTL (BMAL1)	8	*
ARNTL2 (BMAL2)	0	*
PER1	5	*
RORB	0	*
LINGO4/RORC	0	*
FBXL15	0	*

*no KEGG pathway detected at cutoff of r > 0..

The ribosome pathway was the most significant pathway associated with clock gene correlates, particularly in relation to copy number gain of chromosome 17q. We hypothesize that one or more clock genes regulate ribosome subunit gene expression at one or more stages in the ribosome biogenesis pathway. The significant correlations of NR1D1 and CIPC expression with other genes (Table 4) suggest that REV-ERBα and CIPC may regulate transcription, including ribosome subunit genes. FBXL21 and MYCBP2, which encode ubiquitin ligase components, could regulate the degradation of clock proteins involved in feedback on the CLOCK/BMAL complex (Figure 1).

The phototransduction pathway was linked to RORA, USP2, and AANAT correlates, with over-expression seen in MB Group 3 alpha (Table 4). Hooper et al. [55] showed the Group 3 MB expression profiles were similar to that of retinal rod cell precursors. We suggest the hypothesis that RORA may stimulate transcription of phototransduction genes in the Group 3 alpha subtype.

Other significant KEGG pathways associated with clock gene correlates include the GABAergic synapse, synaptic vesicle cycle, Fanconi anemia, WNT signaling, Hippo signaling, and peroxisome pathways (Table 4), indicating potential mechanisms by which clock genes influence gene transcription in MB.

5. Conclusions

We identified clock gene expression differences in public MB datasets, with major variations between MB subgroups and subtypes. Aberrations in clock gene expression and feedback were group-specific, with major differences noted for THRA/NR1D1 and CIPC loci. Key genes encoding proteins involved in clock protein degradation, including FBXL21, FBXL3, FBXL15, MYCBP2, and the deubiquitinase USP2, also showed differential expression.

Copy number gain of chromosome 17q was associated with increased clock gene expression, both for genes on 17q and on other chromosomes, suggesting that the isochromosome 17q syndrome significantly affects clock gene transcription in MB subgroups.

The ribosome pathway was the most significant KEGG pathway associated with clock gene expression, and it also correlated with survival, likely due to 17q copy number gain.

Pathway analysis of clock gene correlates also showed associations with phototransduction, GABAergic, synaptic vesicle, WNT signaling, and Fanconi anemia pathways, with E-box binding being the most significant GO molecular term. MB should be considered a cancer in which E-box regulators play a role.

Kaplan-Meier analysis identified clock genes CRY1, USP2, CLOCK, MYCBP2, and TIMELESS as survival-related, with high expression of CRY1, USP2, MYCBP2, and TIMELESS linked to poor survival, and CLOCK to survival protection. These genes could serve as potential therapeutic targets in MB subgroups and ubtypes. Methods targeting E-box regulators, including protein degradation strategies, should be explored for those currently lacking specific inhibitors.