Circadian Clock Genes Act as Diagnostic and Prognostic Biomarkers of Glioma: Clinic Implications for Chronotherapy

Gliomas are the most common primary intracranial tumors and closely related to circadian clock. Due to the high mortality and morbidity of gliomas, exploring novel diagnostic and early prognostic markers is necessary. Circadian clock genes (CCGs) play important roles in regulating the daily oscillation of biological processes and the development of tumor. Therefore, we explored the influences that the oscillations of circadian clock genes (CCGs) on diagnosis and prognosis of gliomas using bioinformatics. In this work, we systematically analyzed the rhythmic expression of CCGs in brain and found that some CCGs had strong rhythmic expression; the expression levels were significantly different between day and night. Four CCGs (ARNTL, NPAS2, CRY2, and DBP) with rhythmic expression were not only identified as differentially expressed genes but also had significant independent prognostic ability in the overall survival of glioma patients and were highly correlated with glioma prognosis in COX analysis. Besides, we found that CCG-based predictive model demonstrated higher predictive accuracy than that of the traditional grade-based model; this new prediction model can greatly improve the accuracy of glioma prognosis. Importantly, based on the four CCGs' circadian oscillations, we revealed that patients sampled at night had higher predictive ability. This may help detect glioma as early as possible, leading to early cancer intervention. In addition, we explored the mechanism of CCGs affecting the prognosis of glioma. CCGs regulated the cell cycle, DNA damage, Wnt, mTOR, and MAPK signaling pathways. In addition, it also affects prognosis through gene coexpression and immune infiltration. Importantly, ARNTL can rhythmically modulated the cellular sensitivity to clinic drugs, temozolomide. The optimal point of temozolomide administration should be when ARNTL expression is highest, that is, the effect is better at night. In summary, our study provided a basis for optimizing clinical dosing regimens and chronotherapy for glioma. The four key CCGs can serve as potential diagnostic and prognostic biomarkers for glioma patients, and ARNTL also has obvious advantages in the direction of glioma chronotherapy.


Introduction
Glioma is the most common primary intracranial tumor, accounting for more than 70% of malignant brain tumors, and is associated with high mortality and morbidity [1]. Gliomas commonly occur in the cerebellum, brainstem, and cortex [2,3] and are currently classified into low-grade glioma (LGG, grade II and III) and glioblastoma (GBM, grade IV) [4]. Notably, a considerable proportion of LGG develops rapidly in a short time and deteriorates into GBM. Analysis of the isocitrate dehydrogenase status and 1p/19q codeletion can contribute to personalized diagnosis and prognosis [5].
Despite these findings, the poor prognosis of glioma still urgently requires the discovery of novel biomarkers for the early detection of glioma and to improve patient survival.
The circadian clock is a molecular timekeeping mechanism that exists in all creatures and regulates the daily oscillations of biological processes and behavior [6]. Circadian clock genes (CCGs) form feedback loops to maintain normal function of the oscillating system and to keep it in sync with the environmental cycle. To date, 13 CCGs have been well identified, including ARNTL/BMAL1 [7], CLOCK [8], CRY1, CRY2 [9], DBP, NPAS2 [10], NR1D1, NR1D2, PER1, PER2, PER3 [11], RORA [12], and TIMELESS [13]. Increasing evidence has revealed that CCGs play critical roles in DNA damage and repair, cell proliferation and metastasis, immunity, and tumorigenesis [14,15]. Other studies have shown that intervening with the circadian clock could effectively influence stem cell growth, tumor aggressiveness, and drug delivery to improve therapeutic outcomes [16][17][18][19]. For instance, activation of PER1 effectively inhibited the progression of pancreatic cancer [20]. NR1D2 promoted the proliferation, migration, and invasion of GBM cells through specific targets [21]. Although a recent study showed that CCGs have significant prognosis value in glioma [22], the influence of CCGs rhythmic fluctuations on prognosis has not been conclusively determined. Therefore, our study focused on the effect of genes rhythm on glioma prognosis.
In this study, we systematically analyzed the circadian expression patterns of CCGs in brain tissues. The differential expression analysis was used to identify possible diagnostic markers. Kaplan-Meier survival curves, Cox proportional hazards model, and nomogram were used to assess the prognostic value of CCGs. Besides, CCG-involved cancer-related pathways and coexpression network analysis were also launched. Furthermore, CCG-related cellular sensitivity to temozolomide was demonstrated, which may facilitate glioma chronotherapy.

