Association of Circadian Clock Gene Expression with Glioma Tumor Microenvironment and Patient Survival

Simple Summary Gliomas are the most common type of malignant primary brain tumors and are classified according to the cell of origin and genetic features, which can help predict the prognosis and treatment sensitivity. Improving the prognosis remains a challenge; however, chronobiology is a promising field for future works, as circadian clock genes are linked to the tumor biology and outcomes in multiple cancers, including glioma. Here, we examined the relationship of circadian clock genes, IDH mutational status, and prognosis in glioma patients by using unsupervised clustering of the expression of 13 clock genes. We further explored the expression of the clock genes across the tumor regions and cell subpopulations, highlighting the importance of the tumor microenvironment in researching circadian rhythms in cancer. Our research is important for understanding how best to target circadian rhythms to improve patient outcomes in neuro-oncology. Abstract Circadian clock genes have been linked to clinical outcomes in cancer, including gliomas. However, these studies have not accounted for established markers that predict the prognosis, including mutations in Isocitrate Dehydrogenase (IDH), which characterize the majority of lower-grade gliomas and secondary high-grade gliomas. To demonstrate the connection between circadian clock genes and glioma outcomes while accounting for the IDH mutational status, we analyzed multiple publicly available gene expression datasets. The unsupervised clustering of 13 clock gene transcriptomic signatures from The Cancer Genome Atlas showed distinct molecular subtypes representing different disease states and showed the differential prognosis of these groups by a Kaplan–Meier analysis. Further analyses of these groups showed that a low period (PER) gene expression was associated with the negative prognosis and enrichment of the immune signaling pathways. These findings prompted the exploration of the relationship between the microenvironment and clock genes in additional datasets. Circadian clock gene expression was found to be differentially expressed across the anatomical tumor location and cell type. Thus, the circadian clock expression is a potential predictive biomarker in glioma, and further mechanistic studies to elucidate the connections between the circadian clock and microenvironment are warranted.


Introduction
Gliomas are the most common malignancy of central nervous system origin, whose broad classification includes both astrocytomas (which include glioblastomas) and oligodendrogliomas [1,2]. Despite advances in clinical research, the prognosis for these diseases remains dismal, with the survival of patients with low-grade gliomas (LGG) (grades 2 and 3) averaging 5-10 years [3] and the survival of patients with high-grade gliomas (grades 4) averaging 1 to 2 years [4,5]. Multiple markers have been identified as critical for predicting the prognosis and treatment sensitivity, such as Isocitrate Dehydrogenase 1/2 (IDH) mutational status, MGMT promoter methylation status, codeletion of chromosome arms 1p and 19q, and EGFR amplification [6][7][8][9][10]. IDH mutations, in particular, confer a better prognosis across the glioma grades. For example, patients with IDH-mutated grade 4 astrocytoma have a median survival of 31 months, whereas patients with grade 4 IDH wild-type astrocytomas (glioblastomas) only 15 months [11]. These survival differences reflect fundamental differences in tumor biology that are directly related to the presence or absence of an IDH mutation. IDH is a key metabolic enzyme of the tricarboxylic acid (TCA) cycle, where it catalyzes the conversion of isocitrate into alpha-ketoglutarate (alpha-KG) [12]. The IDH mutations encountered in gliomas result in both an impairment of the ability of IDH to perform this conversion and acquisition of a novel catalytic function that enables IDH to produce 2-hydroxyglutarate (2-HG) from alpha-KG [13,14]. The accumulation of 2-HG caused by the IDH mutations alters the cellular processes responsible for gliomagenesis, including DNA methylation, gene expression levels, stem cell differentiation, cellular redox states, and DNA repair [11][12][13][14][15][16]. The dampening of the wild-type metabolic function in mutant IDH cells leads to the formation of tumors that are often initially indolent and not reliant on aerobic glycolysis (non-Warburg-like) [17]. However, many lower-grade IDHmutated astrocytomas will eventually transform into high-grade gliomas that have many of the morphologic characteristics of glioblastoma (IDH wild-type WHO Grade 4), such as vascular proliferation and necrosis, and typically exhibit an alteration of the metabolism from the neomorphic mutant IDH-driven metabolic program to the Warburg-like or aerobic glycolysis [18].
Many of the cellular processes responsible for gliomagenesis and progression are also regulated by circadian clock genes and, therefore, have daily rhythms (24-h cycles) of expression in healthy tissues [19][20][21][22]. Circadian rhythms modulate cellular functions via a transcriptional-translational feedback loop (TTFL) comprised of core clock genes: transcriptional activators CLOCK, and ARNTL and transcriptional repressors PER1/2/3 and CRY1/2 [23]. Within the central pacemaker of the brain, the suprachiasmatic nucleus (SCN), the PER gene expression is synchronized and stimulated by the onset of light. PER and CRY proteins, once transcribed, are translocated into the nucleus to interact with CLOCK and ARNTL, which then suppress the expression of PER and CRY [24]. Additional TTFLs complement the core clock components, including RORα/β/γ, REV-ERBα/β, DBP, TEF, and HLF [25]. These feedback loops impact the transcription of ARNTL and ensure an antiphase oscillation of the gene relative to PER and CRY [26]. Together, all the TTFLs generate rhythmic fluctuations in the transcriptome and modulate cellular processes via the regulation of clock-controlled genes (CCGs). Expression of these clock genes and regulation of CCGs are dependent on the tissue and cell types examined [27,28]. Across multiple cancers, disruption of the biological clock by environmental factors or through mutations in the circadian pathway can lead to an increased risk of tumorigenesis [29]. Differential expression of the core clock genes in cancer leads to alterations in the activation and/or inhibition of important oncogenic and tumor suppressive pathways [30]. Our laboratory reported a relationship between polymorphisms in ARNTL and PER2 and treatment effects in patients with brain tumors [31]. Most of these studies examined the impact of clock genes on bulk samples of tumors. However, gliomas are highly heterogeneous and display a mixture of cell subpopulations with distinct transcriptomic profiles, raising concerns that these prior results may not provide a complete picture of clock gene effects on tumor biology [32,33]. Additionally, given the infiltrative nature and prominent vascular proliferation in glioblastoma, noncancer cells in the microenvironment may further complicate the analysis. A better understanding of how circadian clock genes may impact all components of the tumors across the glioma landscape will be key to translating the findings into the clinical setting.
In recent literature, the core clock gene expression level was found to be an effective predictor of overall survival in glioma patients [30,34]. These studies, however, did not account for the IDH mutational status of tissue or cell samples when studying the effects of core clock genes in glioma. In this study, we assessed the relationship of circadian pathway genes, IDH mutational status, and prognosis in glioma patients by analyzing the molecular profiling and clinical data from The Cancer Genome Atlas (TCGA) [35] and published data from Gao and colleagues [36]. Additionally, we explored questions about the importance of the tumor microenvironment and cell types on the clock gene expression levels using publicly available RNA-seq datasets in (1) bulk samples from different anatomical locations within a tumor [37] and (2) cell types identified and profiled at the single-cell level within tumors [38,39]. Our results provide an analysis of the clock genes across the glioma landscape by accounting for the IDH mutational status while highlighting the importance of circadian rhythms in the heterogeneous microenvironments and cell types that compose glioma.

