A circadian clock gene-related signature for predicting prognosis and its association with sorafenib response in hepatocellular carcinoma

Introduction

Hepatocellular carcinoma (HCC), which is one of the leading causes of cancer-related death globally, is characterized by inflammation and dysfunctional metabolism (1,2). Although there are several therapies for HCC, such as surgery, ablation, arterial chemoembolization, and therapy, few patients benefit from these regimens, especially patients at an advanced stage (3,4). In order to improve the efficacy of these treatments for HCC, some clinical parameters (e.g., alpha fetoprotein, tumor size, stage) are commonly used for determining prognosis. However, these parameters fail to provide satisfactory prediction because of the low specificity and accuracy. Therefore, exploring more prognostic markers for facilitating and improving the overall survival (OS) of patients with HCC is needed.

The circadian rhythm is driven by endogenous circadian clock genes and exists in almost all organisms and cells and participates in a variety of physiological functions (5). Mechanistically, the circadian rhythm is a transcription and translation feedback loop (TTFL) that mainly includes 15 core circadian clock genes, CLOCK, ARNTL, CRY1, CRY2, PER1, PER2, PER3, RORA, RORB, RORC, NPAS2, TIMELESS, NR1D1, NR1D2, and CSNK1E (6,7). TTFL governs 43% of all genes in at least one tissue and oscillates over 24 hours (8). Almost all people experience circadian rhythm disruption mainly due to situations of modern life, such as night shift work (9). Accumulating epidemiological and genetic studies have revealed that chronic disruption of circadian rhythms is associated with a few diseases involving metabolic disorders and tumor (9-11). Furthermore, previous research has confirmed the link between the circadian clock and the initiation and progression of HCC (12,13). Cancer chronotherapy is a potentially effective treatment and involves antitumor agents being administered at special times based on the circadian rhythms of therapeutic efficacy and those of the side effects of normal cells (14). David and colleagues revealed that immune checkpoint inhibitors administered in the evening are less effective than those administered during the daytime (15). Therefore, the role of the circadian clock in the disease process of HCC warrants further investigation.

In this study, we systematically evaluated the expression pattern, copy number variation (CNV), and mutation of core circadian clock genes in HCC and identified prognostic circadian clock genes. Moreover, we constructed a messenger RNA signature (mRNASig) for predicting the OS of patients with HCC and divided patients with HCC into high- and low-risk groups. We also analyzed functional enrichment, the immune landscape, the immune checkpoints, and the tumor mutational burden (TMB) between the two groups. The high-risk group was associated with cancer-related pathways, metabolism, and immunosuppressive features. Additionally, we noticed that patients in the low-risk group had more enriched angiogenesis and were more sensitive to sorafenib. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-217/rc).

Methods

Data acquisition and processing