Datasets and Data
Availability. Profiles of patients with glioma were downloaded from The Cancer Genome Atlas (TCGA) and the Chinese Glioma Genome Atlas (CGGA). TCGA is a project launched by US in 2006 that contains data on various human cancers, while CGGA is a biobank of glioma samples for the Chinese population. A total of 698 cases from TCGA and 1,018 cases from the CGGA were collected. Glioma patients with missing OS values or OS < 30 days were excluded to reduce statistical bias in our analysis. The circadian rhythm data of CCGs in brain tissues of mice were collected from GSE54652 of the GEO database.

Pan-Cancer Analysis and Gene Correlation Analysis.
GSCA (genomic cancer analysis platform) was used for pan-cancer analysis. The Sangerbox was used to show coexpression patterns among CCGs in brain tissues. The clock correlation distance (CCD), an algorithm that have been constructed in the literature [23], was used here to infer circadian clock progression in a group of samples based on the coexpression of 13 clock genes. Gene expression data were obtained from GSE54652 and TCGA. Red notes indicate that the effect of high expression on survival risk is high and vice versa is indicated by the green notes. Pearson's correlation was used to assess the internal relationship of CCGs.

CCG Rhythmic Expression in Mouse Brain
Tissues. CCG expression data in the brain tissue downloaded from GEO's GES54652. This dataset was collected from different organs of mice at two hourly intervals over a 48-hour period. Microarray data were used to show daily changes in gene expression. To determine if the data is rhythmic, R language package called "JTK_CYCLE" was used to detect the rhythm of each gene. The JTK_CYCLE package, with parameters set to fit time-series data to exactly 24 h periodic waveforms and statistical analysis, was used. When P < 0:05, the gene was considered to be rhythmically expressed [24]. GraphPad Prism (v7.0) was used to draw graphs representing cyclical fluctuations. The rhythm graphs of the CCGs were also constructed using GraphPad Prism.

Differential Expression Analysis.
Gene expression profiles of brain tissues in glioma patients with clinical grade were extracted from two datasets (698 tumor samples from TCGA and 1,018 tumor samples from the CGGA). Gene expression levels were normalized using the edgeR and limma packages for TCGA and CGGA data [25,26], respectively, in R software. P value < 0.05 was considered differentially expressed genes (DEGs).