Circadian Clock Genes Define Prognostically Relevant Transcriptomic Subtypes in Glioma
To study the transcriptomic landscape of circadian clock genes in glioma and any coordinated patterns associated with prognosis, we utilized RNA-seq data from (668 grades 2-4) glioma patient samples (428 IDH-mutant, 233 IDH wild-type, and 7 of unconfirmed IDH status) from The Cancer Genome Atlas (TCGA) [40,41]. We performed an unsupervised consensus clustering of the samples based on the similarity of the transcriptomic profile of the (13) clock genes, members of the Kyoto Encyclopedia Genes and Genomes (KEGG) circadian rhythm pathway [42], which allowed us to appraise the patterns in the clock gene transcriptomic profiles and their association with known glioma prognostic markers, including age, IDH status, and sex. We determined the optimal number of clusters by considering a split of the samples into two to seven clusters and determined that four clusters (named Circadian1-Circadian4) had an optimal relative stability of the cluster assignment and composition (according to age and IDH mutational status) over multiple runs of the algorithm. This resulted in two pairs of clusters, each of which was characterized by similar age and sex distributions and IDH mutation compositions (Table 1). Circadian1 and Circadian2 were enriched for the IDH-mutated samples, whereas Circadian3 and Cir-cadian4 were enriched for the IDH wild-type samples. We further annotated the samples with their TCGA methylation classes (LGm1-LGm6) [7], which have been shown to have a strong correspondence with the primary central nervous system tumor methylation classes increasingly used alongside the histopathological diagnosis [43,44] (Table 1).
A coordinated pattern of expression can be observed within the clusters, with PER2, CRY2, NR1D1, and PER3 forming a core of positively correlated genes ( Figure S1). PER1/2/3, individually, were observed to be distinctly different across the clusters, so their expression levels were averaged for each sample to obtain a "meta" PER gene expression value reflective of the overall PER expression for that sample. Clusters Circadian1 and Circadian3 were characterized by higher levels of PER expression compared to Circadian2 and Circadian4, respectively ( Figure 1A). BHLH40/BHLH41, whose products can interact with ARNTL to suppress the transactivation of PER, had lower levels of expression in Circadian1 versus Circadian2 and Circadian3 versus Circadian4 ( Figure S2). A further analysis compared the individual clock gene expression between the TCGA clusters (comparing Circadian2 versus Circadian1 and Circadian4 versus Circadian3; Table S1).
The transcriptomic dichotomies between clusters led to hypothesizing that these subtypes represent different disease states and, thus, analyzing whether there were differences in survival between them. The Kaplan-Meier analysis revealed pairwise differences in the overall survival across the clusters, except between Circadian2 and Circadian3 (Table S2). When comparing the clusters according to their enrichment for IDH-mutated or wild-type samples, we found an association between a higher expression of PER and better overall survival ( Figure 1B,C); this association was found both for the groups enriched for IDHmutated samples (Circadian1 and Circadian2) (log-rank p < 0.0001), as well as the groups enriched for IDH wild-type samples (Circadian3 and Circadian4) (p = 2 × 10 −4 ). To further explore whether similar transcriptomic subtypes based on clock gene expression can be observed in low-grade glioma patients, we examined an additional independent dataset comprised of gene expression profiles from (195) low-grade (2) glioma samples (166 IDH-mutant, and 14 IDH wild-type) from Gao and colleagues [36]. Using unsupervised consensus clustering, we determined the optimal number of clusters to be five (Gao CC1 to CC5, Figure S3A). Cluster Gao CC1 had the highest PER expression, while Gao CC3 and Gao CC4 had the lowest levels of expression. In the absence of overall survival information, we performed an analysis of the progression-free survival on these groups. An association of the progression free survival with the overall PER expression was found in this cohort: patients in Gao CC1 had longer progression-free survival compared to those in both Gao CC3 (p = 0.009, Figure S3B) and Gao CC4 (p = 0.027, Figure S3B). , one with enriched for IDH-mutated samples and low PER expression (Circadian2, blue), one enriched for IDH wild-type samples and high PER expression (Circadian3, green), and one enriched for IDH wild-type and low PER expression (Circadian4, red). The heatmap shows the upregulation (red) and downregulation (blue) of clock genes for each cluster (Circadian1-4), as well as annotation bars for patient sex and TCGA methylation classes of the samples. Inlet graph demonstrates age distribution across the clusters. (B,C) Kaplan-Meier curves for clusters enriched for IDH-mutant (B) and IDH wild-type (C). Favorable prognosis was associated with a higher PER expression regardless of the IDH mutation composition.
To further explore whether similar transcriptomic subtypes based on clock gene expression can be observed in low-grade glioma patients, we examined an additional independent dataset comprised of gene expression profiles from (195) low-grade (2) glioma samples (166 IDH-mutant, and 14 IDH wild-type) from Gao and colleagues [36]. Using unsupervised consensus clustering, we determined the optimal number of clusters to be five (Gao CC1 to CC5, Figure S3A). Cluster Gao CC1 had the highest PER expression, while Gao CC3 and Gao CC4 had the lowest levels of expression. In the absence of overall survival information, we performed an analysis of the progression-free survival on these groups. An association of the progression free survival with the overall PER expression was found in this cohort: patients in Gao CC1 had longer progression-free survival compared to those in both Gao CC3 (p = 0.009, Figure S3B) and Gao CC4 (p = 0.027, Figure  S3B). with enriched for IDH-mutated samples and low PER expression (Circadian2, blue), one enriched for IDH wild-type samples and high PER expression (Circadian3, green), and one enriched for IDH wild-type and low PER expression (Circadian4, red). The heatmap shows the upregulation (red) and downregulation (blue) of clock genes for each cluster (Circadian1-4), as well as annotation bars for patient sex and TCGA methylation classes of the samples. Inlet graph demonstrates age distribution across the clusters. (B,C) Kaplan-Meier curves for clusters enriched for IDH-mutant (B) and IDH wild-type (C). Favorable prognosis was associated with a higher PER expression regardless of the IDH mutation composition.