The workflow of this study is depicted in Figure S1. Datasets with transcript files and CNV and mutation data were derived from The Cancer Genome Atlas (TCGA) database, International Cancer Genome Consortium (ICGC; https://dcc.icgc.org), and a Gene Expression Omnibus (GEO) dataset (GSE109211). The corresponding clinical data were also obtained. With TCGA used as a training cohort, we eliminated cases with <1 month of survival or those without survival data because they might die of other complications (e.g., hemorrhage and heart failure) other than HCC. The transcriptome data were standardized with log₂ [transcripts per million (TPM)+1] transformation for subsequent analysis (16). The ICGC dataset as a validation set was processed in the same manner. The transcription profile and clinical data of patients with HCC treated with sorafenib were obtained from a GEO dataset (GSE109211).

Prognostic circadian clock gene-related signature construction and validation

For predicting the survival probability of patients with HCC, we identified the mRNASig based core circadian clock genes by performing least absolute shrinkage and selection operator (LASSO) analysis with the “glmnet” package in R (The R Foundation of Statistical Computing) and multivariate Cox regression analysis (17). The mRNASig was built by using the regression coefficients derived from multivariable Cox regression analysis and the expression of the optimal prognostic biomarkers. The risk score formula was established as follows:

Riskscore=∑iCoefficientof(i)×Expressionofgene(i)[1]

where the coefficient of gene (i) is the regression coefficient of gene (i), and the expression of gene (i) is the expression value of gene (i) for each patient. We divided patients with HCC in TCGA cohort into high- and low-risk groups using the median risk score. To acquire the same formula, the expression of selected circadian clock genes was normalized using standard deviation (SD) =1 and mean value =0 in TCGA, ICGC, and GEO cohorts. Kaplan-Meier (K-M) analysis and the logarithmic rank test were used to evaluate the survival time of the two groups. Next, the predictive accuracy of the signature was confirmed with the time-dependent receiver operating characteristic (TDROC) curve, the area under the curve (AUC), and the concordance index (C-index).

Identification of independent prognostic factors

Univariate and multivariate Cox regression analyses were performed to identify the independent prognostic value of mRNASig and clinical variables (including age, gender, race, grade, stage, and albumin). To evaluate the predictive performance of the mRNASig in different clinical subgroups, K-M analysis was used to estimate the survival difference between different risk assessment groups of the different subgroups.

Construction and validation of a mRNASig-based nomogram

To predict the 1-, 3-, and 5-year OS of patients with HCC in TCGA and ICGC databases, a nomogram based on the risk score and clinicopathological parameters was developed using a stepwise Cox regression model.

Differentially expressed gene (DEG) analysis and functional enrichment analysis

The DEGs between the high- and low-risk groups were selected using the “limma” package in R software. The genes with an adjusted P value <0.05 and |log₂ fold change (FC)| >1 were considered robust DEGs. To investigate the potential functions of mRNASig mRNA, we performed multiple methods, including Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, gene set enrichment analysis (GSEA), and gene set variation analysis (GSVA). Moreover, with the KEGG and Hallmark databases being used as the reference gene set, we also identified a set of gene sets involved in specific biological pathways as described by Mariathasan et al. (18), including (I) angiogenesis, (II) DNA replication, (III) DNA damage repair 1 and 2, (IV) antigen-processing machinery, (V) cell cycle, (VI) Fanconi anemia, (VII) mismatch repair, and (VIII) nucleotide excision repair (the supplementary table is available at https://cdn.amegroups.cn/static/public/tcr-23-217-1.xlsx).

Tumor immunity analysis of the signature mRNAs

The single-sample GSEA (ssGSEA) algorithm was used to quantify the relative abundance of cell infiltration in patients with HCC. Charoentong et al. constructed the gene set for marking each infiltration immune cell type, including regulatory T cells, activated CD4 T cells, macrophages, nature killer (NK) cells, and T helper type 2 (Th2) cells (the supplementary table is available at https://cdn.amegroups.cn/static/public/tcr-23-217-1.xlsx) (19,20).

Statistical analysis

R software (version 4.1.2; https://ww.r-project.org) was used for all computational and statistical analyses. The chi-squared test or Fisher exact test was used to determine the statistical significance of differences between two groups for categorical variables. The relationship between two groups for continuous variables was evaluated with the Mann-Whitney test. P values <0.05 were considered statistically significant.

Ethical statement

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Results

The landscape of core circadian clock genes in HCC

We investigated core circadian clock genes at the transcription level in normal and tumor samples from patients with HCC (Figure 1A). The expression of these genes was dysregulated. Among the 15 circadian clock genes, only RORA and PER1 were significantly downregulated in tumor tissues in comparison with normal tissues. The expressions of PER2, NR1D1, ARNTL, and PER3 in tumor samples were not significantly different from those in normal tissues, while the others were significantly upregulated in tumor tissues. To further ascertain whether the transcription pattern of the circadian clock genes is influenced by genetic variations, we summarized the CNVs and the incidence of somatic mutations in HCC. The locations of CNV alterations in the circadian clock genes on chromosomes are shown in Figure 1B. In tumor samples, we noticed a prevalent CNV alteration in core circadian clock genes. Most genes with higher expression were also amplified in copy number (e.g., RORC, CSNK1E), while PER1, with lower expression, had a high frequency of CNV deletion (Figure 1C). Moreover, we found that patients with CNV amplification of some core circadian clock genes had higher gene expression compared to those with CNV deletion and those with CNV without alteration (Figure S2), which suggested that alterations of CNV could perturb the transcription pattern of the core circadian clock genes. Furthermore, we investigated somatic mutation of circadian clock genes in HCC (Figure 1D). Among 329 samples, only 27 samples showed mutations of circadian clock genes, with a total frequency of 8.21%. Moreover, PER1, PER2, and TIMELESS demonstrated the highest mutation frequency, while ARNATL and RORC did not show any mutations in HCC samples, which indicated that the somatic mutation of core circadian clock genes was rare. These findings suggest that patients with HCC have disrupted circadian clock genes.

Figure 1 Landscape of expression, CNVs, and mutations of core circadian clock genes in HCC. (A) In TCGA, the expression of core circadian clock genes of normal and tumor samples. (B) The location of CNV alteration of core circadian clock genes on chromosomes in TCGA cohort. (C) The CNV variation frequency of core circadian clock genes in TCGA cohort. The height of the column represents the alteration frequency. (D) The mutation frequency of core circadian clock genes in patients with HCC in TCGA cohort. Each column represents an individual patient. The upper bar plot indicates TMB. The number on the right represents the mutation frequency in each core circadian clock gene. The fraction of conversions in each sample is shown in the stacked bar plot below. Colors indicate the types of base substitution. (E,F) LASSO Cox regression analysis for the investigation of core circadian clock genes related to HCC prognosis. *, P<0.05; **, P<0.01, ***, P<0.001; ****, P<0.0001; ns, no significance. CNV, copy number variation; HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; TMB, tumor mutation burden; LASSO, least absolute shrinkage and selection operator.

Development of mRNASig and the assessment of prediction efficacy

According to LASSO and stepwise multivariate Cox regression analysis, NPAS2, NR1D1, PER1, RORC, and TIMELESS were selected from the 15 circadian clock genes for establishing a prognostic model (Figure 1E,1F). The formula for the risk score was as follows: risk score = [NPAS2 expression × 0.2387950] + [NR1D1 expression × 0.2431583] + [PER1 expression × (−0.3110501)] + [RORC expression × (−0.1861446)] + [TIMELESS expression × (0.1558095)]. The patients in the training cohort were stratified into high- and low-risk groups depending on the median value of the risk score. The K-M analysis revealed that the high-risk group showed remarkably reduced OS (P<0.0001) relative to the low-risk group (Figure 2A). The risk score and survival status of patients with HCC and the expression levels of five selected genes in TCGA cohort are shown in Figure 2B. The AUCs of the ROC curve in training cohort were 0.784, 0.691, and 0.651 at 1, 3, and 5 years, respectively, suggesting the mRNASig had excellent predictive performance (Figure 2C). Furthermore, we found that the AUC of mRNASig at 1 year was higher than that of the clinicopathological parameters, suggesting the excellent performance of the risk score signature (Figure 2D). Moreover, the prognostic effects of mRNASig on OS was investigated using a prespecified subgroup analysis. High-risk group showed a poorer prognosis compared to the low-risk group within the subgroups categorized according to age (≤60/>60 years; Figure 3A,3B), sex (Figure 3C,3D), grade (grades 1–2/grades 3–4; Figure 3E,3F), stage (stages I–II /stages III–IV; Figure 3G,3H), and Child-Pugh score (Child-Pugh score A and B; Figure 3I).

Figure 2 Survival analysis and prognostic performance of mRNASig in TCGA cohort. (A) Overall survival curve of mRNASig. (B) Risk score and survival time of patients and the transcript profile of selected circadian clock genes in mRNASig. (C) TDROC curve of mRNASig for predicting HCC. (D) ROC curves of clinical characteristics and mRNASig for predicting HCC. mRNASig, messenger RNA signature; TCGA, The Cancer Genome Atlas; TDROC, time-dependent receiver operating characteristic; HCC, hepatocellular carcinoma; ROC, receiver operating characteristic; AUC, area under the curve.

Figure 3 The mRNASig is available for most of the different subgroups including age (A,B), gender (C,D), grade (E,F), stage (G,H), and Child-Pugh score (I). mRNASig, messenger RNA signature.

Validation of mRNASig

mRNASig was evaluated with the ICGC cohort being used for validation. The K-M analysis of the ICGC cohort revealed that patients from the high-risk group had poorer prognosis than those from the low-risk group (Figure 4A,4B), which was consistent with the results from the training cohort. Analysis of the validation cohort also revealed that mRNASig performed excellently as a prognostic predictor in HCC, with AUCs of 0.774, 0.705, and 0.513 at 1, 3, and 5 years, respectively (Figure 4C). The AUC value of the mRNASig at 1 year was also higher than that of the clinicopathological parameters in the ICGC cohort (Figure 4D). The above results indicated that the mRNASig was a good predictor of the OS of patients with HCC.

Figure 4 Validation of mRNASig in the validation cohort. The ICGC dataset was used as the validation cohort. (A) Kaplan-Meier analysis of mRNASig. (B) Risk score and survival time of patients and gene expression of the selected circadian clock genes in mRNASig. (C) TDROC curve of mRNASig. (D) ROC curves of clinical characters and mRNASig for predicting HCC. mRNASig, messenger RNA signature; ICGC, International Cancer Genome Consortium; ROC, receiver operating characteristic; TDROC, time-dependent ROC; AUC, area under the curve; HCC, hepatocellular carcinoma.

Identification of independent prognostic variables in HCC and construction of the nomogram

To investigate whether the mRNASig was an independent prognostic factor in the OS of patients with HCC, we conducted univariate and multivariate Cox analyses in the training cohort and validation cohort. The mRNASig and stage served as independent prognostic predictors in both cohorts (Figure 5A and Figure S3A). A nomogram based on two independent prognostic factors (including risk score and stage) was constructed for predicting the OS of patients with HCC at 1, 3, and 5 years (Figure 5B and Figure S3B). The C-index of the nomogram was 0.713 (95% confidence interval: 0.685–0.741), and the calibration curves plots demonstrated a high consistency between the actual and predicted OS of the nomogram in patients with HCC (Figure 5C and Figure S3C).

Figure 5 Construction and evaluation of a nomogram for predicting the OS of HCC in the training cohort. (A) Univariate and multivariate Cox regression analyses identified the prognostic factors. (B) The nomogram of a prognostic model for HCC. (C) The calibration curve of the nomogram. ***, P<0.001. OS, overall survival; HCC, hepatocellular carcinoma.

Function enrichment analysis of mRNASig

To identify the molecular mechanisms of mRNASig, we performed a comprehensive functional enrichment analysis of DEGs between the high- and low-risk groups (Figure 6A). GSEA revealed that several tumorigenesis pathways, including E2F targets, G2/M checkpoint, and KRAS signaling, were significantly enriched in high-risk group (Figure 6B). KEGG analysis showed that metabolism pathways (e.g., fatty acid metabolism) were enriched in the low-risk group (Figure 6C). Proliferation, angiogenesis, and genome instability are the main hallmarks of cancer (21). Therefore, we constructed gene sets related to base excision repair, cell cycle, cell cycle regulation, DNA damage repair 1, DNA damage repair 2, DNA replication, homologous recombination, mismatch repair, and nucleotide excision repair. GSVA analysis showed that these gene sets were significantly elevated in the high-risk group (Figure 6D). Surprisingly, angiogenesis was substantially enriched in patients with HCC with a low-risk score (Figure 6D). Sorafenib was approved as a first-line treatment of HCC due to its ability to mitigate angiogenesis and suppress tumor cell proliferation (22,23). Thus, we investigated whether the mRNASig could predict the sorafenib response of patients with HCC. In the sorafenib cohort (GSE109211) in which 67 patients with HCC accepted sorafenib treatment, we found significant therapeutic advantages of sorafenib in the low-risk group in comparison to the high-risk group (Figure 6E,6F).

Figure 6 The biological characteristics and sorafenib response of high- and low-risk patients with HCC. (A) DEGs between the 2 groups. (B,C) GSEA between the 2 groups in the Hallmark and KEGG gene sets. (D) GSVA showing the activation stations of special biological processes in the 2 groups. (E,F) The proportion of patients responding to sorafenib treatment in the low- and high-risk group. Responder/nonresponder: 42%/58% in the low-risk group and 21%/79% in the high-risk group. **, P<0.01; *** P<0.001; ****, P<0.0001. HCC, hepatocellular carcinoma; NS, nonsignificant; DEGs, differentially expressed genes; GSEA; gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSVA, gene set variation analysis.

The immune landscape in the high- and low-risk prognostic groups

HCC is characterized by a heterogeneous tumor immune environment (4). Using the ssGSEA method, we identified the infiltration abundances of 23 tumor-infiltrating lymphocytes (TILs) in each HCC sample. The high-risk group was characterized by an abundance of activated CD4 T cells, immunosuppressive cells, and type 2 T helper (Th2) cells and significantly low infiltration levels of activated CD8 T cells, eosinophils, macrophages, mast cells, NK cells, neutrophil, and Th1 cells (Figure 7A). Given the essential role of checkpoint immunotherapy in cancer, we further investigated the expression levels of immune checkpoints. The high-risk group demonstrated increased expression levels of inhibitory markers (including CTLA4, PDCD1, TIGIT, and TIM3) when compared to the low-risk group (Figure 7B).

Figure 7 Immune landscape of the high- and low-risk patients with HCC. (A) Immune infiltration of 23 TILs in the 2 groups. (B) CTLA4, PDCD1, TIGIT, and TIM3 gene expression differences between the low- and high-risk groups. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. TME, tumor-associated microenvironment; HCC, hepatocellular carcinoma.

TMB in the high- and low-risk groups

Given that TMB is a biomarker of cancer, we compared the somatic mutation of the high- and low-risk groups. The TMB in somatic cells in the high-risk groups was different from that in the low-risk group (Figure 8A,8B). Interestingly, we found that the TP53 mutation occurred at rate of 45% in high-risk group and only at rate of 15% in the low-risk group. The high-risk group showed higher TMB than did the low-risk group (Figure 8C). Moreover, we discovered that patients with HCC with a high-risk score and high TMB had a poorer prognosis than those with a low-risk score and low TMB (Figure 8D).

Figure 8 Association between mRNASig and TMB and the predictive performance of TMB for OS. (A,B) The differences of TMB between the high- and low-risk patients with HCC. (C) The comparison of TMB between the high- and low-risk groups. (D) Kaplan-Meier analysis of the patients with different risk scores and different TMB. **, P<0.01. mRNASig, messenger RNA signature; TMB, tumor mutational burden; OS, overall survival; HCC, hepatocellular carcinoma.

Discussion

HCC is a common cancer globally, with rapid progression, a variety of treatments, but an extremely poor prognosis. The risk of HCC is increased by several known risk factors, such as hepatic virus, chronic alcoholic abuse, and nonalcoholic steatohepatitis (24), the latter of which has recently been increasing the incidence of HCC in Western countries (25). The circadian rhythm, an evolutionarily conserved system, regulates many physiological functions, including proliferation, metabolism, and immune response (26,27). It is reported that chronic perturbations (e.g., sleep disruption, jet lag, disruption of feeding and fasting) of the circadian rhythm are associated with metabolic disease and increased risk of cancer, including nonalcoholic fatty liver disease, breast cancer, and HCC. Thus, chronotherapy based on the circadian rhythm, with its limited toxicity, has recently received heightened attention from oncologists (14). Nevertheless, clinical trials of chronotherapy have still not shown improved efficacy in cancer compared with conventionally timed therapies (28). Therefore, it is necessary to investigate the role of the chronic disruption of the circadian rhythm in the prognosis, metabolism, and tumor-associated microenvironment (TME) of HCC.

We investigated the heterogeneity in the expression, CNVs, and mutation patterns of core circadian clock genes between HCC samples and noncarcinoma samples. The expression of the core circadian clock genes was disrupted in HCC, with upregulation of RORC, CLOCK, and CRY2 and downregulation of PER1 and RORA. Using pancancer analysis, Ye and colleagues found that the expression of several clock genes is disrupted in various cancers (29). Interestingly, RORC, which is expressed in multiple tissues, including the liver, and is known to be a key regulator of various biological processes, exhibits low expression in breast cancer, melanoma, and bladder cancer (30-33). However, our study demonstrated that the expression of RORC was higher in HCC samples compared to normal liver tissue. It is important to note that gene expression can be influenced by transcription factors and copy number alterations (34). In our study, we observed that certain core clock genes with copy number amplification exhibited a higher transcriptional expression. Specifically, RORC, which shows greater copy number amplification, had higher expression levels in HCC. These findings point to the important contribution of CNV in disrupting circadian clock gene expression in HCC.

To develop novel signatures to predict the prognosis of patients with HCC, we used LASSO analysis and stepwise multivariate Cox regression analysis on circadian clock genes to construct the mRNASig, which included NPAS2, NR1D1, PER1, RORC, and TIMELESS. The excellent predictive performance of mRNASig was confirmed in the validation cohort. Compared to people with low-risk scores, those with high-risk scores had significantly poorer prognosis. Furthermore, the mRNASig could potentially evaluate the different prognoses of the high- and low-risk groups classified according to subclinical parameters, such as age, gender, grade, clinical stage, and Child-Pugh score. The mRNASig was an independent prognostic factor in HCC and markedly increased the predictive capability of clinical models. Therefore, we constructed a nomogram combining the risk score and prognostic factors to predict OS and facilitate clinical decision-making.

To elucidate the underlying molecular mechanisms of the mRNASig and their potential biological function, KEGG enrichment analysis, GSEA, and GSVA were performed. As expected, the high-risk group was associated with malignant pathway, oncogenesis, proliferation, and genome instability. Clock genes have been found to exhibit strong associations with the activation or inhibition of various oncogenic pathways, including apoptosis, cell cycle, DNA damage response, and epithelial-mesenchymal transition (29). The high-risk group was also associated with metabolism. Some experiments showed that NPAS2 has a significant role in the reprogramming of glucose metabolism (35). Mechanistically, HIF-1α, a direct transcriptional target of NPAS2, upregulates glycolytic genes and downregulates mitochondrial biogenesis in HCC cells. Interestingly, we also noticed that angiogenesis was highly enriched in patients with a low-risk score and that these patients were sensitive to sorafenib. Studies have revealed that sorafenib inhibits HCC via antiangiogenesis (36-38). These findings suggest that the strategy of combined therapies including chronotherapy or targeting circadian clock genes and sorafenib may improve the therapeutic response in patients with HCC.

HCC forms an immunosuppressive microenvironment through several mechanisms (39). The immunosuppressive microenvironment mainly contains an accumulation of cells with negative regulatory immune activity [e.g., regulatory T cells (Tregs), tumor-associated macrophages (TAMs), and Th2 cells] and a high expression of coinhibitory lymphocyte signals and tolerogenic enzymes (e.g., PDCD1, CTLA-4, and IDO) (40). Notably, the Treg cells and TAMs of cancers highly express PDCD1, inhibiting the function of effector T cells and promoting T-cell exhaustion (41,42). We thus speculated that the high- and low-risk groups would display a unique immune landscape. Our study revealed that a poor infiltration of activated CD8 T cells and NK cells and a higher abundance of Th2 cells were observed in the high-risk group. During the initiation and progression of HCC, the percentage and absolute number of liver NK cells were greatly decreased (43). Meanwhile, the high-risk group exhibited a high expression of CTLA4, PDCD1, TIGIT, and TIM3, suggesting that patients with a high-risk score may be more sensitive to checkpoint immunotherapy. Additionally, TMB is a biomarker for predicting clinical response to immunotherapy in some cancers, and a higher TMB is associated with better OS (44). Mechanistically, higher TMB means a higher number of somatic mutations in cancer, and a greater number of neoantigens increases the likelihood of T-cell recognition, thus improving the antitumor response (45). In this study, the mRNASig was significantly correlated with a higher TMB. Moreover, patients with a high-risk score and high TMB had a poorer prognosis. These finding suggest that the high-risk group may be more sensitive to immunotherapy.

However, there are a few limitations to our study that should be mentioned. We examined the expression, CNVs, and mutation of core circadian clock genes in HCC and developed the mRNASig based on the genes to predict the OS of patients with HCC using only bioinformatics technology. The potential relationship between mRNASig and sorafenib treatment has not been well investigated. Our hypothesis, supported by transcriptomic evidence, necessitates further comprehensive investigations encompassing both in vivo and in vitro studies to ascertain its validity. Despite these limitations, the outcomes of this study were attained with meticulous data analysis, ensuring the accuracy of the results. Our findings point to a new research direction that can advance our comprehension of the underlying mechanisms linking circadian rhythm and HCC.

Conclusions

We explored the expression, CNVs, and mutations of core circadian clock genes of HCC and developed and validated the mRNASig based on five genes (NPAS2, NR1D1, PER1, RORC, and TIMELESS) which demonstrated an independent prognostic significance for HCC. Our study provides a novel nomogram combining mRNASig and clinical characters for predicting the OS of patients with HCC. Moreover, we also found that the high-risk group was associated with a variety malignant processes, including in pathways, metabolism, and immunosuppressive environments. Further enhancing the molecular comprehension of the intricate interplay between the circadian clock and its implications on carcinogenesis, metabolism and the TME of HCC, might help to find fundamental approaches in the treatment of HCC. Finally, we found that the mRNASig based on core circadian clock genes is a potential therapeutic strategy for identifying the patients more sensitive to sorafenib treatment.

Acknowledgments

Funding: This study was supported by grants from the Science and Technology Program of Guangdong Province (No. 2017ZC0222) and the Guangzhou Science and Technology Research Program (Nos. 202102010251, 201803010065, and 201704020215).

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-217/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-217/prf

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-217/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Yang JD, Hainaut P, Gores GJ, et al. A global view of hepatocellular carcinoma: trends, risk, prevention and management. Nat Rev Gastroenterol Hepatol 2019;16:589-604. [Crossref] [PubMed]
2. Vogel A, Meyer T, Sapisochin G, et al. Hepatocellular carcinoma. Lancet 2022;400:1345-62. [Crossref] [PubMed]
3. Sidali S, Trépo E, Sutter O, et al. New concepts in the treatment of hepatocellular carcinoma. United European Gastroenterol J 2022;10:765-74. [Crossref] [PubMed]
4. Llovet JM, Castet F, Heikenwalder M, et al. Immunotherapies for hepatocellular carcinoma. Nat Rev Clin Oncol 2022;19:151-72. [Crossref] [PubMed]
5. Patke A, Young MW, Axelrod S. Molecular mechanisms and physiological importance of circadian rhythms. Nat Rev Mol Cell Biol 2020;21:67-84. [Crossref] [PubMed]
6. Jiang Y, Shen X, Fasae MB, et al. The Expression and Function of Circadian Rhythm Genes in Hepatocellular Carcinoma. Oxid Med Cell Longev 2021;2021:4044606. [Crossref] [PubMed]
7. Fu L, Kettner NM. The circadian clock in cancer development and therapy. Prog Mol Biol Transl Sci 2013;119:221-82. [Crossref] [PubMed]
8. Zhang R, Lahens NF, Ballance HI, et al. A circadian gene expression atlas in mammals: implications for biology and medicine. Proc Natl Acad Sci U S A 2014;111:16219-24. [Crossref] [PubMed]
9. Pariollaud M, Lamia KA. Cancer in the Fourth Dimension: What Is the Impact of Circadian Disruption? Cancer Discov 2020;10:1455-64. [Crossref] [PubMed]
10. Reinke H, Asher G. Crosstalk between metabolism and circadian clocks. Nat Rev Mol Cell Biol 2019;20:227-41. [Crossref] [PubMed]
11. Papagiannakopoulos T, Bauer MR, Davidson SM, et al. Circadian Rhythm Disruption Promotes Lung Tumorigenesis. Cell Metab 2016;24:324-31. [Crossref] [PubMed]
12. Filipski E, Subramanian P, Carrière J, et al. Circadian disruption accelerates liver carcinogenesis in mice. Mutat Res 2009;680:95-105. [Crossref] [PubMed]
13. Crespo M, Leiva M, Sabio G. Circadian Clock and Liver Cancer. Cancers (Basel) 2021;13:3631. [Crossref] [PubMed]
14. Ruan W, Yuan X, Eltzschig HK. Circadian rhythm as a therapeutic target. Nat Rev Drug Discov 2021;20:287-307. [Crossref] [PubMed]
15. Lévi F. Daytime versus evening infusions of immune checkpoint inhibitors. Lancet Oncol 2021;22:1648-50. [Crossref] [PubMed]
16. Wagner GP, Kin K, Lynch VJ. Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory Biosci 2012;131:281-5. [Crossref] [PubMed]
17. Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw 2010;33:1-22. [Crossref] [PubMed]
18. Mariathasan S, Turley SJ, Nickles D, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 2018;554:544-8. [Crossref] [PubMed]
19. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
20. Charoentong P, Finotello F, Angelova M, et al. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep 2017;18:248-62. [Crossref] [PubMed]
21. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell 2011;144:646-74. [Crossref] [PubMed]
22. Llovet JM, Ricci S, Mazzaferro V, et al. Sorafenib in advanced hepatocellular carcinoma. N Engl J Med 2008;359:378-90. [Crossref] [PubMed]
23. Wilhelm SM, Carter C, Tang L, et al. BAY 43-9006 exhibits broad spectrum oral antitumor activity and targets the RAF/MEK/ERK pathway and receptor tyrosine kinases involved in tumor progression and angiogenesis. Cancer Res 2004;64:7099-109. [Crossref] [PubMed]
24. Llovet JM, Kelley RK, Villanueva A, et al. Hepatocellular carcinoma. Nat Rev Dis Primers 2021;7:6. [Crossref] [PubMed]
25. Anstee QM, Reeves HL, Kotsiliti E, et al. From NASH to HCC: current concepts and future challenges. Nat Rev Gastroenterol Hepatol 2019;16:411-28. [Crossref] [PubMed]
26. Sahar S, Sassone-Corsi P. Metabolism and cancer: the circadian clock connection. Nat Rev Cancer 2009;9:886-96. [Crossref] [PubMed]
27. Wang C, Lutes LK, Barnoud C, et al. The circadian immune system. Sci Immunol 2022;7:eabm2465. [Crossref] [PubMed]
28. Sancar A, Van Gelder RN. Clocks, cancer, and chronochemotherapy. Science 2021;371:eabb0738. [Crossref] [PubMed]
29. Ye Y, Xiang Y, Ozguc FM, et al. The Genomic Landscape and Pharmacogenomic Interactions of Clock Genes in Cancer Chronotherapy. Cell Syst 2018;6:314-328.e2. [Crossref] [PubMed]
30. Jetten AM. Retinoid-related orphan receptors (RORs): critical roles in development, immunity, circadian rhythm, and cellular metabolism. Nucl Recept Signal 2009;7:e003. [Crossref] [PubMed]
31. Oh TG, Bailey P, Dray E, et al. PRMT2 and RORγ expression are associated with breast cancer survival outcomes. Mol Endocrinol 2014;28:1166-85. [Crossref] [PubMed]
32. Brożyna AA, Jóźwicki W, Skobowiat C, et al. RORα and RORγ expression inversely correlates with human melanoma progression. Oncotarget 2016;7:63261-82. [Crossref] [PubMed]
33. Cao D, Qi Z, Pang Y, et al. Retinoic Acid-Related Orphan Receptor C Regulates Proliferation, Glycolysis, and Chemoresistance via the PD-L1/ITGB6/STAT3 Signaling Axis in Bladder Cancer. Cancer Res 2019;79:2604-18. [Crossref] [PubMed]
34. PCAWG Transcriptome Core Group. Genomic basis for RNA alterations in cancer. Nature 2020;578:129-36. [Crossref] [PubMed]
35. Yuan P, Yang T, Mu J, et al. Circadian clock gene NPAS2 promotes reprogramming of glucose metabolism in hepatocellular carcinoma cells. Cancer Lett 2020;469:498-509. [Crossref] [PubMed]
36. Liu L, Cao Y, Chen C, et al. Sorafenib blocks the RAF/MEK/ERK pathway, inhibits tumor angiogenesis, and induces tumor cell apoptosis in hepatocellular carcinoma model PLC/PRF/5. Cancer Res 2006;66:11851-8. [Crossref] [PubMed]
37. Xu M, Zheng YL, Xie XY, et al. Sorafenib blocks the HIF-1α/VEGFA pathway, inhibits tumor invasion, and induces apoptosis in hepatoma cells. DNA Cell Biol 2014;33:275-81. [Crossref] [PubMed]
38. Morse MA, Sun W, Kim R, et al. The Role of Angiogenesis in Hepatocellular Carcinoma. Clin Cancer Res 2019;25:912-20. [Crossref] [PubMed]
39. Sangro B, Sarobe P, Hervás-Stubbs S, et al. Advances in immunotherapy for hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol 2021;18:525-43. [Crossref] [PubMed]
40. Prieto J, Melero I, Sangro B. Immunological landscape and immunotherapy of hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol 2015;12:681-700. [Crossref] [PubMed]
41. McLane LM, Abdel-Hakeem MS, Wherry EJ. CD8 T Cell Exhaustion During Chronic Viral Infection and Cancer. Annu Rev Immunol 2019;37:457-95. [Crossref] [PubMed]
42. Togashi Y, Shitara K, Nishikawa H. Regulatory T cells in cancer immunosuppression - implications for anticancer therapy. Nat Rev Clin Oncol 2019;16:356-71. [Crossref] [PubMed]
43. Sun C, Sun H, Zhang C, et al. NK cell receptor imbalance and NK cell dysfunction in HBV infection and hepatocellular carcinoma. Cell Mol Immunol 2015;12:292-302. [Crossref] [PubMed]
44. Samstein RM, Lee CH, Shoushtari AN, et al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet 2019;51:202-6. [Crossref] [PubMed]
45. Jardim DL, Goodman A, de Melo Gagliato D, et al. The Challenges of Tumor Mutational Burden as an Immunotherapy Biomarker. Cancer Cell 2021;39:154-73. [Crossref] [PubMed]