Analysis of Methylation, Copy Number Variation (CNV)
, and Single-Nucleotide Variation (SNV). GSCA was used to reveal the role of CNV, methylation, and SNV in the regulation of CCG expression. Significance was set at FDR ≤ 0:05. The Pearson correlation coefficient was used to analyze the correlation between CNV, methylation, and gene expression.
2.6. Survival Analysis. The CCG expression in patient tumor samples from the two datasets with survival data (TCGA and CGGA) was included in the survival analysis. Survival curves were plotted by the Kaplan-Meier method using R software, and the median value of gene expression was set as the group cut off separating the high-expression and low-expression groups. Univariate and multivariate Cox regression analyses were calculated by "survival" package. The hazard ratio (HR) and P value of each DEG of CCGs were determined based on gene expression and overall survival of patients using the univariate Cox regression model with the survival package in R software. HR > 1 and HR < 1 indicate that the higher expression of the gene is associated with worse or better overall survival, respectively. 2.7. Prediction Models. R software was used to create Cox proportional hazards prediction models based on CCG expression levels and overall survival. Accordingly, the following formula determines the risk score of each patient: where n is the number of involved genes, Coef is the coefficient of every gene, and Exp is the expression level (log2    [27] was used to assess the correlation coefficients between the CCGs and other hub genes. The coexpressed genes with correlation > 0:5 were extracted. Cytoscape (v 3.8.0) visualized the correlation between key CGGs and coexpressed genes in the network.
2.10. Immune Infiltrate Analysis. Six immune cell types, B cells, CD4+ T cells, CD8+ T cells [28], neutrophils, dendritic cells (DCs) [29], and macrophages [30], in the tumor micro-environment have been previously reported. We used the Tumor Immune Estimation Resource (TIMER), an established algorithm, to estimate the correlation between immune cells and gene expression. Spearman's correlation and estimated statistical significance between CCG expression and immune cell signature infiltration are shown in scatter plots.  18     The HRs and 95% confidence intervals (CIs) were calculated to identify genes associated with overall survival. Pearson's correlation analysis was used to calculate the correlation coefficient. P value < 0.05 or FDR < 0:05 were considered statistically significant.

Pan-Cancer Analysis and Coexpression Patterns of CCGs.
The prognostic influence of CCGs on the survival risk of several cancers was performed (Figure 1(a)).
LGG exhibited the highest P value in pan-cancer analysis, with 10 significant prognostic CCGs in total. The results indicated that among all cancers, glioma was affected mostly by CCGs.   Next, we revealed the coexpression of CCGs in brain tissues (Figures 1(b)-1(d)), according to the coexpression correlations in brain, CCGs were apparently divided into two clusters: cluster I includes ARNTL, NPAS2, CLOCK, and TIMELESS, and another includes CRY1, CRY2, DBP, NR1D1, NR1D2, PER1, PER2, PER3, and RORA. Similar coexpression relationship were found by using the CCD algorithm [23], indicating that the rhythmic oscillation patterns might be also alike in normal brain tissue and tumor samples. To further confirm the coexpression relationship in brain tissues, gene correlation maps were performed, which indicates that there is a high positive correlation between CCGs of the same cluster (Figure 1(e)).

Rhythmic Expression of CCGs in Brain
Tissues. Based on the day-night expression data of mouse genes in the GSE54652 dataset, we performed JTK analysis using bioinformatics methods in R software, in which P < 0:05 was considered a gene with rhythmic expression. The results showed the rhythmic expression of CCGs in consecutive 48 hours. Apparently, most of cluster CCGs exhibited a significant rhythmic expression pattern in the cerebellum (Figures 2(a) and 2(c)), and in brainstem (Figures 2(b) and 2(d)), respectively. The results covered two circadian cycles. There was little difference in CCG rhythmic expression between the cerebellum and brainstem, suggesting that CCGs rhythmic patterns in brain tissue are stable [31,32]. The fluctuation patterns of the two clusters were opposite, with cluster I CCGs reaching the peak in the morning and cluster II CCGs reaching the peak in the evening.

CCGs Were Differentially Expressed in the Grade of Glioma.
To investigate the role of CCGs in all grades of glioma, including grade 2-4, we firstly analyzed the differential expression of CCGs in glioma grade and found that 11/13 CCGs, except for CLOCK and PER1, were recognized as DEGs in TCGA (Figures 3(a) and 3(b)) and CGGA (Figures 3(c) and 3(d)), respectively. Interestingly, two clusters' genes had opposite expression trends, which expression levels of cluster I and CRY1 were increased with the grades, whereas cluster II and CLOCK decreased. We then examined the possible regulatory mechanism of CCG expression by analyzing methylation and copy number variation (CNV). All DEGs were negatively regulated by methylation (FDR < 0:05, Figure 4(a)), except for PER3. Besides, except for PER1, NPAS2, CLOCK, and ARNTL, the other 9 DEGs were positively regulated by copy number variation (FDR < 0:05, Figure 4(b)). The single-nucleotide variation of CCG was not observed in glioma (Figure 4(c)). The Venn diagram clearly demonstrated that 8 CCGs, including most cluster II CCGs, were regulated by both methylation and CNV (Figure 4(d)). The combination of methylation and CNV mediated most of the downregulation of cluster II CCGs. Collectively, these results implied that methylation and CNV mainly contributed to the differential expression of CCGs in gliomas.

CCGs as Prognostic Biomarkers for Glioma Patients.
The prognostic abilities of a single CCG were analyzed and survival curves showed that 10 differential expressed CCGs were significantly associated with the survival rates of glioma patients (P < 0:05). Importantly, high expression of cluster I and CRY1 was associated with poor survival of glioma patients in both the TGCA and CGGA databases (Figures 5(a) and 6(a)). On the contrary, lower expression of cluster II and CLOCK was related to low survival time (Figures 5(b) and 6(b)). ROC curves were plotted to assess CCG predictive accuracy (AUC value ≥ 0:6 as predictive) [33], which showed that most of CCGs had accurate predictive ability in TCGA ( Figure 5(c)) and CGGA (Figure 6(c)). Therefore, these CCGs were potential prognostic biomarkers for patients with glioma.

Prediction Model of CCGs by Cox Regression Analysis.
Univariate and multivariate Cox regression analyses were performed to investigate the prognostic role of CCG combinations. As shown in Tables 1 and 2, univariate Cox regression analysis demonstrated that the poor survival of patients was related to grade, cluster I, and CRY1, with HR > 1 (HR > 1 represents a high-risk factor) [34], whereas 6 CCGs in cluster II, except for PER1, were considered low-risk factors. In addition, multivariate Cox regression analysis indicated that the combination of 6 CCGs, including ARNTL, NPAS2, CRY2, DBP, RORA, and TIMELESS, was significant for patient survival in both the TCGA and CGGA databases (Tables 1 and 2). Accordingly, a 6-CCG-based predictive model was constructed, and the survival curve and ROC curve clearly showed that the 6-CCG-based prediction model had higher accuracy than that of the prediction model based on grade in TCGA (Figure 7(a)) and CGGA (Figure 7(b)). The accuracy of ROC analysis was the similar for all grades of glioma ( Fig S1). Furthermore, nomogram analysis was performed, which revealed that 4 of the 6 CCGs in    d)). Based on the finding that the expression peak phase of CCGs shifted by~12 hours between the mouse and baboon [35], CCG expression fluctuations throughout the day in mice and humans, and appropriate sampling time were redisplayed (Figure 7(e)). For instance, ARNTL and NPAS2 showed peak of expression in mouse in the morning, whereas in human, the peak phase occurred in beginning of the evening (about 8:00 pm). On the contrary, CRY2 and DBP exhibited opposite expression patterns.