Period Gene Expression Is a Biomarker of Survival in Glioma
The unsupervised clustering analysis demonstrated the existence of prognostic transcriptomic subtypes derived from all circadian clock genes and suggested that PER expression might be responsible for the survival differences. However, since the unsupervised clustering analysis resulted in groups with a mixed IDH mutation status, we next performed a survival analysis on the overall cohort stratified by the IDH status and median PER expression ( Table 1). The above median PER expression was confirmed to be beneficial for both IDH-mutated (p = 0.024, Figure 2A) and IDH wild-type gliomas (p = 0.0035, Figure 2B). Within the IDH-mutant group, we interrogated the relationship between the 1p/19q co-deletion status and the impact of the PER expression levels on survival. The Kaplan-Meier analysis of cohorts stratified by the 1p/19q co-deletion status and around the median of the PER expression levels showed that the association of PER expression with overall survival was statistically significant for tumors with intact 1p/19q (p = 0.01) (astrocytic lineage) ( Figure 2C,D) but did not achieve statistical significance in the 1p/19q co-deleted (codel) cohort (molecularly defined as oligodendroglioma). To account for age at diagnosis and sex, we performed a Cox-proportional hazards analysis that included continuous variables of PER expression, age at diagnosis, sex, and IDH mutational status. The PER expression was found to be an independent prognostic factor in glioma (HR = 0.73, p = 0.005). As anticipated from known associations, an increase in age was found to be associated with a worse prognosis (HR = 1.04, p < 0.001), as was the IDH wild-type status (HR = 5.62, p < 0.001) compared to the IDH-mutant status; sex was not found statistically significant (male versus female HR = 0.89, p = 0.473). These results indicate that the PER expression is an important predictor of survival in glioma independent of the IDH mutational status.

