Circadian Regulation Patterns With Distinct Immune Landscapes in Gliomas Aid in the Development of a Risk Model to Predict Prognosis and Therapeutic Response

Circadian disruption in tumorigenesis has been extensively studied, but how circadian rhythm (CR) affects the formation of tumor microenvironment (TME) and the crosstalk between TME and cancer cells is largely unknown, especially in gliomas. Herein, we retrospectively analyzed transcriptome data and clinical parameters of glioma patients from public databases to explore circadian rhythm-controlled tumor heterogeneity and characteristics of TME in gliomas. Firstly, we pioneered the construction of a CR gene set collated from five datasets and review literatures. Unsupervised clustering was used to identify two CR clusters with different CR patterns on the basis of the expression of CR genes. Remarkably, the CR cluster-B was characterized by enriched myeloid cells and activated immune-related pathways. Next, we applied principal component analysis to construct a CRscore to quantify CR patterns of individual tumors, and the function of the CRscore in prognostic prediction was further verified by univariate and multivariate regression analyses in combination with a nomogram. The CRscore could not only be an independent factor to predict prognosis of glioma patients but also guide patients to choose suitable treatment strategies: immunotherapy or chemotherapy. A glioma patient with a high CRscore might respond to immune checkpoint blockade, whereas one with a low CRscore could benefit from chemotherapy. In this study, we revealed that circadian rhythms modulated tumor heterogeneity, TME diversity, and complexity in gliomas. Evaluating the CRscore of an individual tumor would contribute to gaining a greater understanding of the tumor immune status of each patient, enhancing the accuracy of prognostic prediction, and suggesting more effective treatment options.

Introduction

Gliomas represent 80% of all primary brain tumors and are a heterogeneous group of lethal tumors of the central nervous system (CNS) (1). In 2016, the newly established World Health Organization (WHO) classification of CNS tumors added molecular genetic features including isocitrate dehydrogenase (IDH) and chromosomal 1p/19q status to the previous version, which solely histologically classified gliomas into four grades (2). Despite very comprehensive treatment protocols consisting of debulking surgery, chemotherapy, and concomitant radiotherapy (3), the median survival of grade IV glioblastoma multiforme (GBM) is only 12–14 months (4). A remarkable level of genetic, epigenetic, and environmental heterogeneity existing within each individual glioma explains multiple mechanisms of therapeutic resistance and forms a highly resilient disease (5). A comprehensive model based on gliomas integrated with TME will be of benefit to predicting prognosis and selecting drugs and meet the demand of individualized treatment.

Cancer immunotherapies advancing from the trial stage to established first- or second-line indications have revolutionized the treatment of solid tumors. However, clinical trials of glioma immunotherapy are not proceeding well. The limited response of glioma to immunotherapy may be dependent on the tumor-intrinsic characteristics and the brain microenvironment. The specific TME of gliomas, characterized by the low frequency of T cells, the high frequency of myeloid cells, and the lower expression of cell-surface inhibitory markers, makes it challenging for immune checkpoint blockade (ICB) to exert a strong therapeutic effect (6–8). Meanwhile, a biomarker to optimize patient selection is most urgently needed due to a lack of broad ICB approval.

Circadian rhythm is a conserved phenomenon that governs a large array of physiological and metabolic functions. The mammalian circadian rhythm depends on a time-delayed transcription translation feedback loop (TTFL). The aryl hydrocarbon receptor nuclear translocator-like protein 1 (ARNTL; also named brain and muscle ARNT-like protein 1, BMAL1) and circadian locomotor output cycles kaput (CLOCK) constitute the positive arm (9, 10), which can promote the expression of cryptochrome (CRY1 and CRY2) and period (PER1, PER2, and PER3). As the negative arm of the TTFL, CRYs and PERs form a complex to suppress the BMAL1–CLOCK complex (11, 12). As for the second feedback loop, the BMAL1–CLOCK complex regulates the expression of nuclear receptors REV-ERBα/β (also known as nuclear receptor subfamily 1, group D, members 1/2, NR1D1/2) and retinoic acid receptor-related orphan receptors (RORs), which in turn repress and activate BMAL1, respectively (13, 14). Over the past several decades the connection between circadian clocks and tumorigenesis has been well studied. Abnormal rhythms are associated with high-grade brain tumors (15, 16). Specifically, disruption of the circadian clock pharmacologically (REV-ERB agonists) or genetically (CLOCK and BMAL1 short hairpin RNAs) impairs glioma stem cell (GSC) stemness in GBM and causes GSC cell-cycle arrest and apoptosis (17, 18). As for TME, high CLOCK levels in GSCs correlate with a high frequency of microglia via the regulation of olfactomedin-like 3 and correlate with a decreased level of activated CD8+ T cells (18). Therapeutically, in view of cyclical expression patterns of DNA repair genes (19), circadian rhythms influence therapeutic sensitivity and thus play prominent roles in regulating the antitumor efficiency of chemotherapy (20). These findings suggest targeting clock-regulated crosstalk between TME and cancer cells will help to develop a novel prognosis signature and effective clock-oriented immunotherapies that may have synergistic effects with conventional therapies to increase the effectiveness of glioma treatment.

This bioinformatics analysis is intended to explore circadian rhythm-controlled tumor heterogeneity and characteristics of TME in gliomas. The dataset of glioma in The Cancer Genome Atlas (TCGA) was used as the training set and that in the Chinese Glioma Genome Atlas (CGGA) was used as the validation set throughout the study. Unsupervised clustering identified two stable CR patterns with different characteristics of infiltrated immune cells, activated signaling pathways and immune response processes. PCA algorithm was used to construct the CRscore that could not only be an independent factor to predict the prognosis of glioma patients but also guide them to choose suitable treatment strategies: immunotherapy or chemotherapy. In this study, we revealed that glioma immune landscapes regulated by two distinct circadian patterns renewed our cognition to improving individualized therapy based on clock within gliomas.

Methods

Data Acquisition and Preprocessing