Pathway Enrichment Analysis of the Four Key CCGs.
GSCA and GSEA were further used to analyze the possible CCG-involved molecular mechanisms that the pie chart of GSCA demonstrated several pan-cancer pathways (Figure 8(a)), while GSEA showed key signaling pathways (Figure 8(b)). In most cases, cluster I genes (ARNTL, NPAS2) had opposite effects against cluster II genes (CRY2, DBP) in the same pathway. For instance, ARNTL activated Wnt signaling pathway which was inhibited by CRY2 and DBP. Besides, cluster II mainly activated Wnt, mTOR pathway and inhibited cell cycle, DNA damage pathway to promote proliferation and tumor development. Moreover, key genes in these classic cancer-related pathways were found and their expression levels were analyzed, which also exhibited highly similar rhythmic fluctuations to CCGs ( Figure 8C). Together, these results implied that CCGs regulated these cancerrelated pathways in progression of glioma.

Coexpression Network Analysis of the Four Key CCGs.
WGCNA was performed to identify the genes that were coexpressed with the four key CCGs (ARNTL, NPAS2, CRY2, and DBP) in glioma (Figure 9(a)). CCGs were displayed as red nodes and the coexpressed genes were marked as blue nodes. Besides, several genes having rhythmic expression were circled in red. It was found that cluster I and cluster II genes were tightly linked to each other, respectively. The correlations between the red circle genes (representative rhythmic genes) and CCGs were shown by scatter plots (Figure 9(b)). Correspondingly, these highly positively correlated genes also demonstrated rhythmic fluctuations in the cerebellum and brainstem (Figures 9(c) and 9(d)). Among them, TEF was tightly associated with CRY2, DBP, and NR1D2 and thus had a similar rhythmic expression as cluster II genes. Other arrhythmic genes, such as PDK2, have prognostic abilities in lung adenocarcinoma [36]. Taken together, these findings suggested that CCGs can regulate the rhythm of other genes to interfere with the progression of glioma.