Low PER Expression Is Associated with Upregulation of Immune Signaling Pathways
In recent publications, low PER expression has been postulated to lead to an increased risk of tumorigenesis for multiple cancer types [45][46][47][48]. The association between PER expression and prognosis in our results led us to search for potential mechanisms responsible for worsened prognosis with dropping PER expression levels. We performed differential gene expression analysis and functional pathway enrichment analysis between TCGA clusters derived from either the unsupervised consensus clustering or stratification by IDH mutation status and PER expression level. Comparisons of clusters Cir-cadian1 versus Circadian2, and Circadian3 versus Circadian4, revealed that Circadian2 The consensus clustering analysis revealed multigene patterns characterizing some of the clusters, leading us to also assess the individual contribution of the clock genes in a multivariate Cox model of the 13 clock gene expression as continuous covariates and other key factors (age, sex, and IDH mutational status) (Table S3, full model). As before, an increase in age and IDH wild-type status were found to be associated with an increased risk of death (age HR = 1.05, p < 0.001; IDH wild-type versus mutant HR = 3.89, p < 0.001), but the sex was not statistically significant (male versus female HR = 0.83, p = 0.257). Of the circadian clock genes, PER2 (HR = 0.66, p = 0.004) and CSNK1D (HR = 2.37, p = 0.005) were associated with survival; in particular, a higher PER2 expression was associated with a decreased risk of death, whereas a higher CSNK1D expression was associated with an increased risk of death. These two genes' contributions remained significant in a model that included the patient age, sex, and IDH mutational status (Table S3, reduced model). Interestingly, the expression of CSNK1D was not found to be different between Circadian4 and Circadian3 but was decreased in Circadian2 compared to Circadian1 (Table S1).

Low PER Expression Is Associated with Upregulation of Immune Signaling Pathways
In recent publications, low PER expression has been postulated to lead to an increased risk of tumorigenesis for multiple cancer types [45][46][47][48]. The association between PER expression and prognosis in our results led us to search for potential mechanisms responsible for worsened prognosis with dropping PER expression levels. We performed differential gene expression analysis and functional pathway enrichment analysis between TCGA clusters derived from either the unsupervised consensus clustering or stratification by IDH mutation status and PER expression level. Comparisons of clusters Circadian1 versus Circadian2, and Circadian3 versus Circadian4, revealed that Circadian2 and Circadian4 (characterized by lower PER expression compared to Circadian1 and Circadian3, respectively) overexpressed immune signaling pathway genes ( Figure 3A). Examples of these pathways include granulocyte adhesion and diapedesis, dendritic cell maturation, crosstalk between dendritic cells and natural killer cells, and IL-10 signaling, among others. Functional enrichment analysis performed on the samples first stratified by the IDH mutational status then the PER expression levels similarly resulted in low PER expression being associated with an enrichment of immune signaling regardless of the IDH mutational status ( Figure 3B). Within the stratified analysis, we further assessed the functional enrichment by 1p/19q co-deletion status. We found the enrichment of immune signaling in 1p/19q intact tumors with low PER expression ( Figure 3C, right panel) but not in the 1p/19q co-deleted tumors ( Figure 3C, left panel).
The discovery that samples with low PER expression were characterized by upregulation of immune pathways led us to consider whether these differences could be an artifact of the "bulk" nature of the profiled samples, which included a varied mixture of microenvironment components and cells. To investigate this hypothesis, we downloaded computed xCell scores-estimates of the contribution of different cell types in bulk tumors derived through a transcriptomic deconvolution algorithm [49]-and compared the immune scores and stromal scores across the clusters Circadian1-Circadian4 ( Figure S4). Differences in these scores could be observed when comparing clusters with different IDH mutation status composition; however, comparisons of clusters with similar composition (Circadian1 versus Circadian2 and Circadian3 versus Circadian4) showed no differences in stromal scores, and only the pair of Circadian3 versus Circadian4 (enriched for IDH wild-type samples) was significantly different in the immune score (Table S4). Since Circadian1 and Circadian2 were not significantly different, we concluded that the cellular composition of the samples could not be fully responsible for the upregulation of immune signaling that we observed in low PER clusters.
We also investigated whether the differences in PER expression could be due to changes in DNA methylation, hypothesizing that a lower PER expression might be a result of higher levels of DNA methylation in these genes. We compared the gene level methylation of the PER genes across the identified clusters Circadian1-Circadian4. However, we did not find gene-level methylation differences that could explain the lower PER expression through the epigenetic silencing of these genes ( Figure S5).

