Circadian Clock Genes Are Correlated with Prognosis and Immune Cell Infiltration in Colon Adenocarcinoma

Background Colon adenocarcinoma (COAD) is a malignancy with a high incidence and is associated with poor quality of life. Dysfunction of circadian clock genes and disruption of normal rhythms are associated with the occurrence and progression of many cancer types. However, studies that systematically describe the prognostic value and immune-related functions of circadian clock genes in COAD are lacking. Methods Genomic data obtained from The Cancer Genome Atlas (TCGA) database was analyzed for expression level, mutation status, potential biological functions, and prognostic performance of core circadian clock genes in COAD. Their correlations with immune infiltration and TMB/MSI score were analyzed by Spearman's correlation analysis. Pearson's correlation analysis was performed to analyze their associations with drug sensitivity. Lasso Cox regression analysis was performed to construct a prognosis signature. Moreover, an mRNA-miRNA-lncRNA regulatory axis was also detected by ceRNA network. Results In COAD tissues, the mRNA levels of CLOCK, CRY1, and NR1D1 were increased, while the mRNA levels of ARNTL, CRY2, PER1, PER3, and RORA were decreased. We also summarized the relative genetic mutation variation landscape. GO and KEGG pathway analyses demonstrated that these circadian clock genes were primarily correlated with the regulation of circadian rhythms and glucocorticoid receptor signaling pathways. COAD patients with high CRY2, NR1D1, and PER2 expression had worse prognosis. A prognostic model constructed based on the 9 core circadian clock genes predicted the COAD patients' overall survival with medium to high accuracy. A significant association between prognostic circadian clock genes and immune cell infiltration was found. Moreover, the lncRNA KCNQ1OT1/hsa-miRNA-32-5p/PER2/CRY2 regulatory axis in COAD was also detected through a mRNA-miRNA-lncRNA network. Conclusion Our results identified CRY2, NR1D1, and PER2 as potential prognostic biomarkers for COAD patients and correlated their expression with immune cell infiltration. The lncRNA KCNQ1OT1/hsa-miRNA-32-5p/PER2/CRY2 regulatory axis was detected in COAD and might play a vital role in the occurrence and progression of COAD.


Introduction
Colon adenocarcinoma (COAD) has become a severe health burden that causes over half a million cancerrelated deaths annually [1]. COAD is the most common type of colon tumor. In many developing countries, patients with COAD are usually diagnosed at an advanced stage due to the lack of obvious symptoms at an early stage, related screening projects, and screening instruments [2]. The quality of life of advanced COAD patients is very low due to severe digestive symptoms [3]. Although surgical treatment, targeted treatment, chemotherapy, and immunotherapy have been flexibly utilized in the treatment of COAD, the prognosis of COAD patients is still disappointing [1]. It has been reported that the overall five-year survival rate of COAD patients is less than 40%, mainly due to postoperative recurrence and metastasis [4,5]. Recently, aberrant expression of genes in tumor tissues and their potential applications in clinical practice have been a focus of clinicians. Thus, exploring potential gene signatures for therapeutic and prognostic assessments of COAD is significant clinically.
Physiological circadian rhythms are vital for various biological functions, such as metabolism, regeneration, immunology, and endocrinology [6]. All human tissues are modulated by the circadian clock through extremely complicated processes. The molecular circadian clock consists of a core clock gene loop, which mainly includes CLOCK, ARNTL, PER1, PER2, PER3, CRY1, CRY2, NR1D1, and RORA [6,7]. In recent years, an increasing number of studies have attempted to explore the association between the circadian clock and cancer. Previous studies have found that disruption of normal rhythms and dysfunction of core clock gene loops contribute to the occurrence and development of many cancer types and influence the functions of the tumor immunity cycle [8][9][10][11]. For cancer therapy, previous studies revealed that cancer chronotherapy, i.e., timing chemotherapy administration based on circadian rhythms, might reduce the toxicity of drugs [12]. Some studies based on animal models have reported that drugs targeting biological clock targets, such as RORγ synthetic agonists and RORα synthetic agonists, could promote the activation of anticancer immunity [13,14]. Nevertheless, the role of the circadian clock in prognosis evaluation and its clinical significance in COAD are rarely discussed.
With the rapid development of bioinformatics databases and the enrichment of genomic data, comprehensive analysis of the role of the circadian clock in COAD has become feasible. The present study explores the expression profiles, prognostic value, and related regulatory axis of the circadian clock in COAD, which may offer novel insights into prognostic biomarkers and biological clock targets for COAD.