Immune Infiltrate Analysis of the Four Key CCGs.
Immune infiltration is a pivotal biological characteristic of tumors. To investigate whether the prognosis of CCGs is related to immune infiltration, the correlations between these CCGs and immune cells in glioma were explored using the online tool TIMER. Four immune cell types with the close correlation to key CCGs (expression level adjusted by log2 TPM) were shown in Figure 10. The result showed that ARNTL and NPAS2 was positively correlated with CD8+ T cells, whereas CRY2 and DBP were negatively correlated with neutrophils and macrophages. This revealed that same cluster CCGs affect patient prognosis in a similar way in immune infiltration. CRY2 is highly correlated with immune cells (jcorrelationj > 0:4), suggesting that CRY2 may affect prognosis through immune infiltration more than other CCGs. Table 3 showed the important internal biomarkers of immune cells, in which CCGs affect the state and function of immune cells through activation or inhibition of these biomarkers, thus affecting patient survival. Spearman's correlation in the table indicated that CRY2 had prominent effect on these biomarkers.

Drug Sensitivity Analysis of the Four Key CCGs.
Drug sensitivity analysis provided a bubble plot of the relationship between the four key CCGs and various cancer drugs. (Figure 11(a)). We used bioinformatics to find the drug sensitivities most associated with the four CCGs in glioma based on GDSC, which is a database containing data on drug sensitivity and gene expression profiles of hundreds of tumor cell lines. The results showed that four CCGs were associated with the sensitivities of several drugs, among which temozolomide was found to be an existing clinical drug for glioma [37]. The blue circle showed negative correlation meant that the higher the expression of ARNTL, the lower IC 50 values (represented higher sensitivity) of temozolomide in cancer cells. Although other CCGs presented similar rhythms, they were not significantly associated with drug sensitivity in the bubble plot. Previous Notes: characteristics with P < 0:05 in the univariate analysis were further screened in the multivariate analysis. Exp: expression; HR: hazard ratio; CI: confidence interval; * P < 0:05, * * P < 0:01, and * * * P < 0:001.   14 BioMed Research International literature suggesting the therapeutic effect of temozolomide on glioma varies with ARNTL rhythmic expression, the effect was the best in peak and the worst in the trough [38]. Accordingly, the time-of-temozolomide sensitivity curve was plotted (Figure 11(b)), which clearly demonstrated the different therapeutic potentials of temozolomide at different time points of administration.

Discussion
In this study, we systematically analyzed and explored the role of CCGs in the diagnosis, prognosis, and chronotherapy of glioma. Key CCGs, including ARNTL, NPAS2, CRY2, and DBP, acted as potent diagnosis and prognosis biomarkers of glioma patients and rhythmically regulated cancer-related  signaling pathways and coexpressed genes to affect drug sensitivity and possible clinic outcomes. The CCG-based predictive model demonstrated higher predictive accuracy than that of the traditional grade-based model; this new prediction model can greatly improve the accuracy of glioma prognosis. Another study also confirmed the ability of CCGs as prognostic markers for gliomas, suggesting that CCGs' expression affect immunity and cell cycle [22]. However, the unique rhythm expression of CCGs has not been analyzed in the paper. Our study found that the oscillating expression of CCGs has a great impact on diagnosis and prognosis. Rhythmic expression was firstly introduced into the CCG-based predictive model. Appropriate sampling time could greatly improve the ability of early diagnosis and obtain a better prognosis, while inappropriate sampling time may lead to misdiagnosis and delay treatment. In addition, we also