Clock Genes Are Differentially Expressed across Anatomical Tumor Locations
The enrichment of immune signaling pathways in samples with low PER expression ( Figure 3) indicates that PER expression might affect or be affected by the immune microenvironment in glioma. The tumor microenvironment (TME) in glioma is heterogeneous [50][51][52] and characterized by a number of environmental pressures that can directly impact gene expression [53]. To assess the potential relationship between the tumor microenvironment and PER gene expression, we analyzed the clock gene expression in multiple regions within the tumor using the IvyGap dataset [37]. IvyGap contains transcriptomic profiles of samples taken from multiple distinct anatomic regions within the same patient tumor. Specifically, the dataset includes gene expression data representing 41 patients with samples taken from the following anatomic regions (with well-characterized histopathologic features) present in IDH wild-type glioblastoma (GBM): cellular tumor (CT), leading edge (LE), infiltrating tumor (IT), pseudopalisading region around necrosis (CTpan), and microvascular proliferation (CTmvp). We found that the intratumor location affects the core clock gene expression ( Figure 4A). The core clock gene expression, including the PER, were found to be higher in the CTpan and LE regions and lower in the CT and CTmvp regions ( Figure 4A). These results suggest that the tumor microenvironment impacts the clock gene expression.

Clock Gene Expression Varies across Cellular Subpopulations within Glioma Tumors
The higher immune score as defined by overexpression of immune signaling pathways in low PER TCGA clusters (Circadian2 and Circadian4) as compared to the high PER counterparts (Circadian1 and Circadian3) and the differential expression of PER based on tumor location raised the question of whether different cellular types within the glioma tumors differ in their expression of clock genes. We investigated this hypothesis using data from Darmanis and colleagues [38], who performed single-cell RNA-sequencing to assess the glioma microenvironment landscape including tumor and stromal cells from the core and periphery of (4 IDH wild-type) glioblastoma patient tumors. Supervised clustering of the profiles of neoplastic and immune cell populations revealed a pattern of upregulation of core clock genes in immune cells compared to neoplastic cells; these genes included PER1-3. Interestingly, a gradient of expression for these genes was observed within the immune cells, with cells from the periphery displaying a stronger upregulation of the core clock genes compared to the cells from the core of the tumor ( Figure 4B). This was harder to assess for neoplastic cells given the relatively small number of neoplastic cells in the tumor periphery that were profiled. Together with the IvyGap dataset, these findings suggest both the location of the cells in relation to the tumor and the classification of the cell type are important factors that should be further studied.

Clock Gene Expression Varies across Cellular Subpopulations within Glioma Tumors
The higher immune score as defined by overexpression of immune signaling p ways in low PER TCGA clusters (Circadian2 and Circadian4) as compared to the h PER counterparts (Circadian1 and Circadian3) and the differential expression of based on tumor location raised the question of whether different cellular types within glioma tumors differ in their expression of clock genes. We investigated this hypoth using data from Darmanis and colleagues [38], who performed single-cell RNA-sequ ing to assess the glioma microenvironment landscape including tumor and stromal from the core and periphery of (4 IDH wild-type) glioblastoma patient tumors. Superv clustering of the profiles of neoplastic and immune cell populations revealed a patter upregulation of core clock genes in immune cells compared to neoplastic cells; these g included PER1-3. Interestingly, a gradient of expression for these genes was obser within the immune cells, with cells from the periphery displaying a stronger upregula of the core clock genes compared to the cells from the core of the tumor ( Figure 4B). was harder to assess for neoplastic cells given the relatively small number of neopla cells in the tumor periphery that were profiled. Together with the IvyGap dataset, t findings suggest both the location of the cells in relation to the tumor and the classifica of the cell type are important factors that should be further studied. We expanded our investigation of cell-specific effects on the expression of core clock genes in tumors ( Figure S5) by taking into account the IDH status using the data published by Klemm and colleagues [39]. This published dataset represents tumor samples from 24 patients, seven of which had IDH-mutated tumors, and the remainder had IDH wild-type tumors. We examined the impact of PER2 expression in these cells, as well as ARNTL, as it is often antiphasic to the period expression in vitro ( Figure 5A). We assessed the PER2 expression levels in microglia (MG), myeloid-derived macrophages (MDM), and CD45-negative cells (CD45n) between the IDH-mutant and IDH wild-type tumors ( Figure  5). The PER2 expression is lower across all cell types in the IDH wild-type tumors; however, the difference is only statistically significant in the MG (p = 0.03).
patients, seven of which had IDH-mutated tumors, and the remainder had IDH wild-type tumors. We examined the impact of PER2 expression in these cells, as well as ARNTL, as it is often antiphasic to the period expression in vitro ( Figure 5A). We assessed the PER2 expression levels in microglia (MG), myeloid-derived macrophages (MDM), and CD45negative cells (CD45n) between the IDH-mutant and IDH wild-type tumors ( Figure 5). The PER2 expression is lower across all cell types in the IDH wild-type tumors; however, the difference is only statistically significant in the MG (p = 0.03).