Materials and Methods
2.1. Datasets. On May 20, 2021, COAD patient genomic data were downloaded from The Cancer Genome Atlas (TCGA) database. TCGA-COAD dataset (N = 455) was selected to perform the analyses. The corresponding clinical data, mainly including sex, age, tumor grade, and survival data, were also obtained from TCGA database on May 20, 2021. Furthermore, the copy number variation (CNV) and somatic mutation data in COAD were also extracted from TCGA database using UCSC Xena (https://xena.ucsc.edu/), an easy-to-use web portal for data files derived from TCGA. Setting the P value cutoff as 0.05, statistical analyses were conducted using R software V4.0.3. Before subsequent processing, the expression data were standardized to transcripts per kilobase million (TPM). To further analyze the expression of circadian clock genes, another gene expression profile of COAD was obtained from Oncomine (https://www .oncomine.org) [15], a cancer microarray database and integrated data-mining platform.

Identification of Differentially Expressed Core Circadian
Clock Genes. A total of 9 genes were selected to perform further analysis, which included CLOCK, ARNTL, PER1, PER2, PER3, CRY1, CRY2, NR1D1, and RORA [6,7]. The profiles of 24-hour oscillations of the core circadian clock genes were downloaded from CircaDB (https://circadb.hogeneschlab .org) [16]. For each tissue that was cyclic ordering by periodic structure (CYCLOPS) ordered, cosinor regression was used to test if individual genes are rhythmic or not [17]. In the 24-hour oscillation plots, the x-axis represented the CYCLOPS-estimated sample phase in radians (from 0 to 2?). The y-axis represented the expression level (TPM). The time ? represented E-box phase (i.e., time of peak expression of E-box target genes NR1D1, NR1D2, and PER3). The datasets concerning the expression levels of the 9 core circadian clock genes were also excavated from the Oncomine (https://www.oncomine.org) database [15]. To further verify the expression level of the 9 core circadian clock genes, we also analyzed their expression in TCGA-COAD dataset. Student's t-test was utilized to compare the transcriptional levels of the core circadian clock genes in COAD tissues with those in normal colon tissues using the "reshape2" and "limma" packages in R software.

PPI Network and Mutation Analysis of Core Circadian
Clock Genes. To study the correlation between each member of the 9 core circadian clock genes, we constructed a PPI network with the Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org/), a database of predicted functional associations between genes and proteins [18]. Using R software, the overall mutation status of core circadian clock genes in COAD, including mutation categories and mutation frequencies of CLOCK, ARNTL, PER1, PER2, PER3, CRY1, CRY2, NR1D1, and RORA as well as their waterfall plot, were constructed using the "maftools" package. The locations of variants of the 9 genes on 23 chromosomes were shown using the "RCircos" package.

Functional Enrichment Analysis for Core Circadian
Clock Genes. To identify the potential molecular mechanisms and biological functions of core circadian clock genes in COAD, functional enrichment analysis, including gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, was performed using the "ggplot2" package in R software. In addition, GO analysis included biological process (BP), cellular component (CC), and molecular function (MF) analyses.

Construction of a Prognostic Model for Core Circadian
Clock Genes. To assess the potential prognostic value of core circadian clock genes, Kaplan-Meier plots, including overall survival (OS) plots, disease-free survival (DFS) plots, disease-specific survival (DSS) plots, and progression-free survival (PFS) plots, were drawn. The hazard ratio (HR) and its 95% confidence interval (CI) are also listed in the plots. Based on the 9 core circadian clock genes, a LASSO Cox regression model was used to predict the influence of the expression of core circadian clock genes on the prognosis of COAD patients. According to the median risk score, patients were divided into the low-and high-risk groups, and the OS curves of the two groups were compared. The predictive accuracy of this prognostic model was assessed by time ROC analysis. Moreover, core circadian clock genes 2 Computational and Mathematical Methods in Medicine with potential prognostic value were selected for further analyses.
2.6. Immune Infiltration, TMB, MSI, and Drug Sensitivity Analyses. TIMER (https://cistrome.shinyapps.io/timer/) is a bioinformatics platform that is aimed at visualizing the association between immune cell infiltration and target gene expression in cancer cells [19]. The "Gene" module in TIMER was utilized to determine the correlation between the expression of prognostic circadian clock genes and the immune infiltration levels of B cells, CD8+ T cells, CD4+ T cells, neutrophils, dendritic cells, and macrophages. Based on the "ggstatsplot" package in R, tumor mutation burden (TMB) and microsatellite instability (MSI) analyses were used to assess the relationship between the expression of prognostic circadian clock genes and TMB and MSI scores, respectively. Spearman's correlation test was selected for the above analysis with a P value cutoff of 0.05. Moreover, to study the effect of prognostic circadian clock gene expression on drug resistance, drug sensitivity data in the Genomics of Drug Sensitivity in Cancer (GDSC) database were extracted.
An mRNA-miRNA-lncRNA network was constructed to determine the potential regulatory axis of prognostic circadian clock genes in COAD. The miRNA targets binding to prognostic circadian clock genes were identified by the miR-TarBase (http://mirtarbase.mbc.nctu.edu.tw/) [20] and star-Base (https://starbase.sysu.edu.cn/index.php) [21] databases. miRTarBase and starBase are bioinformatics tools that can be used to predict miRNA targets of mRNAs. The mRNA-miRNA regulatory network was constructed using Cytoscape software. The Cytohubba plug-in was used to select the most significantly connected miRNAs. Next, according to the selected miRNAs, the lncRNA targets interacting with miRNAs were extracted from starBase (https:// starbase.sysu.edu.cn/index.php) and lncBase (https://diana .e-ce.uth.gr/lncbasev3) [22], and the associations between the miRNAs and the lncRNAs were visualized. starBase and lncBase are bioinformatics tools that can be used to predict lncRNA targets of miRNAs. To further confirm the significant role of miRNAs and lncRNAs, we also evaluated the expression levels and prognostic performances of the targeted miRNAs and lncRNAs using TCGA-COAD dataset.

Work Flow of the Current Study.
Work flow of the present study is shown in Figure 1. Firstly, the core circadian clock genes in human were identified based on the previous studies [6,7]. To reveal their rhythmicity in human colon tissues, the 24-hour oscillations of these genes were listed based on the expression profiles in CircaDB. Next, to find out the differentially expressed genes, we downloaded the RNA-sequencing data in Oncomine and TCGA-COAD dataset, respectively. The datasets in Oncomine were from the previous publications. Through comparing the expression profiles in the two databases, the expression status of core circadian clock genes in COAD was more trustworthy. And then, based on clinical data in TCGA-COAD database, the potential prognostic value of these genes in COAD was also analyzed by Kaplan-Meier methods, which included OS, DFS, PFS, and DSS. After screening out the prognostic circadian clock genes in COAD, we were also interested on the correlation between prognostic circadian clock genes expression and immune cell infiltration; these results might assist researchers to explore potential immunotherapeutic target for COAD patients. Immune infiltration analysis was conducted by TIMER. MSI and TMB could predict the efficacy of cancer immunotherapy. In theory, cancer patients with high TMB/MSI score might get better therapeutic effect. Thus, we detected the influence of prognostic circadian clock genes expression on TMB/MSI score. Meanwhile, the association between gene expression and drug sensitivity of colon tissues was also analyzed. To explore potential genes that might influence the development of COAD, we performed the subgroup analysis to show the correlation between gene expression and clinical stage. In addition, we attempted to explore a potential mRNA-miRNA-lncRNA 3 Computational and Mathematical Methods in Medicine regulatory axis concerning these prognostic circadian clock genes in COAD. To comprehensively investigate these 9 genes in COAD, the following analyses were performed. Mutation analysis helped readers to intuitively observe the overall mutation status of circadian clock genes in COAD. As the core genes in circadian clock loop, the potential biological functions and molecular mechanisms of these 9 genes were investigated by functional enrichment analysis. These results would provide novel insights on the functional studies of circadian clock genes. And in recent years, machine learning and big data technique were used to predict prognosis. Thus, LASSO Cox regression model based on the core circadian clock genes was constructed. Finally, the PPI net-work was constructed to visualize the correlation between each member of the 9 core circadian clock genes. Correlation analysis was performed to show the coexpression among these 9 genes in COAD, thus revealing synergistic interaction on expression among these genes.   Computational and Mathematical Methods in Medicine in CircaDB. Based on the Oncomine database, the number of datasets that contained the expression profiles of CLOCK, ARNTL, PER1, PER2, PER3, CRY1, CRY2, NR1D1, and RORA is shown in Figure 3. In colorectal cancer, the transcriptional levels of CLOCK, CRY1, and NR1D1 were upregulated in 1, 1, and 2 datasets, respectively. The transcriptional levels of ARNTL, CRY2, PER1, PER3, and RORA were downregulated in 2, 6, 2, 5, and 5 datasets, respectively. Using TCGA-COAD dataset, the mRNA levels of the 9 core circadian clock genes in TCGA-COAD tissues and the corresponding mRNA levels in normal colon tissues were identified and are shown in Figure 4(a). Specifically, the mRNA levels of CLOCK, CRY1, and NR1D1 were significantly higher in COAD tissues, while the mRNA levels of ARNTL, CRY2, PER1, PER3, and RORA were significantly decreased in COAD tissues (all P values < 0.001). At the transcriptional level, these core circadian clock genes were greatly altered in COAD tissues. Moreover, a protein-protein interaction (PPI) network was constructed to illustrate the associations among the 9 core circadian clock genes (Supplementary Figure 1(a)). Correlation analysis of core circadian clock genes showed significant positive correlations with each other apart from NR1D1 in TCGA-COAD dataset, suggesting that these core circadian clock genes had synergistic interactions on expression (Supplementary Figure 1(b)).

Landscape of the Genetic Variation of Core Circadian
Clock Genes in COAD Tissues. Next, we summarized the somatic variations (SNVs) and copy number variations (CNVs) of the 9 circadian clock genes in COAD. As shown in Figures 4(b) and 4(c), genetic alterations were detected in 59 (88.06%) of 67 COAD samples. The categories of genetic mutations included missense mutation, nonsense mutation, frame-shift insertion, frame-shift deletion, in-frame deletion, splice-site variation, and multi-hit mutation. The most common variant classification was missense mutation. In addition, SNPs were the most common variant type, and C>T ranked as the top SNV class. Of these genes, PER3 was the gene with the highest mutation frequency. The locations of the CNV variations on chromosomes are also presented in Figure 4(d). The CNV variation frequencies of the 9 core circadian clock genes in COAD are also summarized in Figure 4(e). Copy number amplification was shown in NR1D1 and CLOCK, while copy number deletion was shown in PER3, PER2, RORA, ARNTL, CRY1, and PER1.

Functional Enrichment Analysis of Core Circadian Clock
Genes in COAD Tissues. To clarify the potential biological functions and molecular mechanisms of core circadian clock genes in COAD, GO and KEGG pathway analyses were performed ( Figure 5). GO analysis, including BP, CC, and MF analyses, illustrated that these 9 genes mainly participated in the regulation of circadian rhythm, circadian regulation of gene expression, rhythmic process, glucocorticoid receptor signaling pathway, corticosteroid receptor signaling pathway, transcription corepressor binding, and E-box binding ( Figure 5(a)). KEGG pathway analysis revealed that the 9 core circadian clock genes were mainly involved in the regulation of circadian rhythm, circadian entrainment and dopaminergic synapses ( Figure 5(b)).

Prognostic Model Construction of Core Circadian Clock
Genes in COAD Tissues. Kaplan-Meier plots, including OS, DFS, PFS, and DSS plots, were used to assess the potential prognostic value of the 9 core circadian clock genes. The P         Table 1. COAD patients with high CRY2 expression showed worse OS (Figure 6 54)). In addition, COAD patients with high PER2 expression showed worse OS ( Figure 6(b)) (P value = 0.047, HR with 95%CI = 1:49 (1.00-2.21)). The above results demonstrated that CRY2, NR1D1 and PER2 might be of great significance in the prognostic evaluation of COAD patients. Moreover, a prognostic model including the 9 core circadian clock genes was constructed by Lasso Cox regression analysis (Figure 7(a)-7(b)): risk score = ð−0:1332Þ * CLOCK + ð0:4144Þ * CRY2 + ð0:0934Þ * NR1D1: According to the risk score, patients were divided into the high-and low-risk groups. The score distribution, survival status, and expression levels of PER2, NR1D1, and CLOCK are shown in Figure 7(c). When the risk score increased, the risk of death of COAD patients increased, while the survival time decreased (Figure 7(c)). Compared with the low-risk group, COAD patients with high risk showed a worse OS (median survival time = 5:4 vs. 8.2 years, P value = 0.015, Figure 7(d)). Furthermore, the AUCs of the 5-and 10-year ROC curves were 0.604 (95% CI: 0.512-0.697) and 0.696 (95% CI: 0.536-0.856), respectively (Figure 7(e)), which showed medium to high accuracy.
3.6. Immune Infiltration Analysis for Prognostic Circadian Clock Genes in COAD Tissues. The above study revealed that CRY2, NR1D1, and PER2 were potential biomarkers for COAD. We then selected CRY2, NR1D1, and PER2 for further analysis because of their prognostic performance. The association between the infiltration level of immune cells and the expression level of prognostic circadian clock genes in COAD was analyzed. As shown in Figure 8 (Figure 8(c)). We were also interested in the correlation between the common immune biomarkers and the expression levels of prognostic circadian clock genes in COAD. As shown in Table 2, CRY2 expression was positively correlated with most of the gene markers on immune cells. However, PER2 and NR1D1 expression levels were significantly correlated with only 24 and 29 kinds of gene markers, respectively. These results revealed that CRY2 might be a potential immunotherapeutic target for COAD patients.       The above results demonstrated that the prognostic circadian clock genes were associated with immune cell infiltration. Therefore, we conducted TMB and MSI analyses to assess the potential clinical applications of CRY2, NR1D1, and PER2 in COAD immunotherapy. In TMB analysis, negative correlations were detected between TMB and CRY2 expression (Figure 9(a), P value = 0.002) and NR1D1 expression (Figure 9(c), P value = 0.048), while no significant correlation was found between TMB and PER2 expression (Figure 9(b), P value = 0.210). In MSI analysis, unexpectedly, no significant associations between MSI and CRY2 expression (Figure 9(d), P value = 0.769),   (Figure 9(e), P value = 0.948), or PER2 expression (Figure 9(f), P value = 0.588) were detected. Next, drug sensitivity analysis was performed to detect the potential functions of prognostic circadian clock genes in drug screening. Based on data from the GDSC database, high NR1D1 expression in COAD tissues was mainly correlated with drug resistance to 38 kinds of drugs and drug sensitivity to 6 kinds of drugs. High PER2 expression was primarily associated with drug sensitivity in 58 kinds of drugs and drug resistance in 13 kinds of drugs, while high CRY2 expression was mainly associated with drug sensitivity in 60 kinds of drugs and drug resistance in 18 kinds of drugs (Figure 9(g)), indicating that PER2 and CRY2 were potential biomarkers for drug screening in COAD.  Figure 10(c), P value = 0.0089) were associated with the clinical stage of COAD patients, which demonstrated that these 3 genes might be involved in the development of COAD. However, for potential applications in drug screening and prognostic performance, PER2 and CRY2 were selected to construct an mRNA-miRNA-lncRNA network. The miRNAs targeting PER2 and CRY2 identified by starBase and miRTarBase are shown in Figure 11(a). Among these miRNAs, hsa-miRNA-32-5p, hsa-miRNA-340-5p, hsa-miRNA-20b-5p, and hsa-let-7b-5p were identified as the top highly connected miRNAs (Figure 11(b)), which were upregulated or downregulated in COAD tissues. We then detected the expression profiles and prognostic value of these four genes. The results showed that only hsa-miRNA-32-5p was differentially expressed in tumors and was significantly associated with patient prognosis. More specifically, the expression of hsa-miRNA-32-5p was upregulated in COAD tissues compared with normal colon tissues (P value = 1:3E − 06, Figure 11(c)). COAD patients with high hsa-miRNA-32-5p levels showed worse OS (P value = 0.03453, Figure 11(d)). Thus, hsa-miRNA-32-5p was selected as the most promising miRNA target of PER2 and CRY2. Next, we further detected the upstream lncRNA targets using lncBase and starBase. According to the data from the above databases, a total of 11 lncRNA targets were identified (Figure 11(e)), including NDUFA6-AS1, SNHG14, CTBP1-AS2, LINC01128, INTS6-AS1, XIST, LINC00365, KCNQ1OT1, OIP5-AS1, DAPK1-IT1, and PITPNA-AS1. The expression levels of the above lncRNA targets in COAD tissues and normal colon tissues were analyzed. Similarly, the results showed that only lncRNA KCNQ1OT1 was differentially expressed in tumors and was significantly associated with patient prognosis. The data demonstrated upregulation of the lncRNA KCNQ1OT1 in COAD tissues compared with normal colon tissues (P value = 3:2E − 16, Figure 11(f)). COAD patients with high KCNQ1OT1 expression had a worse OS than those with low KCNQ1OT1 expression (median survival time = 5.3 years vs. 7.9 years, P value = 0.039, Figure 11(g)). Thus, the lncRNA KCNQ1OT1/hsa-miRNA-32-5p/PER2/CRY2 regulatory axis might influence the occurrence and progression of COAD.

Discussion
In 1971, scientists first found three variations in one gene on the X chromosome of Drosophila melanogaster; one caused physiological response and behavior to become completely asynchronous with the circadian clock, and the other two caused the circadian clock to change to 19 and 28 hours, respectively [23]. By analyzing the entire X chromosome, researchers identified the location of the gene and named the gene PER. After this groundbreaking study, many studies attempted to expand the circadian clock gene family and determine the potential roles of circadian clock genes in the occurrence and development of disease and their potential applications in disease diagnosis, treatment and prevention. In the past 20 years, the circadian clock has been     [6,7]. Previous studies reported that dysfunction and abnormal expression of these genes were involved in the occurrence and progression of various human cancer types [8,9,24,25]. Some studies have focused on the potential of core circadian clock genes to be cancer biomarkers [26][27][28]. Nevertheless, the specific functions of core circadian clock genes in the prognosis and therapy of COAD are rarely discussed. Thus, we adopted bioinformatics methods to illustrate the potential role of core circadian clock genes in COAD.
To a certain extent, the expression of core circadian clock genes reflects the status of the circadian clock and influences its functions [29]. In mammals, interlocked transcriptional-translational feedback loops regulate the circadian clock. In normal tissues, the positive limb of this loop is composed of CLOCK and ARNTL, which induce the expression of CRY1-2 and PER1-3. In turn, overexpressed PER and CRY proteins suppress their own transcription, thus forming the circadian rhythm [30,31]. Compared to those in normal colon tissues, the transcriptional levels of CLOCK, CRY1, and NR1D1 in COAD tissues were significantly higher, while the mRNA levels of ARNTL, CRY2, PER1, PER3, and RORA in COAD tissues were downregulated in the current study. The dysexpression of these genes in COAD tissues might be involved in the disruption of cir-cadian rhythms, thus leading to the occurrence and progression of COAD.
Researchers are also interested in the potential prognostic value of core circadian clock genes in COAD. In the present study, COAD patients with high CRY2 levels showed worse OS, DFS, DSS, and PFS than those with low CRY2 levels. COAD patients with high NR1D1 expression had worse DFS and PFS. In addition, high PER2 expression predicted a worse OS for COAD patients. A previous study based on human colorectal cancer tissues indicated that CRY2 expression was decreased and that high CRY2 expression predicted worse OS, which supported the results of the current study [26]. Na et al. indicated that NR1D1 was considered a biomarker predicting good prognosis in breast cancer patients [27]. PER2 is considered a prognostic predictor in several human cancer types, including lung cancer, colorectal cancer, and gastric cancer [28,32,33]. Moreover, the overall prognostic performances of CLOCK, PER1, PER3, and CRY1 were also determined in a previous metaanalysis [34].
Tumor-infiltrating immune cells play pivotal roles in eliminating cancer cells and hindering cancer progression. Exploring the association between core circadian clock genes and tumor-infiltrating immune cells is of significant importance. A previous study concluded the role of the circadian clock and its core genes in cancer immunity [35]. Xu et al. reported an M1-like proinflammatory phenotype of macrophages in mice with relatively disrupted PER1 and PER2     [37]. However, more high-quality studies are needed to illustrate the mechanisms and specific processes of core circadian clock genes in cancer immunity. In addition, we detected the associations between the expression of core circadian clock genes in COAD and biomarkers in immune cells, which might assist researchers in determining potential immunotherapy targets for COAD. RORγ synthetic agonists could increase the expression levels of IL-17A, IL-17F, GM-CSF, and IL-22, thus boosting Th17 cells and Tregs to block immunosuppression [38]. The expression of BMAL1, CLOCK, Rev-Erbα, and PER2 in intraperitoneal macrophages regulated the expression levels of F4/80 and CD11c in tumor tissues [39]. However, immunotherapy based on core circadian clock gene targets is still at the initial stage. Thus, more complementary studies aimed at core circadian clock genes are urgently needed. GO and KEGG pathway analyses could predict the biological functions and molecular mechanisms of the target gene set. In our study, the main functions and mechanisms of the 9 core circadian clock genes were involved in the regulation of circadian rhythms and some signaling pathways. The corticosteroid receptor signaling pathway includes the glucocorticoid receptor signaling pathway and mineralocorticoid receptor signaling pathway. In a mouse model, So et al. identified glucocorticoid response elements in the molecular circadian clock and found that PER2 was an integral component of the glucocorticoid regulatory pathway and was directly regulated by the glucocorticoid receptor [40]. Dickmeis et al. described multiple mechanisms by which molecular circadian clock genes regulate rhythms of glucocorticoid release and modulate glucocorticoid signaling [41]. For the mineralocorticoid receptor signaling pathway, circadian gene CLOCK signaling contributed to mineralocorticoid receptor-mediated cardiac inflammation and fibrosis [42]. The above findings revealed the significance of circadian clock genes in the corticosteroid receptor signaling pathway.
Through mRNA-miRNA-lncRNA network construction, the lncRNA KCNQ1OT1/hsa-miRNA-32-5p/PER2/ CRY2 regulatory axis was identified. Previous studies have reported that miRNA-32-5p regulates the progression of gastric cancer and epithelial-mesenchymal transition and metastasis in lung adenocarcinoma [43,44]. However, studies concerning the role of miRNA-32-5p in COAD are lacking. The upstream lncRNA KCNQ1OT1 was also detected and is considered a potential diagnostic biomarker for colon and rectal cancer and enhanced drug resistance in colon cancer cells [45,46]. In the present study, we also studied the prognostic performance of miRNA-32-5p and lncRNA KCNQ1OT1 and found that they were correlated with COAD patient prognosis. The above findings revealed that the lncRNA KCNQ1OT1/hsa-miRNA-32-5p/PER2/CRY2 regulatory axis might be of significant importance in the occurrence and progression of COAD and may have potential applications in the diagnosis and treatment of COAD. Nevertheless, in vivo and vitro studies are urgently needed to clarify the role of this regulatory axis in COAD.
Some limitations still exist in the current study. First, our study and analysis were mainly performed at the transcriptional level, and the results may not be applicable in studies based on the protein level. Second, related fundamental and clinical studies that focus on the molecular mechanisms of circadian clock genes in colon or rectal cancer are rare. Third, the genetic background and etiology of COAD patients are influenced by many factors, such as patient race, sex, and age. Thus, more in-depth studies are necessary to validate the role of circadian clock genes in COAD.

Conclusion
In this study, a comprehensive bioinformatics analysis was used to identify the prognostic circadian clock gene signature in COAD, which included CRY2, NR1D1, and PER2. The lncRNA KCNQ1OT1/hsa-miRNA-32-5p/PER2/CRY2 regulatory axis was identified and might be of great significance in the occurrence and progression of COAD.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.