16
BioMed Research International predicted the optimal administration time of temozolomide, which provides an effective reference for the chronotherapy of glioma. The differential expression of CCGs in glioma grading confirmed that cluster I genes (ARNRL, NPAS2, and TIME-LESS) and CRY1, with cluster II genes (CRY2, DBP, NR1D1, NR1D2, PER2, PER3, and RORA) showed significantly opposite expression trends in brain tissues ( Figure 3). Besides, the expression of eight CCGs (CRY2, DBP, NR1D1, NR1D2, PER2, PER3, RORA, and TIMELESS) was negatively regulated by methylation and positively regulated by CNV in glioma, which is consisting with previous studies that DNA methylation and CNV can promote abnormal gene expression [39][40][41]. Furthermore, we found that single CCG has high independent predictive accuracy and may be a potential prognostic biomarker for glioma ( Figures 5 and 6). Our CCG-based prediction model had higher accuracy than the traditional grade-based model (Figures 7(a) and 7(b)) [42].
Four key CCGs, including ARNTL, NPAS2, CRY2, and DBP, showed opposite behavior in differential expression and prognosis (Figures 7(c) and 7(d)). Previous studies have also confirmed the opposite effects of NPAS2 and CRY2 on cancer prognosis [43]. Moreover, the four key CCGs-based risk scores varied along the expression fluctuations. Patients sampled at different times had great impacts on the total predictive points and survival rates (Figure 7(e)). A recent study also indicated that ARNTL and PER1 showed peaks of expression in mice in the morning and evening, respectively, whereas in baboon, the peak phases occurred in opposite time point [35]. Therefore, it was reductive that sampling glioma patients at beginning in the evening (about 8:00 pm) could obtain higher predictive sensibility whereas sampling in the morning might result in a misleading diagnosis.
In addition, the four key CCGs were related to several pathways, including cell cycle, DNA damage, Wnt, MAPK, and mTOR pathways (Figure 8(b)), which was consistent with previous reports showing that NPAS2 activated MAPK signals to be associated with poor survival of primary tumors [44,45]. CRY2 has also been demonstrated to activate the mTOR pathway to regulate the differentiation and function of immune cells [46]. ARNTL activated Wnt to promote tumor cells development [47]. On the contrary, CRY2 and DBP inhibited DNA damage and cell cycle pathways to enhance tumor cell growth. These results implied that abnormal expression of CCGs could affect cancer pathways to promote glioma progression. Besides, several genes involved in these pathways were also rhythmically expressed    (Figure 8(c)). Among closely coexpressed genes, the rhythmic expression of TEF was exactly the same as that of cluster II genes. TEF was controlled by the circa-dian rhythm and affected the expression of other rhythmic and functional genes [48], suggesting the modulation of TFE-mediated downstream effectors by CRY2 and DBP.

BioMed Research International
What is more, due to correlation of CRY2 and immune cells, immune infiltration was mainly contributed by CRY2, thus influencing the occurrence and progression of glioma.
In general, not all drugs are suitable for chronotherapy, most chronotherapy focuses on drugs with a half-life of less than 15 hours [49]. The plasma concentration of temozolomide peaks within 1 hour (T max < 1 h) after oral dosing and has a plasma half-life of 1.8 h, which makes the drug suitable for chronotherapy [49][50][51]. In fact, previous studies have shown that temozolomide is most effective on glioma cells at the ARNTL peak and least effective at the trough [38]. When considering the critical roles of pharmacokinetics and pharmacodynamics in drug metabolism and therapy, it is reasonable that administrating temozolomide at about 7-8:00 pm is likely to obtain a better curative effect. In our work, we showed the therapeutic relationship between temozolomide and ARNTL and pointed out a better time to administer temozolomide was about 7-8:00 pm (the peak expression of ARNTL) in the evening (Figure 11(b)), suggesting precise chronotherapy for glioma patients.

Conclusion
In summary, we provided a CCGs-based accurate prediction model and showed that ARNTL, NPAS2, CRY2, and DBP had great effects on diagnosis and prognosis of glioma. Besides, if the patient was sampled at night (means to have a higher expression of CCG), it can lead to early diagnosis of glioma, whereas sampling in the morning may cause misdiagnosis and thus delay the treatment of glioma. Furthermore, according to expression of ARNTL, appropriate timing of temozolomide administration can effectively improve the efficacy, which provides a reference for chronotherapy of glioma. Our works strengthen the importance of CCGs as prognosis biomarkers and the clinic implications for glioma chronotherapy. False discovery rate HR: Hazard ratio ROC: Receiver operating characteristic.