Discussion
The discovery of targetable universal mechanisms that affect cancer growth remains an elusive challenge but would provide a promising vehicle for identifying treatments of

Discussion
The discovery of targetable universal mechanisms that affect cancer growth remains an elusive challenge but would provide a promising vehicle for identifying treatments of aggressive cancers, such as glioma. Alterations in circadian clock genes have been directly linked to clinical outcomes such as treatment toxicities [31] and overall survival [30,34] in glioma patients. Even though IDH mutational status has distinct biologic and prognostic importance in malignant gliomas, to our knowledge, no prior clock gene/glioma studies have taken into consideration the IDH mutational status. Analyzing molecular profiles from patient tumors available in several public datasets, we discovered distinct transcriptional subtypes based on 13 circadian clock genes and that these subtypes could not be explained by the IDH mutation status alone but may be reflective of different disease states. This hypothesis was supported by the finding that the prognosis of the patients whose tumors displayed different clock subtypes were also found to be different. Low expression of the period genes (PER) was associated with poor prognosis in both IDH-mutant and wild-type groups of astrocytic lineage, although the association did not rise to statistical significance for the oligodendroglial (1p/19q co-deleted) lineage. Importantly, absolute PER expression levels were not found to be directly related to prognosis: for example, in the TCGA clustering analysis, patients from cluster Circadian3 enriched for IDH wild-type samples had a higher expression of PER but worse prognosis than patients from cluster Circadian1 enriched for IDH-mutated gliomas ( Figure S2). While this result is expected due to the known advantage of IDH mutations on survival and the lower age of the patients harboring IDH-mutant tumors, it also highlights the need to assess PER expression relative to the IDH mutant class of the tumors.
Our functional pathway analysis of cohorts with lower PER expression and worse relative survival showed the upregulation of immune signaling pathways and suggested that immunological mechanisms may be responsible for these differences in survival. In this context, circadian rhythms directly modulate the immune system in the periphery, regulating a myriad of immune functions, such as cytokine and chemokine secretion [54,55]. Disruption of these rhythms caused by chronic jetlag [56,57] or the genetic deletion of core clock genes [58,59] has a negative impact on the functionality of immune cells and can lead to premature death. Age is another factor that can disrupt clock gene expression and circadian rhythms at the organism level leading to premature death [60]. Although our analysis indicates that the PER expression level is an additional independent prognostic factor, alongside age and IDH mutation status, the untangling of any interaction between these factors deserves further examination through time series studies. An intriguing connection between age and period length at the cellular level has been recently proposed by Li et al., who examined the circadian period length in single cells and clonal cell lines and found that a longer period is associated with the increased intercellular heterogeneity and transcriptomic noise often observed in aging cells and cancer [61]. Within solid cancers [62][63][64], including glioma [30], circadian misalignment is associated with increased tumor growth rate and decreased survival. These findings have been hypothesized to be caused by suppression of the circadian-related immune function. A recent publication tested the hypothesis in a melanoma mouse model [65], demonstrating a loss of daily rhythms in macrophage infiltration within the tumor microenvironment and spleen under chronic jetlag conditions. By analyzing the IvyGap dataset, we found that there is higher expression of core clock genes, including PER, in specific anatomical tumor microenvironments: the leading edge and the pseudopalisading region around the necrosis of the tumor. These subregions are associated with different immune signatures and oxygenation levels [66]. Greater interrogation of the specific immune cells within these regions would provide insight into the negative circadian effects caused by the hypoxia conditions [67] or tumor cell proximity.
Most cells within the body express circadian rhythms in the transcriptome; however, the genes that are regulated and their pattern of expression are dependent on the tissue type [28]. Within the brain, the regions can express different and sometimes antiphasic patterns of clock and clock-controlled genes compared to the master biological clock found in the suprachiasmatic nucleus [28,68,69]. In vitro, the different cell types that compose the brain, including neurons [70], glial [71], epithelial [72], and immune cells [73], have been shown to express rhythms of clock genes. When examining circadian rhythms in the glioma, tumor cells have been the primary and only focus of cell culture works. These studies have demonstrated the importance of clock genes on tumor growth [30] and the optimal timing of treatment, including temozolomide [74] and radiation [75,76]. Additionally, our findings here demonstrate that clock genes of immune cells found in the tumor microenvironment are different based on IDH mutational status, specifically microglia. Our data do not include a time series, so alterations in circadian rhythms cannot be deduced. The bulk tissue isolated from patients often does not have time stamp information on the collection times, and this is a major limitation of genomic studies using publicly available datasets. However, the rewiring of circadian rhythms by other solid tumors has been shown within support cells found in the local microenvironment [77] and systemically, impacting circadian metabolic function in the liver [78]. Glioma cells have been shown to directly suppress the immune function in the periphery, which demonstrates the capability of tumors to impact distant foci and suggests a mode by which circadian rhythms may be altered within immune cells [79]. Recently, the relationship between circadian disruption in tumors and microglia invasion into the glioma was explored in mice and cell cultures by Chen et al. [80]. The authors demonstrated that the suppression of CLOCK and BMAL in tumor cells improved the survival and reduced migration of the microglia. These findings, however, only focused on the tumor rhythms and did not examine the impact of clock genes in the microglia. Understanding the circadian rhythms within and beyond the neoplastic cells will be the key to developing clock gene-based strategies to reduce the treatment side effects and improve survivorship. Future works should broaden the focus of clock genes in glioma research to include the microenvironment and its relationship to tumor cells.