The RNA sequencing data of TCGA-LGG (Lower Grade Glioma; WHO grade II–III), TCGA-GBM (glioblastoma multiforme; WHO grade IV), and GTEx-brain in transcripts per kilobase million (TPM) format uniformly processed by TOIL were downloaded from UCSC xena (https://xenabrowser.net/datapage/). The corresponding clinicopathological data were obtained from the cBioPortal website (http://www.cbioportal.org/). Patients without survival information were excluded in a further study. The genomic mutation data (including somatic mutation and copy number variation) of TCGA-LGG and TCGA-GBM were downloaded from the Genomic Data Commons (https://portal.gdc.cancer.gov/) using the “TCGAbiolinks” R package (21). Somatic mutation data were analyzed using “maftools” R package, and significant amplifications or deletions of copy number were detected using GISTIC 2.0. To obtain the CGGA validation set, the RNA sequencing data and related clinicopathological data were downloaded from the CGGA website (https://www.cgga.org.cn). Two independent datasets, including one with 325 individuals and another with 693 individuals, were acquired. Batch effects from non-biological technical biases were corrected using the “ComBat” algorithm of the “sva” R package (22).

Unsupervised Clustering for Circadian Rhythm Genes

Six CR gene sets were concluded from the Molecular Signatures Database (MSigDB) (23), KEGG pathways (24), REACTOME database (25), GO biological processes (GO: 0007623 CIRCADIAN RHYTHM), WikiPathways (26), and review literatures (27–29). The gene annotated in not less than two CR gene sets was defined as a CR gene for further integrated analysis and modeling. A total of 91 acknowledged CR genes were curated and analyzed to identify distinct CR patterns (Table S1 and Figure S1). The CR gene widely described in reviews is defined as a core CR gene. Unsupervised clustering analysis was applied to identify CR patterns based on the expression of 91 CR genes from TCGA or CGGA and classify patients for further analysis. A consensus clustering algorithm was performed using the “ConsensuClusterPlus” R package and was repeated 1,000 times for guaranteeing the stability of classification (30).

Gene Set Variation Analysis and Gene Ontology Annotation

To investigate the difference in biological process between CR patterns, we performed GSVA enrichment analysis using the “GSVA” R package (31). The gene sets of “h.all.v7.4.symbols” and “c2.cp.kegg.v7.4.symbols” were downloaded from the MSigDB database to run GSVA analysis (http://www.gsea-msigdb.org/gsea/downloads.jsp). The GSVA scores among different circadian patterns were determined using the “limma” R package. The threshold was set as the adjusted p-value < 0.05 and the absolute (log₂ Fold change, FC) > 0.1. The “clusterProfiler” R package was used to perform functional annotation for CR-related genes, and the top 20 GO annotations were shown (32).

Estimation of TME Cell Infiltration and Signatures

As we did in our previous study, the “ESTIMATE” algorithm was used to calculate Tumor Purity, ESTIMATE Score, Immune Score, and Stromal Score (33, 34). Furthermore, we used the ssGSEA (Single-sample Gene Set Enrichment Analysis) algorithm, CIBERSORT algorithm (35), and xCELL algorithm (36) to quantify the relative abundance of each type of cells infiltrating the glioma TME. The gene sets of “TME-associated signatures” were downloaded and analyzed using the “IOBR” R package (37). The ssGSEA method was chosen in the process of signature score evaluation. TIP (Tracking Tumor Immunophenotype) is a meta-server that systematically integrates two existing third-party methods “ssGSEA” and “CIBERSORT” for tracking, analyzing, and visualizing the status of anti-cancer immunity and the proportion of tumor-infiltrating immune cells across a seven-step cancer-immunity cycle using RNA-seq or microarray data (38). The correlations between the CRscore and the steps of the cancer-immunity cycle were analyzed using the “ggcor” R package.

Generation of the CRscore

The construction of CRscore was performed as follows: firstly, the differentially expressed genes (DEGs) associated with the CR cluster phenotype were determined using the “limma” R package. Specifically, gene expression data were normalized by voom and then fed to lmFit and eBayes functions to calculate the differentially expressed statistics. The significance filtering criteria of DEGs were set as an adjusted p-value < 0.05 and |log₂(FC)| ≥1.5. Then, the prognostic analysis was performed for each DEGs using a univariate Cox regression model. A total of 719 genes with significant prognosis were extracted for further analysis using the “randomForestSRC” R package. We also used the “randomSurvivalForest” algorithm to rank the importance of 719 genes, and the gene with a relative importance >0.2 was selected to construct the CR signature (39). We then curated the expression profile of the final 50 determined genes to perform PCA, and principal components 1 and 2 were extracted and served as the signature score. We then adopted a formula similar to previous studies (40), CRscore = Σ(PC1_i + PC2_i), where i is the expression of 50 prognostic-related CR cluster DEGs.

Immuno-/Chemotherapeutic Response Prediction

The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (41) and Immune Cell Abundance Identifier (ImmuCellAI) algorithm (42) were used to predict response to ICB therapy. We also used immunotherapeutic cohorts with complete clinical information to predict patients’ response to ICB therapy. In the GSE78220 cohort, the data of patients with metastatic melanoma treated with pembrolizumab, an anti-PD-1 antibody, were downloaded from Gene Expression Omnibus (GEO), and the FPKM data of gene expression profiles were converted to the more comparable TPM value among samples. We also used the “Prophetic” R package to predict the chemotherapy response of each sample based on Genomics of Drug Sensitivity in Cancer (GDSC) (https://www.cancerrxgene.org/). The correlation between estimated IC₅₀ and CRscore was established using Spearman correlation analysis with thresholds as |Spearman R| > 0.5 and p-value < 0.05.

Statistical Analyses

For comparison between two groups, statistical significance for normally distributed variables was estimated by Student’s t-tests, and nonnormally distributed variables were analyzed by the Wilcoxon rank-sum test. For comparison among three or more groups, one-way ANOVA and Kruskal–Wallis test were used for normally or nonnormally distributed variables, respectively. Fisher’s exact test was performed for categorical data. Correlation coefficients were computed by Spearman and distance correlation analyses. The “surv-cutpoint” function from the “survival” R package was applied to stratify samples into CR-high and -low groups. Kaplan–Meier curves and the Log-Rank test were adopted to assess whether there were differences in overall survival among groups. Univariate and multivariate Cox regression analyses were utilized to evaluate the independent prognostic value of the CRscore regarding OS, which were revealed by forest maps. We next performed multivariate Cox regression to establish a nomogram; the survival predictive accuracy of prognostic models was assessed based on the calibration plots. All statistical p-values were two-sided, with p < 0.05 considered as statistically significant. All statistical tests were performed in R statistical software (v4.05, R Core Team, R Foundation for Statistical Computing, Vienna, Austria).

Results

Landscape of Genomic Variations and Transcription Profile of Circadian Rhythm Genes in Glioma

To investigate whether there were specific expression patterns of CR genes in gliomas, we determined several differentially expressed CR genes between glioma and normal samples using the “limma” R package (Figure 1A). The genes involved in negative transcription regulation (ID3, ID2, and KLF10), cell metabolism (NACLU and AHCY), and cell survival (MAGED1 and TIMELESS) were upregulated in tumors when |log₂FC| ≥ 1 and adj.p ≥ 0.05 were set as thresholds. Furthermore, the comprehensive landscape of interaction among CR genes and feature groups to which they belong was depicted in the Figure 1B network, and we noticed that most of the CR genes were positively correlated with each other. To investigate whether genomic variation contributed to abnormal CR gene expression in gliomas, we described mutations and copy number variations (CNVs) of CR genes, respectively. In the top 10 mutations of CR genes, only PTEN exerted 11% mutant frequency, suggesting that genetic mutations of CR genes might not be the main cause of rhythm perturbations (Figure 1C). The vast majority were mutations in PETN, a common tumor suppressor gene, which could be included here because of its expression governed by circadian rhythms. Besides, the CNV frequency showed that CLOCK, the positive arm of TTFL, gained copy number, while PER3, a part of negative arm of TTFL, had a frequency of CNV deletion (Figure 1D). Thus, we have come to know that CNV partly conduced to disturbed rhythms of tumors. The location of CR genes with CNVs on chromosomes is shown in Figure 1E. To explore whether the core regulators of circadian TTFLs had unique expression patterns in gliomas, we analyzed the mRNA levels of core CR genes among normal, LGG, and GBM samples (Figure 1F). It was worth noting that major factors (CRYs and PERs) of negative arms were upregulated in LGG compared to normal, but with relatively low expression in GBM. Given all this, the disruption of circadian rhythms due to the high heterogeneity of genomic variations and transcription profile covertly contributed to tumorigenesis and progression of gliomas.

FIGURE 1

Figure 1 Landscape of genomic variations and transcription profile of circadian rhythm genes in gliomas. (A) Volcano plot of differentially expressed CR genes between normal and glioma (TCGA-LGGGBM: TCGA-LGG and TCGA-GBM) samples (Wilcoxon test: adjust p < 0.05, and |log2FC| > 1). The overexpressed genes in gliomas were highlighted in red. (B) The interaction between CR genes in gliomas. The CR genes with different features were depicted by circles in different colors. The lines connecting CR genes represented their interaction with each other. The size of each circle represented the differential expression of CR genes in gliomas compared to normal. (C) The mutation frequency of CR genes in 584 patients from the TCGA cohort. Each column represented an individual patient. The upper bar plot indicated mutation-accumulation. The bar plot on the right indicated the proportion of each variant type with the number above representing mutation frequency. Only top 10 mutations were included. (D) The CNV frequency of CR genes in glioma patients from the TCGA cohort. The height of each column represented the alteration frequency. The amplification frequency, pink dot; The deletion frequency, blue dot. (E) The location of CR genes with CNV alteration on chromosomes. (F) The mRNA expression level of core CR genes among normal, TCGA-LGG, and TCGA-GBM samples. Normal, blue; LGG, pink; GBM, red. The asterisks represented the statistical p-value (Kruskal–Wallis test: *p < 0.05; ***p < 0.001).

Identification of Circadian Patterns Mediated by CR Genes

Given the aberrant expression of certain circadian regulators in gliomas, we wondered whether there were distinctively different CR patterns, on the whole, among individual tumors. The unsupervised hierarchical cluster analysis was used to uncover circadian patterns, and the TCGA cohort population was separated into two distant clusters, termed CR cluster-A and -B, by the “ConsensusClusterPlus” R package (Figures 2A, S2A–C). Additionally, similar results obtained in the CGGA cohort verified our findings above (Figures S3A–C). Moreover, principal component analysis confirmed that CR genes could distinguish two CR clusters perfectly, in both TCGA and CGGA cohorts (Figure 2B). Obvious discrepancy between two CR clusters was embodied in histopathology features and patients’ prognoses. Patients in CR cluster-A were characterized by lower grades and harmless molecular genetic features, which coincided with particularly prominent survival advantage (Figures 2C, S2D, E). The same results were obtained in LGG and GBM groups, respectively (Figures S2D, S3D). We also added the unsupervised clustering results in LGG and GBM, respectively. Unsupervised clustering could also divide LGG patients into two groups, and there was a significant difference in prognosis between the two groups (Figures S2F–I). As for core CR genes, the distribution of their transcriptome levels with remarkable differences between two clusters corresponded with the previous clustering (Figure 2D). ARNTL (BMAL1) and CLOCK together constitute the positive arm of TTFLs, but they took advantages of expression level in different patterns. Similarly, the components of the negative arm were expressed preferentially.

FIGURE 2

Figure 2 Identification of two circadian patterns mediated by CR genes in TCGA cohort. (A) Unsupervised clustering of 91 CR genes for 596 glioma patients in the TCGA cohort resulted in two CR clusters. Age, gender, tissue, WHO grade, histopathology, IDH status, ATRX status, MGMT promoter status, TERT promoter status, and survival status are shown as patient annotations. Pink represented the relatively high expression of CR genes and blue represented the relatively low expression. (B) Principal component analyses for the transcriptome profiles of CR patterns in TCGA (left) and CGGA (right) cohorts, respectively, showing a remarkable difference on transcriptome between two CR patterns. (C) Survival analyses for CR patterns in the TCGA cohort (left) including 379 cases in CR cluster-A and 217 cases in CR cluster-B (Log-Rank test: p < 0.0001), and in the CGGA cohort (right) including 614 cases in CR cluster-A and 356 cases in CR cluster-B using Kaplan–Meier curves (Log-Rank test: p < 0.0001). (D) The transcriptome of core CR genes with remarkable differences in two CR patterns corresponded with the previous clustering. Pink represented the relatively high expression of core CR genes and blue represented the relatively low expression. (E, F) GSVA enrichment analyses in two CR clusters showing the activation states of Hallmark pathways (MSigDB) in TCGA (E) and CGGA (F) cohorts. Activated pathways, pink; Inhibited pathways, blue.

Next, to explore the global functions in two CR patterns beyond the individual and single gene expression, we performed GSVA enrichment analyses, which showed the activation states of Hallmark pathways (MSigDB) in TCGA (Figure 2E) and CGGA (Figure 2F) cohorts. Hedgehog Signaling activation, Pancreas beta cells, and Kras Signaling Downregulation were observed in CR cluster-A, while immune response-related pathways consisting of “Interferon Response”, “Interleukin Signaling”, and “Tnfα Signaling via Nfκb” and pathways associated with malignant biological behaviors including “Epithelial-Mesenchymal Transition” and “Angiogenesis” were activated in CR cluster-B. KEGG pathway enrichment analyses in TCGA (Figure S3F) and CGGA (Figure S3G) cohorts also supported the activated immune reaction-related pathways. From above, we identified two internally gathered circadian patterns and functional enrichment analyses recognized CR cluster-B closely related to the immune response.

Characterization of the Immune Cell Infiltration in Distinct Circadian Patterns

To describe immune status in different CR patterns, we utilized four different methods (ESTIMATE, ssGSEA, CIBERSORT, and xCELL) for evaluating immune cell infiltration (Figure 3A). ESTIMATE algorithm showed that CR cluster-B exhibited high Stromal and Immune Scores and low Tumor Purity, which represented a significantly increased immune cell infiltration and generally indicated a relatively hot tumor immune microenvironment. However, in gliomas, components of immune system enriched in CR cluster-B were macrophages (ssGSEA, CIBERSORT, and xCELL), dendritic cells (ssGSEA and xCELL), neutrophils (ssGSEA, CIBERSORT, and xCELL), endothelial cells (xCELL), and Th2 cells (ssGSEA and xCELL), generally considered as immunosuppressive components in tumors. Additionally, CR cluster-B was mainly characterized by relatively high expressions of human leukocyte antigens (HLA) and immune checkpoint molecules (Figure 3B). On the one hand, class I HLA, expressed on the surface of almost all nucleated cells, could increase the possibility of presenting more immunogenic antigens and of benefiting from ICB. High immune checkpoint expression also implied effective anti-immune checkpoint therapy. On the other hand, in line with the characteristics of immune cell infiltration, class II HLA molecules, mainly expressed on the surface of antigen-presenting cells (dendritic cells, B cells and macrophages), and many inhibitory molecules (IL-4, IL-10, and VEGF) were unregulated in CR cluster-B. We then examined the specific correlation between each TME infiltration cell type and each core CR gene using Spearman’s correlation analyses (Figure 3C). The correlation focused on the relationship between immunosuppressive cells (macrophages, neutrophils, and T helper cells) and the core TTFL components (ARNTL, BHLHE40/41, CRY2, and PER2/3). Furthermore, CR clusters could be distinguished by different TME relevant signatures and TME score (Figure 3D). The above results revealed that immunosuppressive components enriched in CR cluster-B suggested an immunosuppressive microenvironment and coincided with poor prognoses of patients.

FIGURE 3

Figure 3 Characterization of the immune cell infiltration in distinct circadian patterns. (A) Heatmap of TME cell infiltration characteristics in two CR patterns assessed by four different methods. (B) Heatmap showing the differences of immune-related gene expression in two CR patterns. Upregulation, pink; Downregulation, blue. (C) The correlation between TME infiltration cell type and each core CR genes using spearman analyses. Negative correlation, blue; Positive correlation, pink. (D) CR clusters were distinguished by different TME relevant signatures and TME score (Wilcoxon test: ***p < 0.001).

Generation of Circadian Gene Clusters and Functional Annotation

To further validate CR regulation patterns and investigate the potential biological behavior of CR patterns, we determined 719 genes, which were DEGs between CR cluster-A and -B and also related to prognosis of glioma patients using the “limma” R package and univariate Cox regression analysis (Figure 4A). Then, we performed unsupervised clustering analyses based on the obtained CR-related DEGs. Consistent with the clustering grouping of previous CR patterns, the unsupervised clustering algorithm also revealed two distinct CR-related signatures termed CR gene cluster-A and -B, respectively (Figure 4B, Figures S4A, B). This confirmed that two distinct CR patterns did exist in gliomas. Survival analysis for CR gene clusters in the TCGA cohort using Kaplan–Meier curves suggested survival advantages in CR gene cluster-A no matter in all glioma samples (Figure 4C) or in TCGA-LGG and TCGA-GBM, separately (Figure S4C). Next, the “clusterProfiler” R package was used to perform GO enrichment analysis for Gene Set 1 that were upregulated in gene cluster-B (Figure 4D), and Gene Set 2 upregulated in gene cluster-A (Figure 4E). GO’s annotations for Gene Set 1 fell into several categories including antigen presentation, inflammation, and immune responses, while the annotations for Gene Set 2 focused on neurotransmitter transmission and other major duties of CR genes controlling the brain. These analyses confirmed that gene cluster-B was significantly associated with immune-relevant signatures, whereas gene cluster-A was associated with regulations of synaptic transmission. Besides, the expression of core CR genes (Figure S4D), the infiltration of immune cells (Figure S4E), and Hallmark (Figure S4F) as well as KEGG (Figure S4G) pathways in two CR gene clusters were consistent with the previous clustering results. The results of secondary clustering verified that there were indeed two CR-based patterns in gliomas.

FIGURE 4

Figure 4 Generation of circadian gene clusters and functional annotations for CR-related genes. (A) Volcano plot of CR pattern-related DEGs between CR cluster-A and -B (Wilcoxon test: adjust p < 0.05, and |log2FC| > 1.5). (B) Unsupervised clustering of CR pattern-related DEGs to classify patients of the TCGA cohort into different subtypes, termed CR gene cluster-A and -B, respectively. The gene clusters, CR clusters, and other parameters were used as patient annotations. (C) Survival analysis for CR gene clusters in the TCGA cohort including 368 cases in gene cluster-A and 228 cases in gene cluster-B using Kaplan–Meier curves (Log-Rank test, p < 0.0001). (D) Functional annotation for CR-related genes upregulated in gene cluster-B (Gene Set 1 in B) using GO enrichment analysis. The length of the bar plots represents the number of genes in that category. The color depth represented q-value. (E) Functional annotation for CR-related genes upregulated in gene cluster-A (Gene Set 2 in B) using GO enrichment analysis. The length of the bar plots represents the number of genes in that category. The color depth represented q-value.

Construction of the CRscore and Exploration of Its Immunological Relevance

Considering the intratumoral heterogeneity and complexity of circadian disruption, based on 50 prognostic-related CR cluster DEGs, we developed a score termed CRscore using the PCA algorithm. The general distribution of CRscore, CR clusters, gene clusters, and molecular pathological parameters in the TCGA cohort were described in the heatmap (Figure 5A). Low CRscore crowded in CR cluster-B and gene cluster-B (Figure 5B). Considering that the immune microenvironment of LGG and GBM might differ, we explored the immunological relevance of CRscore separately in LGG and GBM. The correlation was supported by TME-related pathways and immune cell infiltration analyses and even more pronounced in LGG (Figures 5C, D). Next, we explored the indicating effect of CRscore in the anti-tumor immune process (Figure 5E) and discovered that CRscore was likely to be associated with the recruitment of immune cells. Some steps including priming and activation of immune cells; recruiting Th1, Th22 cells, and macrophages; and tumor killing effect were positively correlated with CRscore in both LGG and GBM. Peculiarly, in steps significantly correlated with the CRscore, a low score in GBM is related to recruitment of immunosuppressive cells (neutrophils, eosinophils, basophils, Tregs, and MDSCs). A low score in LGG indicated enhanced recruitment of anti-tumor immune cells (CD8+ cells, Th1 cells, and DCs).

FIGURE 5

Figure 5 Construction of the CRscore and exploration of its immunological relevance. (A) An overview of the association between CRscores and other patient annotations in TCGA cohort. (B) Comparison of CRscores between gene cluster A and B in the TCGA cohort, p < 0.001 (left); Comparison of CRscores between CR cluster A and B in the TCGA cohort, p < 0.001 (right). (C) The correlations between CRscore and TME relevant signatures in LGG as well as GBM of TCGA cohort, respectively (p > 0.05). Pink represented positive correlation and blue represented negative correlation in LGG; purple represented positive correlation and green represented negative correlation in GBM. (D) The correlations between CRscore and ssGSEA scores of TME cells in LGG as well as GBM of the TCGA cohort, respectively (p > 0.05). Pink represented positive correlation and blue represented negative correlation in LGG; purple represented positive correlation and green represented negative correlation in GBM. (E) The correlations between CRscore and steps of the cancer immunity cycle in LGG (right) as well as GBM (left) of TCGA cohort, respectively.

Evaluation of the Prognostic Potentiality of the CRscore

After having identified the CRscore as an intrinsic pattern closely linked to the immune process, we sought to determine whether the CRscore could accurately predict outcomes of glioma patients. Firstly, the overlap among tissue types, CR clusters, gene clusters, and CR groups was described in Figure 6A. The 596 patients in the TCGA cohort were assigned to two groups (CR-high and -low) based on CRscores using the cutoff value (0.4) obtained with the “survminer” R package. Survival analysis indicated that the CR-high group had prolonged survival time in the TCGA cohort as well as in LGG or GBM (Figure S5A), which was further validated in the CGGA cohort (Figures 6B, S5B). The area under the ROC curves (AUC) at 5 years (Figure 6C) and time-dependent AUC (Figure 6D) of the CRscore compared with other histological or molecular indicators in the TCGA and CGGA cohort, respectively, verified the predictive efficiency of the CRscore. When the CRscore was evaluated as a categorical variable (high or low CRscore) in the Cox regression model, the CRscore was determined to be an independent and robust prognostic factor. Univariate and multivariate analyses of clinicopathological characteristics and CRscore in the TCGA and CGGA cohort are shown in Figures 6E, F, respectively. Furthermore, we established a prognostic nomogram to predict 3-, 5-, and 10-year overall survival based on the stepwise Cox regression model. WHO grade, CRscore, and age were included in the prediction model (Figure 6G). The C-index of the nomogram was 0.864 (95% CI, 0.840–0.888). Nomogram prediction and actual observation in the TCGA cohort reached an excellent agreement at the 3-, 5-, and 10-year survival probability after calibration (Figure 6H). Net decision curve and the net reduction analyses demonstrated the superiority of this nomogram in predicting prognosis (Figures 6I, J).

FIGURE 6

Figure 6 Evaluation of the prognostic potentiality of the CRscore. (A) The overlap among tissue types, CR clusters, gene clusters, and CR groups. (B) Survival analyses for CR groups in the TCGA cohort (left) including 357 cases in the CR-high group and 239 cases in the CR-low group (Log-Rank test, p < 0.0001), and in the CGGA (right) cohort including 507 cases in the CR-high group and 463 cases in the CR-low group using Kaplan–Meier curves (Log-Rank test, p < 0.0001). (C) Predictive accuracy at 5-year of CRscore compared with other histological or molecular indicators in TCGA (left) and CGGA (right) cohort, respectively. The accuracy was equal to the area under the ROC curves (AUC). (D) Time-dependent AUC of CRscore compared with other histological or molecular indicators in the TCGA (left) and CGGA (right) cohort, respectively. (E, F) Univariate and multivariate analyses of clinicopathological characteristics and CRscore with overall survival in the TCGA (E) and CGGA (F) cohort, respectively. (G) Nomograms for predicting the probability of patient mortality based on WHO Grade, CRscore, and age. (H) Plots depicted the calibration of nomograms. (I) Net decision curve analyses demonstrating the benefit for predicting overall survival on nomogram. (J) The net reduction analyses demonstrated in how many patients the nomogram could be avoided without prognosis of miscalculation.

Landscape of Tumor Somatic Mutation in Two CR Groups

Given that the aim of exploring the value of the CR pattern on treatment and the response to ICB are closely related to somatic mutation, we decided to analyze differences in somatic mutation between CR-high and -low groups in the TCGA cohort using the “maftools” R package. As shown in Figure 7A, the CR-low group presented more extensive tumor mutation burden (TMB), on the whole, than the high-score group. Differences in CRscores between wild-type and mutant groups of each gene inversely proved the negative correlation between mutation and score (Figures 7B, C). The mutation occurrence varied between two groups. In detail, genetic mutations of IDH1, TP53, ATRX, and CIC mainly appeared in the CR-high group, while TTN, PTEN, and EGFR appeared in the CR-low group (Figure 7D).

FIGURE 7

Figure 7 Landscape of tumor somatic mutation in two CR groups. (A) Differences of TMB in CR-high and -low groups in TCGA cohort (Wilcoxon test: p < 0.001). (B) Top 20 most frequent mutations in patients from TCGA cohort and the distribution of mutations in CR high and low groups. (C) Differences in CR scores between wild-type and mutant groups of each of the top 20 genes in the TCGA cohort (Wilcoxon test: ns, P>0.05; *p < 0.05; **p < 0.01; ***p < 0.001). (D) Differences in the frequency of top 20 mutations between CR-high and -low groups in the TCGA cohort (Fisher’ exact test: p > 0.05; *p < 0.05; **p < 0.01; ***p < 0.001).

The Role of the CRscore in the Prediction of Therapeutic Benefits

Newly identified predictors, such as TIDE (Figure 8A) and ImmuCellA (Figure 8B), are widely used to evaluate the immune response. Our analysis revealed that a lower CRscore was not only associated with a poorer response to ICB but also more prone to immune escape. In anti-PD-1 immunotherapy cohort (GSE78220), a survival benefit trend was observed in patients with high CRscore (Figures 8C–E). To further understand the effects of the CRscore on drug response, we assessed the association between the CRscore and the response to drugs in cancer cell lines. Using the Spearman correlation analysis, we identified twenty-two drugs targeting various pathways (Figure 8G) significantly correlated between CRscores in the Genomics of Drug Sensitivity in Cancer (GDSC) database (Figure 8F). Among them, twenty agents with lower IC₅₀ in the CR-low group might be beneficial to low CRscore patients, while two agents might be beneficial to patients with high scores (Figure 8H). Together, these results implied that CR pattern played a crucial role in mediating the immune response and was also correlated with drug sensitivity. Thus, the CRscore might be a potential biomarker for establishing appropriate treatment strategies.

FIGURE 8

Figure 8 The role of the CRscore in the prediction of therapeutic benefits. (A, B) The relative distribution of TIDE (A) in TCGA (p > 0.001) and CGGA (p > 0.001) cohorts and ImmuCellAI (B) in TCGA (p = 0.17) and CGGA (p > 0.001) cohorts were compared between CR score high versus low groups, respectively. (C) The proportion of patients’ response to PD-1 blockade immunotherapy in CR-high or -low groups. PD, progressive disease; CR, complete response; PR, partial response (Fisher’ exact test: p = 0.013). (D) Survival analyses for low (8 cases) and high (19 cases) CR score patient groups in the anti-PD-1 immunotherapy cohort (GSE78220) using Kaplan–Meier curves (Log-Rank test: p = 0.0066). (E) The correlation of CRscore with clinical response to anti-PD-1 immunotherapy. (F) The correlation between CRscore and the estimated IC₅₀ for drugs evaluated by the Spearman analysis. Each point represents a drug. (G) Pathways targeted by candidate drugs. (H) The box plots of the estimated IC₅₀ for candidate drugs in CR-low and -high groups (Wilcoxon test: ***p < 0.001).

Discussion

Circadian disruption is associated with tumorigenesis and tumor progression through effects on cancer cell biological properties, including proliferation, DNA repair, apoptosis, metabolism, and stemness. Emerging evidence indicates that circadian clock also plays an influential role in the TME. Abnormalities of clock components (e.g., CLOCK, ARNTL/BMAL1, PER, CRY, RORs, and REV-ERBα) in tumors could reshape TME by regulating the expression and secretion of soluble factors including HIF1α, ARNT, VEGF, OLFML3, and other unidentified factors (27). These factors modulated TME biology, including endothelial cell biology such as promoting angiogenesis and antiangiogenic therapy resistance, the infiltration of myeloid cells, as well as the infiltration and activation/suppression of lymphocytes. Clock components in immune cells, in turn, affected tumor growth. BMAL1 in macrophages inhibited the production of ROS and HIF1α, and affected tumor growth through regulating macrophage alternative polarization (43) RORα and RORγ in T cells (Th17 cells and CD8+ T cells) could modulate their differentiation and activation, which affected tumor growth and the antitumor immune response (44–46). The above studies pointed out the importance and necessity of studying the clock-regulated crosstalk between TME and cancer cells.

On the basis of the CR gene list concluded from six gene sets, we performed the first unsupervised clustering in glioma patients from TCGA and CGGA cohorts, respectively, and identified two CR clusters with different composition and expression of core CR genes, which were considered as different CR patterns. The CR cluster-B showed a unique immune status characterized by high infiltration of immune cells, high expression of immune-related genes, and highly activated immune-related pathways, which suggested that circadian rhythms might affect the biological behavior of tumor itself, as well as the infiltration and function of immune cells, thus affecting the prognosis. Considering that the pro-tumoral pathway was activated and the immunosuppressive cells were infiltrated, it was reasonable that the outcomes of the patient were relatively poorer in CR cluster-B. CSNK1E, shown to phosphorylate PERs, and PERs were all downregulated in the CR cluster-B group, which were all negatively correlated with the infiltration of most immune cells. BHLHE40/41 repressing CLOCK-ARNTL’s transactivation of PER1 corresponds to the opposite expression of BHLHE40/41 and PER1 at mRNA levels. In view of the correlation between RORγ (RORC) and Th17/T cells has been confirmed by experimental studies (45, 46), correlations among genes or between a single gene and a certain type of immune cell might uncover the mechanism of circadian rhythms regulating tumor immune microenvironment, which needed to be explored further.

To construct a CR-based prognostic signature, a second unsupervised clustering was performed based on the DEGs between CR cluster-A and -B. Two distinct CR gene clusters (CR gene cluster-A and -B) proved that two stable CR patterns did exist in gliomas. Next, we constructed CRscore by using PCA algorithm and established the relationship between CRscore and immune signatures. Analyses of immune components and anti-tumor processes provided the reason why prognoses of patients in low score group were poor. CRscore was negatively associated with recruitment and infiltration of immune cells, as well as the pathways related to TME. However, there were some subtle and critical differences between LGG and GBM. In LGG, a lower CRscore was related to enhanced recruitment of immune cells but diminished killing tumor effect, characterized by immune disability. In GBM, a lower CRscore represented recruitment of immunosuppressive cells, characterized by immune suppression. Therefore, in both LGG and GBM, CRscore was negatively correlated with prognosis. More accurate prognostic evaluation method was established by univariate and multivariate regression analysis as well as nomogram and validated.

As neo-antigens result from mutations, the more neo-antigens that are present, the higher the TMB. With more neo-antigens present, there may be an increased probability that some of the neo-antigens presented by MHC proteins will be immunogenic (47, 48). This is the root of the hypothesis that cancers with high TMB are more likely to benefit from ICBs. Yet, approximately 60% of patients with high TMB do not respond to ICBs (49). Similarly, in this study, the CR-low group with more extensive TMB did not gain advantages in ICB treatment, and even the opposite result was obtained. Our results confirmed the theory that TMB had limitations as a predictive biomarker, especially when used in isolation (50). As mentioned above, patients with lower CRscores not only tended not to respond to ICBs but also were more prone to immune escape when using TIDE and ImmuCellA to evaluate the immune response. In other cancer cohorts treated with ICB (GSE78220), the CR-low group showed a higher progression disease ratio. However, in drug sensitivity analysis, the CR-low group was sensitive to more drugs, suggesting another treatment strategy. Recent advances of the glioma TME highlight the complex and immunosuppressive environment within the tumor that underlies the resistance to immunotherapy. According to the characteristics of glioma, ICB strategies different from other solid tumors need to be developed urgently. Macrophages found within the glioma tend to be immunosuppressive and are associated with poor survival, suggesting that tumor-associated macrophage or microglia (TAM)-targeting therapies (alone or in combination with other therapies) altering TAM characteristics or abundance may improve immunotherapeutic outcomes (6, 51). Strategies to enhance therapy response to ICB might thus involve the mechanism-driven combination of ICB and targeting of TAMs.

Although we reviewed the databases and literatures to curate a catalog of 91 recognized CR genes, a series of newly identified or side CR genes need to be incorporated into the model to optimize the accuracy of the CRscore. Besides, the CR patterns and the CRscore were identified by using retrospective datasets; thus, the role of CR score in immunotherapy needs to be confirmed by more indicators, and more prospective cohorts of glioma patients receiving immunotherapy are needed to validate our findings.

In conclusion, we comprehensively evaluated and verified the CR patterns among 1,566 glioma samples based on recognized CR genes, and systematically correlated these patterns with TME characteristics. This integrated analysis indicated that circadian rhythm upset laid a critical foundation for understanding crosstalk between TME and tumor cells in gliomas. More broadly, evaluating the CRscore of the individual tumor will contribute to enhancing our cognition of the characteristics of TME and provide important insight into therapeutic efficacy.

Data Availability Statement

Publicly available datasets were analyzed in this study. These data can be found here: TCGA (The Cancer Genome Atlas); CGGA ( Chinese Glioma Genome Atlas).

Author Contributions

MS, RT, and YL contributed to the study conception. RT was responsible for writing the first draft of the paper. RT and YL conducted the data analysis. MS was responsible for funding acquisition and supervision. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (Project 81872467 and 82073880), the Shanghai Municipal Health Commission (Program 201940233), and the Program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning (Program TP2019047).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank all members in the lab for critical reading and discussion.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.797450/full#supplementary-material

Supplementary Figure 1 | Six CR gene sets were concluded from five datasets and review literatures. The gene annotated in not less than two circadian gene sets was defined as a CR gene. A total of 91 acknowledged CR genes were concluded here. The number of pick circles connected by lines indicated the number of times a CR gene appears in different gene sets; The number of genes sharing the same annotation pattern was marked at the top of the corresponding bar chart. The bar chart on the left showed the number of circadian genes in each gene set.

Supplementary Figure 2 | Identification of two circadian patterns mediated by CR genes in TCGA cohort. (A) CDF plot showing a real random variable of its probability distribution based on consensus scores for each k (from 2 to 6, indicated by different colors) in TCGA cohort. (B) Consensus matrices of TCGA cohort for each k (k = 2–6), displaying the clustering stability using 1000 iterations of hierarchical clustering. (C) Survival analyses for CR patterns in TCGA cohort using Kaplan–Meier curves for OS when k = 3–6 (Log-Rank test: P < 0.0001). (D) Kaplan–Meier curves for OS of TCGA-LGG cohort (upper) including 372 cases in CR cluster-A and 72 cases in CR cluster-B (Log-Rank test: P < 0.0001), and TCGA-GBM (lower) including 7 cases in CR cluster-A and 145 cases in CR cluster-B (Log-Rank test: P = 0.0085). (E) The proportion of representative clinical parameters of patients in CR cluster A and B relevant to Figure 2A (Fisher’ exact test: P <0.001). (F) Consensus matrices of TCGA-LGG cohort for k = 2, displaying the clustering stability using 1000 iterations of hierarchical clustering. (G) Kaplan–Meier curves for OS of patients in TCGA-LGG cohort when k = 2. (H) Consensus matrices of TCGA-GBM cohort for k = 2, displaying the clustering stability using 1000 iterations of hierarchical clustering. (I) Kaplan–Meier curves for OS of patients in TCGA-GBM cohort when k = 2.

Supplementary Figure 3 | Verification of two circadian patterns in CGGA cohort and functional annotations of CR clusters. (A) CDF plot showing a real random variable of its probability distribution based on consensus scores for each k (from 2 to 6, indicated by different colors) in CGGA cohort. (B) Consensus matrices of CGGA cohort for each k (k = 2–6). (C) Kaplan–Meier curves for OS of CGGA cohort when k = 3–6 (Log-Rank test: P < 0.0001). (D) Kaplan–Meier curves for OS of LGG patients in CGGA cohort (upper) including 476 cases in CR cluster-A and 120 cases in CR cluster-B (Log-Rank test: P < 0.0001), and GBM patients in CGGA cohort (lower) including 138 cases in CR cluster-A and 236 cases in CR cluster-B (Log-Rank test: P < 0.0001). (E) Unsupervised clustering of 91 CR genes for 970 glioma patients in CGGA cohort resulted in two CR clusters. Age, gender, tissue, WHO grade, IDH status, 1p/19q codeletion status, MGMT promoter status, and survival status are shown as patient annotations. (F, G) GSVA enrichment analyses in TCGA (F) and CGGA (G) cohorts showing the activation states of KEGG pathways in two CR clusters. Activated pathways, pink; Inhibited pathways, blue.

Supplementary Figure 4 | Generation of circadian gene clusters and functional annotations of CR gene clusters. (A) Consensus matrices of TCGA cohort for k = 2. (B) CDF plot showing a real random variable of its probability distribution based on consensus scores for each k (from 2 to 6, indicated by different colors) in TCGA cohort. (C) Kaplan–Meier curves for OS of TCGA-LGG patients (left) and TCGA-GBM patients (right), respectively. (D) The mRNA expression level of core CR genes in gene cluster A and B. Gene cluster A, blue; Gene cluster B, pink. The asterisks represented the statistical p value (Wilcoxon test: ns, P >0.05; *, P <0.05; ***, P <0.001). (E) The ssGSEA score of TME cells in two CR gene clusters of TCGA cohort. (Wilcoxon test: ns, P>0.05; **, P <0.01; ***, P <0.001). (F, G) GSVA enrichment analyses in TCGA cohort showing the activation states of Hallmark pathways (MSigDB) (F) and KEGG pathways (G) in two CR gene clusters. Activated pathways, pink; Inhibited pathways, blue.

Supplementary Figure 5 | Survival analyses for CR groups. (A) Survival analyses for CR groups in TCGA-LGG and TCGA-GBM using Kaplan-Meier Curves (Log-Rank test, P <0.0001 in LGG; P <0.0001 in GBM). (B) Survival analyses for CR groups in LGG and GBM from CGGA cohort using Kaplan-Meier Curves (Log-Rank test, P <0.0001 in LGG; P =0.0068 in GBM).

References

1. Ostrom QT, Patil N, Cioffi G, Waite K, Kruchko C, Barnholtz-Sloan JS. CBTRUS Statistical Report: Primary Brain and Other Central Nervous System Tumors Diagnosed in the United States in 2013-2017. Neuro Oncol (2020) 22(12 Suppl 2):iv1–96. doi: 10.1093/neuonc/noaa200

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Louis DN, Perry A, Reifenberger G, von Deimling A, Figarella-Branger D, Cavenee WK, et al. The 2016 World Health Organization Classification of Tumors of the Central Nervous System: A Summary. Acta Neuropathol (2016) 131(6):803–20. doi: 10.1007/s00401-016-1545-1

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Bolcaen J, Kleynhans J, Nair S, Verhoeven J, Goethals I, Sathekge M, et al. A Perspective on the Radiopharmaceutical Requirements for Imaging and Therapy of Glioblastoma. Theranostics (2021) 11(16):7911–47. doi: 10.7150/thno.56639

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Touat M, Idbaih A, Sanson M, Ligon KL. Glioblastoma Targeted Therapy: Updated Approaches From Recent Biological Insights. Ann Oncol (2017) 28(7):1457–72. doi: 10.1093/annonc/mdx106

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Nicholson JG, Fine HA. Diffuse Glioma Heterogeneity and Its Therapeutic Implications. Cancer Discov (2021) 11(3):575–90. doi: 10.1158/2159-8290.CD-20-1474

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Goswami S, Walle T, Cornish AE, Basu S, Anandhan S, Fernandez I, et al. Immune Profiling of Human Tumors Identifies CD73 as a Combinatorial Target in Glioblastoma. Nat Med (2020) 26(1):39–46. doi: 10.1038/s41591-019-0694-x

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Iglesia MD, Parker JS, Hoadley KA, Serody JS, Perou CM, Vincent BG. Genomic Analysis of Immune Cell Infiltrates Across 11 Tumor Types. J Natl Cancer Inst (2016) 108(11):djw144. doi: 10.1093/jnci/djw144

CrossRef Full Text | Google Scholar

8. Nduom EK, Wei J, Yaghi NK, Huang N, Kong LY, Gabrusiewicz K, et al. PD-L1 Expression and Prognostic Impact in Glioblastoma. Neuro Oncol (2016) 18(2):195–205. doi: 10.1093/neuonc/nov172

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Koike N, Yoo SH, Huang HC, Kumar V, Lee C, Kim TK, et al. Transcriptional Architecture and Chromatin Landscape of the Core Circadian Clock in Mammals. Science (2012) 338(6105):349–54. doi: 10.1126/science.1226339

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Ripperger JA, Schibler U. Rhythmic CLOCK-BMAL1 Binding to Multiple E-Box Motifs Drives Circadian Dbp Transcription and Chromatin Transitions. Nat Genet (2006) 38(3):369–74. doi: 10.1038/ng1738

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Kume K, Zylka MJ, Sriram S, Shearman LP, Weaver DR, Jin X, et al. Mcry1 and Mcry2 Are Essential Components of the Negative Limb of the Circadian Clock Feedback Loop. Cell (1999) 98(2):193–205. doi: 10.1016/s0092-8674(00)81014-4

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Sangoram AM, Saez L, Antoch MP, Gekakis N, Staknis D, Whiteley A, et al. Mammalian Circadian Autoregulatory Loop: A Timeless Ortholog and Mper1 Interact and Negatively Regulate CLOCK-BMAL1-Induced Transcription. Neuron (1998) 21(5):1101–13. doi: 10.1016/s0896-6273(00)80627-3

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Cho H, Zhao X, Hatori M, Yu RT, Barish GD, Lam MT, et al. Regulation of Circadian Behaviour and Metabolism by REV-ERB-Alpha and REV-ERB-Beta. Nature (2012) 485(7396):123–7. doi: 10.1038/nature11048

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Ueda HR, Chen W, Adachi A, Wakamatsu H, Hayashi S, Takasugi T, et al. A Transcription Factor Response Element for Gene Expression During Circadian Night. Nature (2002) 418(6897):534–9. doi: 10.1038/nature00906

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Chen Z, Liu P, Li C, Luo Y, Chen I, Liang W, et al. Deregulated Expression of the Clock Genes in Gliomas. Technol Cancer Res Treat (2013) 12(1):91–7. doi: 10.7785/tcrt.2012.500250

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Duhart JM, Brocardo L, Caldart CS, Marpegan L, Golombek DA. Circadian Alterations in a Murine Model of Hypothalamic Glioma. Front Physiol (2017) 8:864. doi: 10.3389/fphys.2017.00864

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Dong Z, Zhang G, Qu M, Gimple RC, Wu Q, Qiu Z, et al. Targeting Glioblastoma Stem Cells Through Disruption of the Circadian Clock. Cancer Discov (2019) 9(11):1556–73. doi: 10.1158/2159-8290.CD-19-0215

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Chen P, Hsu WH, Chang A, Tan Z, Lan Z, Zhou A, et al. Circadian Regulator CLOCK Recruits Immune-Suppressive Microglia Into the GBM Tumor Microenvironment. Cancer Discov (2020) 10(3):371–81. doi: 10.1158/2159-8290.CD-19-0400

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Mure LS, Le HD, Benegiamo G, Chang MW, Rios L, Jillani N, et al. Diurnal Transcriptome Atlas of a Primate Across Major Neural and Peripheral Tissues. Science (2018) 359(6381):eaao0318. doi: 10.1126/science.aao0318

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Wagner PM, Prucca CG, Velazquez FN, Sosa Alderete LG, Caputto BL, Guido ME. Temporal Regulation of Tumor Growth in Nocturnal Mammals: In Vivo Studies and Chemotherapeutical Potential. FASEB J (2021) 35(2):e21231. doi: 10.1096/fj.202001753R

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, et al. TCGAbiolinks: An R/Bioconductor Package for Integrative Analysis of TCGA Data. Nucleic Acids Res (2016) 44(8):e71. doi: 10.1093/nar/gkv1507

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. M(6)A Regulator-Mediated Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Gastric Cancer. Mol Cancer (2020) 19(1):53. doi: 10.1186/s12943-020-01170-0

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdottir H, Tamayo P, Mesirov JP. Molecular Signatures Database (MSigDB) 3. 0. Bioinformatics (2011) 27(12):1739–40. doi: 10.1093/bioinformatics/btr260

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K. KEGG: New Perspectives on Genomes, Pathways, Diseases and Drugs. Nucleic Acids Res (2017). doi: 10.1093/nar/gkw1092

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Jassal B, Matthews L, Viteri G, Gong C, Lorente P, Fabregat A, et al. The Reactome Pathway Knowledgebase. Nucleic Acids Res (2020) 48(D1):D498–503. doi: 10.1093/nar/gkz1031

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Martens M, Ammar A, Riutta A, Waagmeester A, Slenter DN, Hanspers K, et al. WikiPathways: Connecting Communities. Nucleic Acids Res (2021) 49(D1):D613–D21. doi: 10.1093/nar/gkaa1024

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Xuan W, Khan F, James CD, Heimberger AB, Lesniak MS, Chen P. Circadian Regulation of Cancer Cell and Tumor Microenvironment Crosstalk. Trends Cell Biol (2021) 31(11):940–50. doi: 10.1016/j.tcb.2021.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Partch CL, Green CB, Takahashi JS. Molecular Architecture of the Mammalian Circadian Clock. Trends Cell Biol (2014) 24(2):90–9. doi: 10.1016/j.tcb.2013.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Patke A, Young MW, Axelrod S. Molecular Mechanisms and Physiological Importance of Circadian Rhythms. Nat Rev Mol Cell Biol (2020) 21(2):67–84. doi: 10.1038/s41580-019-0179-2

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Wilkerson MD, Hayes DN. ConsensusClusterPlus: A Class Discovery Tool With Confidence Assessments and Item Tracking. Bioinformatics (2010) 26(12):1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Hanzelmann S, Castelo R, Guinney J. GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

32. Yu G, Wang LG, Han Y, He QY. Clusterprofiler: An R Package for Comparing Biological Themes Among Gene Clusters. OMICS (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring Tumour Purity and Stromal and Immune Cell Admixture From Expression Data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Tian R, Li Y, Liu Q, Shu M. Identification and Validation of an Immune-Associated RNA-Binding Proteins Signature to Predict Clinical Outcomes and Therapeutic Responses in Glioma Patients. Cancers (Basel) (2021) 13(7):1730. doi: 10.3390/cancers13071730

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust Enumeration of Cell Subsets From Tissue Expression Profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Aran D, Hu Z, Butte AJ. Xcell: Digitally Portraying the Tissue Cellular Heterogeneity Landscape. Genome Biol (2017) 18(1):220. doi: 10.1186/s13059-017-1349-1

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Zeng D, Ye Z, Shen R, Yu G, Wu J, Xiong Y, et al. IOBR: Multi-Omics Immuno-Oncology Biological Research to Decode Tumor Microenvironment and Signatures. Front Immunol (2021) 12:687975. doi: 10.3389/fimmu.2021.687975

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Xu L, Deng C, Pang B, Zhang X, Liu W, Liao G, et al. TIP: A Web Server for Resolving Tumor Immunophenotype Profiling. Cancer Res (2018) 78(23):6575–80. doi: 10.1158/0008-5472.CAN-18-0689

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Yuan Y, Chen J, Wang J, Xu M, Zhang Y, Sun P, et al. Development and Clinical Validation of a Novel 4-Gene Prognostic Signature Predicting Survival in Colorectal Cancer. Front Oncol (2020) 10:595. doi: 10.3389/fonc.2020.00595

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Chong W, Shang L, Liu J, Fang Z, Du F, Wu H, et al. M(6)A Regulator-Based Methylation Modification Patterns Characterized by Distinct Tumor Microenvironment Immune Profiles in Colon Cancer. Theranostics (2021) 11(5):2201–17. doi: 10.7150/thno.52717

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T Cell Dysfunction and Exclusion Predict Cancer Immunotherapy Response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Miao YR, Zhang Q, Lei Q, Luo M, Xie GY, Wang H, et al. ImmuCellAI: A Unique Method for Comprehensive T-Cell Subsets Abundance Prediction and Its Application in Cancer Immunotherapy. Adv Sci (Weinh) (2020) 7(7):1902880. doi: 10.1002/advs.201902880

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Early JO, Menon D, Wyse CA, Cervantes-Silva MP, Zaslona Z, Carroll RG, et al. Circadian Clock Protein BMAL1 Regulates IL-1beta in Macrophages via NRF2. Proc Natl Acad Sci USA (2018) 115(36):E8460–8. doi: 10.1073/pnas.1800431115

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Lee IK, Song H, Kim H, Kim IS, Tran NL, Kim SH, et al. RORalpha Regulates Cholesterol Metabolism of CD8(+) T Cells for Anticancer Immunity. Cancers (Basel) (2020) 12(7):1733. doi: 10.3390/cancers12071733

CrossRef Full Text | Google Scholar

45. Yu X, Rollins D, Ruhn KA, Stubblefield JJ, Green CB, Kashiwada M, et al. TH17 Cell Differentiation Is Regulated by the Circadian Clock. Science (2013) 342(6159):727–30. doi: 10.1126/science.1243884

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Hu X, Liu X, Moisan J, Wang Y, Lesch CA, Spooner C, et al. Synthetic RORgamma Agonists Regulate Multiple Pathways to Enhance Antitumor Immunity. Oncoimmunology (2016) 5(12):e1254854. doi: 10.1080/2162402X.2016.1254854

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Chabanon RM, Pedrero M, Lefebvre C, Marabelle A, Soria JC, Postel-Vinay S. Mutational Landscape and Sensitivity to Immune Checkpoint Blockers. Clin Cancer Res (2016) 22(17):4309–21. doi: 10.1158/1078-0432.CCR-16-0903

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Chalmers ZR, Connelly CF, Fabrizio D, Gay L, Ali SM, Ennis R, et al. Analysis of 100,000 Human Cancer Genomes Reveals the Landscape of Tumor Mutational Burden. Genome Med (2017) 9(1):34. doi: 10.1186/s13073-017-0424-2

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Balar AV, Galsky MD, Rosenberg JE, Powles T, Petrylak DP, Bellmunt J, et al. Atezolizumab as First-Line Treatment in Cisplatin-Ineligible Patients With Locally Advanced and Metastatic Urothelial Carcinoma: A Single-Arm, Multicentre, Phase 2 Trial. Lancet (2017) 389(10064):67–76. doi: 10.1016/S0140-6736(16)32455-2

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Schumacher TN, Scheper W, Kvistborg P. Cancer Neoantigens. Annu Rev Immunol (2019) 37:173–200. doi: 10.1146/annurev-immunol-042617-053402

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Aslan K, Turco V, Blobner J, Sonner JK, Liuzzi AR, Nunez NG, et al. Heterogeneity of Response to Immune Checkpoint Blockade in Hypermutated Experimental Gliomas. Nat Commun (2020) 11(1):931. doi: 10.1038/s41467-020-14642-0

PubMed Abstract | CrossRef Full Text | Google Scholar