Data Sources
RNA-seq expression data from TCGA (668 grades 2-4) primary glioma samples were downloaded using the R TCGAbiolinks library [81,82]. The data were preprocessed, normalized, and filtered as recommended in R's TCGAWorkflow library [83].
Array hybridization mRNA expression data and progression-free survival (PFS) patient information from 195 patients diagnosed with low-grade (2) glioma enrolled in the EORTC22033-26033 clinical trial were downloaded from Gene Expression Omnibus (GEO) (accession: GSE107850 [36]).
scRNA-seq data of 4028 single cells from 4 primary GBM tumors were downloaded from GEO (accession: GSE84465; [38]). Samples from each patient were taken from the tumor core and peritumoral tissue (cortex). The data were normalized using quantile normalization. Individual gene queries were performed using the online data portal (http://gbmseq.org/, accessed on 16 January 2021).
RNA-seq expression data of immune cell populations within the microenvironment of primary glioma samples were downloaded from Klemm et al. [39]. Cell sorting of clinical samples was performed prior to sequencing for identification of immune cell populations. The dataset represented 48 patients diagnosed with brain tumors, 24 of which were diagnosed with glioma grades 2-4. The glioma patient cohort included 7 patients with IDH mutated tumors and 17 with IDH wild-type tumors. (https://joycelab.shinyapps. io/braintime/, accessed on 30 April 2020).
xCell bulk tumor deconvolution scores for the TCGA samples were downloaded from the xCell portal (https://xcell.ucsf.edu/, accessed on 5 August 2020) [49]. The predicted, sample-level Immune Score and Stroma Score were used to assess the overall differences in the contribution of immune and stroma elements to the bulk samples from the discovered cohorts of interest.

Hierarchical Clustering
Data were analyzed using R's TCGAbiolinks and ConsensusClusterPlus libraries. Unsupervised consensus clustering of the samples was performed based on the transcriptomic profiles of 13 clock genes (BHLHE41, BHLHE40, CLOCK, NPAS2, ARNTL, NR1D1, CRY2, PER3, PER2, PER1, CSNK1E, CRY1, and CSNK1D); the genes used are members of the KEGG_CIRCADIAN_RHYTHM_MAMMAL pathway of the Kyoto Encyclopedia of Genes and Genomes (KEGG) knowledgebase and were obtained from the Molecular Signatures database (MSigDB) (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp, accessed on 25 January 2018) [42,85,86]. Clustering into 2-7 groups was considered. Clustering into 4 groups was found to be in the optimal range as assessed using the ConsensusCluster-Plus [87] analysis of the stability of cluster membership over multiple runs of the algorithm, as well as the drop in improvement of the Cumulative Distribution Function with each addition to the number of clusters. A similar analysis was performed on the LGG dataset from Gao et al. [36].
Unless otherwise specified, R's ComplexHeatmap [88] library was used for visualization of all hierarchical clusters, including all heatmaps depicting transcriptomic differences across cohorts, tumor locations, cell types, and xCell scores.

Survival Analysis
Kaplan-Meier analysis was performed to assess any differences in overall survival (OS) between sample clusters discovered through unsupervised clustering of the TCGA samples. For the LGG cohorts, in the absence of OS information, PFS differences were instead assessed. The significance of the differences was determined using log-rank tests; p < 0.05 were considered significant in pairwise comparisons. Analogous analyses were performed for the cohorts stratified by IDH mutation status, 1p19q codeletion status, and PER1/2/3 expression levels dichotomized around the median expression.
Cox-proportional hazards analysis was performed to assess the predictive power of clock gene expression levels on OS of TCGA patients and PFS for the LGG cohort. Models were constructed with covariates age at diagnosis, IDH mutational status (and for TCGA 1p/19q codeletion status), and combined expression of PER(1/2/3) (averaged per each sample), as well as with covariates age at diagnosis, IDH mutational status (and for TCGA 1p/19q codeletion status), and all 13 clock genes as separate covariates. Likelihood ratio test of the overall model was assessed, as well as the Hazard Ratio (HR) of individual covariates. HR with a 95% Confidence Interval not crossing 1 were considered significant (or interchangeably, the individual covariates' contributions to the model were considered significant when associated p < 0.05).

Differential Gene Expression and Functional Pathway Enrichment Analysis
Differential gene expression analysis between TCGA clusters was performed using the TCGAbiolinks library, whose functionality is tailored to the differential gene expression analysis of samples with experimental and batching characteristics of TCGA. Two sets of differential expression analyses were performed: between clusters discovered through consensus clustering (here, the comparisons were done for clusters with similar molecular backgrounds, comparing the clusters enriched for IDH-mutant samples to each other, and comparing the clusters enriched for IDH wild-type samples to each other), as well as between clusters that have the same IDH mutational status but are split according to the median expression of PER. Fold change (FC) > 1.5 and adjusted p < 0.05 were used as the criteria for a significant differential expression of the genes.
Gene lists obtained from the differential gene expression analysis were submitted for functional pathway enrichment analysis using TCGAbiolinks and annotations from the Database for Annotation, Visualization and Integrated Discovery (DAVID) knowledgebase [89,90]. Pathways with adjusted p < 0.05 of enrichment were considered significantly enriched.

Other Comparisons
Comparisons of the Immune and Stroma Scores between the TCGA clusters were assessed using R's limma package [91]: fitting a linear model, computing moderated tstatistics, moderated F-statistic, and log-odds of differential expression. The results were adjusted for multiple comparisons, and an adjusted p < 0.05 was used as the criterion of significance.
Individual gene expression differences were depicted using R's ggplot2 boxplot functionality and assessed using t-tests through the ggpubr library. Barplots of gene expression were generated by running queries at http://gbmseq.org/ (accessed on 16 January 2021) for the Darmanis et al. dataset [38] and at https://joycelab.shinyapps.io/braintime/ (accessed on 30 April 2020) for the Klemm et al. dataset [39]. Where t-tests were used, p < 0.05 was considered to indicate statically significant differences.

Conclusions
Our results highlighted the association of the circadian clock gene expression and glioma patient outcomes across the glioma landscape and independent of the IDH mutational status of glial tumors. Importantly, we showed a complex interplay of the tumor microenvironment and the circadian clock, demonstrating that expression of the clock genes is highly plastic and varied across tumor regions and cell types. Taken together, our findings support that circadian clock transcriptomics are directly reflective of clinically relevant disease states and therefore represent important potential biomarkers and targets for therapeutic interventions.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/cancers13112756/s1: Figure S1: Graphs demonstrating the coordinated patterns of expression of clock genes in the TCGA data and unsupervised consensus clusters. Figure S2: Clock gene expression in the unsupervised clustering analysis. Figure S3: Prognostic relevance of PER in unsupervised clusters for LGG dataset. Figure S4: Contribution of different cell types in the bulk tumor derived through a transcriptomic deconvolution algorithm. Figure S5: Differences in PER2 and ARNTL expression patterns between different cells of the tumor microenvironment. Figure S6: Clock gene expression in the subgroups separated by methylation status. Table S1. Differences in clock gene expression between high and low PER-expressing circadian clusters with predominantly IDH mutant or wild type. Table S2: Kaplan-Meier analysis of the overall survival between unsupervised consensus clusters of TCGA data. Table S3: Summary of the full Cox model of the 13 core clock gene expression and reduced model of significant clock genes for the TCGA dataset. Table S4