Circadian Regulator-Mediated Molecular Subtypes Depict the Features of Tumor Microenvironment and Indicate Prognosis in Head and Neck Squamous Cell Carcinoma

Aye, Ling; Wang, Zhanying; Chen, Fanghua; Xiong, Yujun; Zhou, Jiaying; Wu, Feizhen; Hu, Li; Wang, Dehui

doi:https://doi.org/10.1155/2023/9946911

Journal of Immunology Research

On this page

Abstract Introduction Materials and Methods Results and Discussion Discussion Conclusions Abbreviations Data Availability Ethical Approval Conflicts of Interest Authors’ Contributions Acknowledgments Supplementary Materials References Copyright Related Articles

Special Issue

Roles of Infiltrating Immune Cells in Cancer

View this Special Issue

Research Article | Open Access

Volume 2023 | Article ID 9946911 | https://doi.org/10.1155/2023/9946911

Circadian Regulator-Mediated Molecular Subtypes Depict the Features of Tumor Microenvironment and Indicate Prognosis in Head and Neck Squamous Cell Carcinoma

Ling Aye,¹Zhanying Wang,²Fanghua Chen,³Yujun Xiong,⁴Jiaying Zhou,⁵Feizhen Wu,^6,7Li Hu ,¹and Dehui Wang¹

Academic Editor: Xiangyang Yu

Received12 Aug 2022

Accepted07 Sept 2022

Published12 Jun 2023

Abstract

Introduction. Circadian rhythm is involved in multiple biological activities and implicated in cancer development. However, the role of circadian rhythm in head and neck squamous cell carcinoma (HNSCC) has not been fully interpreted yet. Herein, the present study set out to explore the significance of circadian regulator genes (CRGs) in HNSCC. Materials and Methods. The molecular landscape and clinical significance of 13 CRGs in HNSCC were explored based on The Cancer Genome Atlas (TCGA). The biological functions of PER3, a key CRG, were validated by cellular experiments. The correlation of CRGs with microenvironment, pathway activities, and prognosis was determined by bioinformatic algorithms. A novel circadian score was introduced to evaluate the circadian modification pattern of each patient and further validated in an independent cohort from the Gene Expression Omnibus (GEO) dataset. Results. CRGs presented high heterogeneity in HNSCC at both genomic and transcriptomic levels. Specifically, PER3 indicated a better prognosis and inhibited HNSCC cell proliferation. Moreover, HNSCC tissues displayed three circadian regulator patterns with distinct clinical outcomes, transcriptomic characteristics, and microenvironment features. Circadian score was an independent risk factor and exhibited excellent predictive efficiency in both the training cohort from the TCGA database and the validation cohort from the GEO database. Conclusions. CRGs played an indispensable role in HNSCC development. An in-depth exploration of circadian rhythm would improve the understanding of HNSCC carcinogenesis and confer novel insights for future clinical practices.

1. Introduction

Head and neck squamous cell carcinoma (HNSCC) currently stands as the sixth most common cancer worldwide [1, 2], representing approximately 6% of all new cancer cases [3]. Despite the inspiring efficacy of surgery, chemotherapy, radiotherapy, and immunotherapy in early diagnosed cases, the clinical outcomes and quality of life of HNSCC patients with advanced HNSCC are still far from satisfaction, which is largely due to the inadequate understanding of HNSCC pathogenesis [4–6].

The recent advent of integrative multiomics technologies has manifested the importance of specific molecular pathways in carcinogenesis, among which circadian rhythm was deemed a promising target for cancer therapy [7]. Circadian rhythm refers to the periodic changes of organisms over time, including behavioral patterns and morphological structures, which plays an important role in a variety of physiological and biochemical processes [8]. The material basis of this rhythm is the molecular chronometer, also known as the circadian clock [9]. The circadian clock is regulated by a cell-autonomous and self-sustained oscillation mechanism, which is consisted of three elements: input pathways for entrainment, the pivot oscillator generating the rhythmicity, and output pathways for controlling rhythmic behaviors [10]. The input pathways, usually composed of external environmental factors such as light and temperature, are processed by oscillators to maintain the 24-hour cycle of an organism and ultimately affect the output pathways, such as body temperature and heart rate [11]. The oscillator is essentially a transcription-translation feedback loop (TTFL) [12]. In the canonical model for the mammalian circadian clock, CLOCK-BMAL1 complex activates the transcription of Cryptochrome (CRY) and Period (PER) genes by binding to their promoters during the day [11]. The resulting protein, reversely, dissociates the combination of CLOCK-BMAL1 and inhibits CLOCK-BMAL1-mediated transcriptional activation during the night [13, 14]. Other molecules, such as basic helix-loop-helix-PER-ARNT-SIM (bHLH-PAS) [15], casein kinase 1 (CSNK1) [16], and REV-ERBα/β (NR1D) [17], also participate in the composition of autonomous feedback loops, respectively.

The oscillator is conserved across species and exists in most body cells [18], allowing them to manipulate various biological processes and, therefore, participate in numerous diseases, including diabetes [19], cardiovascular diseases [20], and cancers [21]. Particularly, the disorder of circadian rhythm-related genes (CRGs) has been observed in many types of cancers, including breast cancer [22], lung cancer [23], and glioma [24]. As central regulators of multiple biological processes, CRGs can directly alter cancer cell behaviors by coordinating oncogenic or tumor-suppressor signaling pathways. For example, in vivo screening of key transcription factors shows that CLOCK and BMAL1 are potent regulators of stemness-related pathways and required for the proliferation of acute myeloid leukemia cells [25]. On the other hand, CRGs are indirectly implicated in cancer development by affecting other factors like tumor microenvironment (TME) [26]. CLOCK has been documented to recruit immune-suppressive microglia by directly regulating OLFML3 secretion, ultimately leading to glioblastoma progression [27]. These studies jointly suggest that circadian rhythm plays a critical role in cancer development and progression. However, the molecular landscape and clinical significance of CRGs in HNSCC remain ambiguous, which hamper the understanding of the HNSCC pathogenesis and the development of effective therapies.

Herein, we integrated the genome and transcriptome information of head and neck squamous cell carcinoma (HNSCC) based on The Cancer Genome Atlas (TCGA) to comprehensively evaluate possible circadian regulator patterns (CRPs) in HNSCC. The interplay of CRPs with different TME and signal pathways was determined using novel bioinformatic algorithms. Moreover, we identified 341 differentially expressed genes (DEGs) from three CRPs, among which 67 DEGs were correlated with overall survival (OS) and termed circadian signature genes (CSGs). Through combing machine learning methods, we established and validated a circadian score system based on CSGs to describe the circadian rhythm patterns of individual HNSCC patient. Our results demonstrated that circadian score might serve as a surrogate biomarker for indicating clinical outcomes of HNSCC patients.

2. Materials and Methods

2.1. Collection of HNSCC Datasets and Preprocessing

The main concept and workflow of our study are shown in Supplementary Figure 1. The genetic alteration, mRNA expression profile, and corresponding clinical information of the TCGA cohort were downloaded from the Cancer Genomics Browser of University of California Santa Cruz (UCSC). The expression profile and corresponding clinical information of GSE65858 were downloaded from the Gene Expression Omnibus (GEO) database [28]. The clinical characteristics of these patients are summarized in Table 1.

2.2. Unsupervised Clustering for CRGs

The mRNA sequencing data were transformed into transcripts per kilobase million (TPM) values and normalized for cluster analysis. Unsupervised clustering methods were used to identify different circadian rhythm-related subtypes by using the package ConsensusClusterPlus. A total of 13 CRGs (CSNK1E, PER1, RHYTHM, CRY1, CRY2, BHLHE41, PER2, PER3, ARNTL, NPAS2, BHLHE40, CSNK1D, and NR1D1) were selected based on KEGG_CIRCADIAN_RHYTHM_MAMMAL gene set from Kyoto Encyclopedia of Genes and Genomes (KEGG).

2.3. DEG Analysis and Function Enrichment Analysis

DEGs were identified at a cutoff of value < 0.05 with the Limma package [29]. Enrichment analysis was performed by using the clusterProfiler package with default parameters [30]. KEGG [31], gene set enrichment analysis (GSEA) [32], Gene Ontology (GO) enrichment analysis, and gene set variation analysis (GSVA) [33] were employed as well. Specifically, hallmark gene sets were downloaded from the Molecular Signature Database (MSigDB), and immune-related gene sets were downloaded from the study of Mariathasan et al. [34].

2.4. Estimation of Immune Cell Infiltration

The immune infiltration was analyzed via ESTIMATE [35], CIBERSORT [36], and single-sample gene set enrichment analysis (ssGSEA) algorithm. The gene sets utilized in ssGSEA were downloaded from the study of Charoentong et al. [37].

2.5. Generation of Circadian Score

We determined the prognostic value for DEGs among three CRPs by the univariable Cox proportional hazard regression analysis. The DEGs with significant prognostic value () were termed circadian signature genes (CSGs). Then, the least shrinkage and selection operator (LASSO) and multivariable Cox model were utilized to construct circadian score. The circadian score was constructed using 6 CSGs (DEFB1, MTHFD2, KCTD11, POMP, RPL12, and TNFRSF12A) with a linear combination of their expression values. These inputs were weighted with the regression coefficients from the stepwise Cox regression analyses. The circadian score was shown as follows:

2.6. Cell Culture and Functional Assays

The HNSCC cell lines 6-10B and Tu 686 were purchased from FuHeng Cell Center (Shanghai, China). The cell lines mentioned above were cultured in complete medium: RPMI Medium 1640 (Gibco, USA) supplemented with 10% fetal calf serum (Gibco, USA), penicillin (100 units/ml) and streptomycin (100 μg/ml) at 37°C in a thermostatic incubator with 5% CO₂. All cells were validated using DNA short tandem repeat analysis each six months and tested for mycoplasma contamination using One-Step Mycoplasma Detection Kit (Yeason, Shanghai, China) each four weeks.

Functional assays were conducted in two groups: a negative group (sh-NC) and an shRNA knock-down group (sh-PER3). Cell proliferation was detected by the Cell Counting Kit-8 (CCK-8) and clone formation assays. In brief, cells were inoculated into 96-well plates (1000 cells/well). At each time point (1st, 2nd, 3rd, 4th, and 5th day), 10 μl of CCK-8 solution (Yeason, Shanghai, China) was added to the sextuplicate wells. The plates were incubated for 1.5 h, and the optical density values were detected. For clone formation assays, cells were seeded in a six-well plate (1000 cells/well) with the culture medium refreshed every 3 days for 2 weeks. Following the 2 weeks, the cells were washed with PBS, fixed, stained, and counted.

2.7. Virus Transfection

The vectors of plko-Puromycin-EGFP-shRNA-PER3 and plko-Puromycin-EGFP-NC were purchased from Ruoji Technology (Shanghai, China) and were transfected into 6-10B and Tu 686 cells. The sequences of short-hairpin RNA (shRNA) were as follows:

shRNA1: CCCUCGGAGAGACGCAAUAAA (F) and UUUAUUGCGUCUCUCCGAGGG (R)

shRNA2: AUGACCAUGAAGUUAUCAUUG (F) and CAAUGAUAACUUCAUGGUCAU (R)

shRNA3: CGACAGCCUCUUCUGCGAUAU (F) and AUAUCGCAGAAGAGGCUGUCG (R)

These three target sequences were mixed with the same proportion. After transfection, the cells with suitable fluorescence expression were screened with puromycin at a concentration of 4 μg/ml.

2.8. Real-Time Polymerase Chain Reaction (Real-Time PCR)

Total RNA was purified using Mini BEST Universal RNA extraction KIT (TaKaRa, Japan), and cDNA was synthesized using the Prime-Script RT Master Mix (TaKaRa, Japan) according to the manufacturer’s instructions. Real-time PCR was performed using SYBR Green Realtime PCR Master Mix (Yeasen, China). Samples from each experiment were analyzed in triplicate. The primer sequences used in this study were as follows:

GAPDH: GGACTCATGACCACAGTCCA (F) and CCAGTAGAGGCAGGGATGAT (R)

PER3: AGCGTTCAAGCAAACAGTGAG (F) and CAAGCAGGTCAACAAAGTGAGA (R)

2.9. In Vivo Assay

Five-week-old athymic BALB/c female nude mice were purchased from the Jihui Laboratory Animal Care Co., Ltd. (Shanghai, China) and maintained in pathogen-free environment, in accordance with the stated guidelines of 3Rs (replacement, reduction, and refinement). These mice were randomly divided into two subgroups (/each group): sh-NC group and sh-PER3 group. A total of cells were injected into the back of the mice to generate subcutaneous tumors. Tumor volume was measured regularly until the 24th day and calculated using the formula: . Twenty-four days after injection, the tumor specimens were surgically removed, weighted, and fixed. The animal experiments met the requirements of the Animal Ethics Committee of Shanghai Eye & ENT Hospital, Fudan University.

2.10. Statistical Analysis

The results were expressed as the . Parameter test or rank-sum test was used for comparisons between groups. The Bonferroni test was used for multiple comparisons. Categorical data were analyzed by the chi-square test or Fisher’s exact test. For survival analysis, the Kaplan-Meier method, log-rank test, and Cox regression analysis were used to test the prognostic value. All statistical tests were bilateral with value < 0.05 being statistically significant. R software (4.1.1) and GraphPad Prism 7 were used for data analyses.

3. Results and Discussion

3.1. The Molecular Landscape of CRGs in HNSCC

In this study, we identified 13 CRGs and summarized the incidence of somatic mutations of these CRGs based on the TCGA cohort. Genetic mutations of CRGs occurred in approximately 8.8% of HNSCC patients, among which CSNK1E and PER1 exhibited the highest mutation frequency (Figure 1(a)). The copy number variation (CNV) analysis revealed significant copy number deletions for BHLHE40 and PER2 and copy number amplifications mainly for BHLHE41 (Figure 1(b)). The univariable Cox analysis was used to assess the prognostic value of each CRG (Figure 1(c)). The result showed that the hazard ratio of PER3 was the lowest among 13 CRGs, suggesting its role as a protective factor in HNSCC. In addition, the differential expression patterns of CRGs between HNSCC tumors and adjacent normal tissues were also analyzed. The result manifested that five CRGs were highly expressed in tumors (Figure 1(d)), while opposite trends were observed for CRY2, PER2, and PER3. The aforementioned analyses suggested the high heterogeneity of genetic alteration and expressional landscape of CRGs in HNSCC, implying that dysregulation of CRGs was of great significance in HNSCC.

(a)

(b)

(c)

(d)

Figure 1

Genetic alteration, expression variation, and clinical significance of CRGs in HNSCC. (a) The mutation frequency of CRGs in TCGA HNSCC cohort. (b) The CNV frequency of CRGs in the TCGA HNSCC cohort. Gain, yellow; loss, blue. (c) The forest plot of the univariate Cox regression model depicting the prognostic value of CRGs in TCGA HNSCC cohort. Hazard : risk factors for survival. Hazard : protective factors for survival. (d) Expression level of CRGs between normal and tumor tissues in the TCGA HNSCC cohort. Ns: not significant; : ; : ; : .

3.2. PER3 Knockdown Promoted Proliferation of HNSCC Cell

Considering the significant prognostic value and the distinct expression pattern of PER3, PER3 was selected for deeper investigation. The Wilcoxon test showed that PER3 expression was downregulated in patients with higher T stage () and M stage (, supplementary Figure 2A-D). GSEA analysis manifested that the enrichment of proliferation-related pathways, including DNA replication and the ErbB pathway, was significantly increased in samples with low expression of PER3 (, Figures 2(a) and 2(b)). Thereafter, cellular experiments were performed to verify the effect of PER3 on the biological behaviors of HNSCC cells. Firstly, shRNA was employed to knock down PER3 expression in two typical HNSCC cell lines, and the transfection efficiency was detected by real-time PCR (Figures 2(c) and 2(f)). CCK-8 assay demonstrated that inhibition of PER3 increased HNSCC cell viability in vitro (Figures 2(d) and 2(g)), which was further confirmed by colony formation assay (Figures 2(e) and 2(h)). Subsequently, a subcutaneous xenograft model was established to validate the function of PER3 in vivo (Figure 2(i)). The tumor weight and volume were obviously higher in the sh-PER3 group than those in the sh-NC group (Figures 2(j) and 2(k)). In summary, these results collectively illustrated that PER3 downregulation correlated with proliferation-related pathways and indeed promoted HNSCC cell vitality.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

Figure 2

Downregulation of PER3 promoted HNSCC cell proliferation. (a, b) GSEA GO and GSEA KEGG analyses identified proliferation-related pathways enriched in the low PER3 expression samples. (c, f) The efficiency of PER3 transfection in Tu 686 and 6-10B cells was confirmed by real time-PCR. (d, g) Downregulation of PER3 in Tu 686 and 6-10B cells promoted cell proliferation as determined by CCK-8. (e, h) Colony formation assays showed that inhibition of PER3 in Tu 686 and 6-10B cells promoted cell proliferation. (i) Subcutaneous tumors obtained from nude mice implanted with the sh-NC or sh-PER3 cell lines. (j) The tumor volume is shown to be significantly increased in sh-PER3 group compared with that in the sh-NC group. (k) The tumor weight is significantly higher in sh-PER3 group than in the sh-NC group. Data represent of three independent experiments. : ; : ; : .

3.3. HNSCC Patients Exhibited Three Circadian Rhythm-Related Subtypes with Distinct Prognosis and TME

To explore the expression pattern of CRGs, consensus clustering analysis was conducted based on 13 CRGs to stratify HNSCC patients. Three clusters were achieved and termed CRPs, of which there were 198 cases in CRP-A, 170 cases in CRP-B, and 129 cases in CRP-C (Figure 3(a)). Obviously, CRGs had diverse expression levels among three CRP subtypes (Figure 3(b)). Log-rank analysis revealed that CRP-B presented a particularly prominent survival disadvantage over CRP-A and CRP-C (, Figure 3(c)). Chi-square analysis showed that CRP was closely correlated to T stage () and clinical stage (, Figure 3(d)). Meanwhile, immune infiltration status was also different in three CRPs, with CRP-B showing the lowest immune score (, Figure 3(e)). Specifically, the infiltration of CD8+ T cells, which participated in active antitumor responses and correlated with better prognosis [38], displayed the lowest level in CRP-B (Figure 3(f)). ssGSEA results also validated that CD8+ T cells and cytotoxic T lymphocytes (CTL) were downregulated in CRP-B (Figure 3(g)). Based on these findings, we confirmed that HNSCC samples exhibited three CRPs, characterized by different prognosis and immune infiltration profiles.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 3

HNSCC patients exhibited three circadian rhythm-related subtypes with distinct prognosis and TME. (a) Consensus clustering matrix of HNSCC patients for . (b) Expression of 13 CRGs in three CRPs. (c) Survival analysis for three CRP subtypes. (d) An alluvial diagram depicting the relationship between CRP and clinical features (; ; ; ; ; ; ; ; ; and ). (e) Boxplots depicted the differences in immune scores in the three clusters. (f, g) Different immune cell component among three CRP subtypes as determined by CIBERSORT and ssGSEA. Ns: not significant; : ; : ; : .

3.4. Transcriptomic Characteristics of Different CRP Subtypes

To elucidate the underlying mechanisms involved in different clinical outcomes and immune contextures among three CRP subtypes, potential pathways were explored using bioinformatic algorithms. GSEA analysis revealed that, compared with CRP-A and CRP-C, CRP-B was primarily enriched in stroma-related processes such as collagen binding and N-glycan biosynthesis (, Figures 4(a) and 4(b)). Previous research reported that stromal components could repress antitumor responses in various ways, such as shaping immune cell composition through the formation of geographic barriers and affecting immune cell viability via the release of soluble factors [39]. Therefore, we speculated that stromal components in TME inhibited CD8+ T cell-mediated cytotoxic responses, ultimately leading to the poor prognosis of CRP-B. By combing hallmark gene sets and immune infiltration-related gene sets, we further explored different signal activities by GSVA analysis. The results verified that CRP-B was characterized by the activation of stroma-related pathways, like epithelial-mesenchymal transition (EMT), and the repression of antitumor immune responses, including CD8 effecter, cytolytic activity, and interferon responses (Figures 4(c) and 4(d)). Besides, GSEA and GSVA jointly suggested that stemness-related pathways, such as Wnt and Myc pathways, might be involved in these biological processes (Figures 4(a) and 4(d)).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4

Transcriptomic characteristics of different CRP subtypes. GSEA GO (a) and GSEA KEGG (b) analyses demonstrated stroma- and stemness-related pathways enriched in the CRP-B vs. CRP-A and CRP-C. (c) Boxplot of several immune signatures for three CRPs. (d) GSVA enrichment analysis showing the activation states of possible pathways in different CRPs based on hallmark gene sets. (e) 341 circadian-related DEGs between three CRPs were shown in the Venn diagram. (f) Functional annotation for circadian signature genes using GO enrichment analysis. Ns: not significant; : ; : ; : .

Differential expression analysis was conducted to further explore the heterogeneity of different CRPs, and 341 DEGs were identified (Figure 4(e)). The univariable Cox analysis revealed that 67 DEGs were associated with overall survival (OS) () and termed CSGs. GO analysis exhibited that CSGs mainly centered around mRNA and protein regulation, in line with the biological essence of circadian clock as a TTFL (Figure 4(f)). To further ascertain the possible role of CSGs, unsupervised clustering analysis was conducted based on the obtained CSGs, which yielded three circadian genomic clusters (CGCs), namely, CGC1-3. Similar to CRP, three CGCs distinguished different clinical characteristics, OS, and immune infiltration level, particularly T cell and CD8+ T cell abundance (supplementary Figure 3). Notably, we found that PER3 expression varied from different CGC subtypes and CRP subtypes (supplementary Figure 2E and F), implying that PER3 might be implicated in different clinical outcomes of these subtypes. Altogether, these results depicted the different transcriptomic features among CRPs, which might be involved in the different TME landscape and clinical outcomes.

3.5. Construction of Circadian Score and Exploration of Its Clinical Relevance

Although our aforementioned findings elaborated on the association among CRPs, prognosis, pathway activities, and TME, these analyses were mainly based on the patient population and could not provide a specific circadian rhythm pattern for an individual patient. Therefore, a circadian score system based on CSGs was introduced. To optimize this model, mainly to minimize multicollinearity, LASSO analysis was utilized with six CSGs ultimately selected for the construction of circadian score (Figures 5(a) and 5(b)). Patients were arranged into high-risk or low-risk groups according to the median cutoff value of circadian score (Figure 5(c)). Log-rank analysis revealed that the OS in the high-risk group was significantly shorter than that in the low-risk group (, Figure 5(d)). The 1-, 2-, and 3-year AUC values of the circadian score were 0.613, 0.650, and 0. 675, respectively (Figure 5(e)). Furthermore, the reliability of circadian score was validated in different clinical subgroups and a validation cohort. Subgroup analysis indicated that circadian score was significantly correlated with poor OS in all subgroups except the M-high subgroup (supplementary Figure 4), which might be attributed to the limited sample size of the M-high subgroup (15 out of 497 patients, Table 1). Circadian score also demonstrated extraordinary performance in another independent HNSCC cohort from the GEO database (Figures 5(f) and 5(g)). The correlation of circadian score with clinicopathological parameters, CRP, and CGC was also explored to clarify the features of circadian score. Chi-square analysis displayed that the circadian risk group was significantly correlated with M stage, survival status, CRP, and CGC (Figure 5(h)). The Kruskal-Wallis test revealed that patients in CRP-B showed the highest circadian score (Figure 5(i)), consistent with their worst OS. For CGC, CGC-2 had the lowest circadian score in contrast to other CGC subtypes (Figure 5(j)). The above results jointly revealed that circadian score was an effective biomarker in predicting HNSCC prognosis.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

Figure 5

Construction of the circadian score and exploration of its clinical relevance. (a) LASSO coefficient profiles of the 67 genes selected. The dotted vertical line was drawn at the optimal scores by minimum criteria. (b) Six genes were selected using LASSO Cox regression analysis. The dotted vertical line was drawn at the optimal scores by minimum criteria. (c) The risk scores for 497 patients of the TCGA cohort, the survival of each patient, and gene expression distribution in the TCGA cohort. (d) Kaplan-Meier curves of the OS in the TCGA cohort on the basis of the circadian score. (e) ROC analysis of the circadian score for prediction of survival at 1, 2, and 3 years in the TCGA cohort. (f) Kaplan-Meier curves of the OS in the GEO cohort on the basis of the circadian score. (g) ROC analysis of the circadian score for prediction of survival at 1, 2, and 3 years in the GEO cohort. (h) An alluvial diagram showing the relationship of circadian score with clinical characteristics, CRP, and CGC (; ; ; ; ; ; ; ; ; and ). (i) Violin plot of the circadian score for three CRPs. (j) Violin plot of the circadian score for three CGCs. Ns: not significant; : ; : ; : .

3.6. Establishment of a Circadian Score-Based Nomogram

The univariable and multivariable Cox regression analyses were performed to assess the independence of circadian score. The univariable Cox regression analysis demonstrated that age (), gender (), N (), M (), and circadian score () were associated with the OS of HNSCC (Figure 6(a)). The multivariable Cox regression analysis revealed that the circadian score was still significantly correlated to OS () and an independent prognostic factor (Figure 6(b)). To facilitate the application of our signature in future clinical practice, a nomogram was built based on the results of the univariable Cox analysis (Figure 6(c)). Besides, calibration curves showed an excellent consistency between the actual and nomogram-predicted probabilities of 1-year, 2-year, and 3-year OS (Figure 6(d)). Together, these results emphasized the clinical significance of the circadian score and the accuracy of the novel circadian score-based nomogram.

(a)

(b)

(c)

(d)

4. Discussion

Mounting evidence has reported the critical role of circadian clock in mediating immunity, inflammation, and tumorigenesis [40, 41]. However, the influence of CRGs in HNSCC has not been comprehensively delineated.

Here, we systematically evaluated the role of CRGs in determining different clinical outcomes, pathway activities, and TME features in HNSCC (Figure 7). Specifically, PER3 expression was significantly correlated to a better prognosis among 13 CRGs. Gene enrichment analysis indicated that proliferation-related pathways were enriched in the PER3-low expression samples, suggesting the potential of PER3 in inhibiting cancer cell proliferation. By combining cellular experiments, we found that PER3 depletion could promote HNSCC cell proliferation, consistent with its prognostic value.

Based on 13 CRGs, this study identified three distinct CRPs in HNSCC. Compared to the other two subtypes, samples in the CRP-B subtype exhibited immunosuppressive microenvironment and had the worst OS rate, in line with the critical role of immune surveillance in HNSCC progression [42]. In addition, gene enrichment analysis displayed that CRP-B shared a feature of stroma- and stemness-related pathways, such as EMT and Wnt pathways. It is widely recognized that cancer cells maintain a dynamic interaction with their surrounding stromal components to affect the local immune contexture mutually. For instance, Evan and Vousden [43] discovered that mesothelial cell-derived fibroblasts could induce naive CD4+ T cells into regulatory T cells. Meanwhile, stemness is reported to be associated with immunologically “cold” tumors. Specifically, Wnt signaling pathways can dampen the functions of dendritic cells and effector T cells via secretion of IL-10 and IL-12, leading to an immunosuppressive microenvironment [44, 45]. Consistent with these previous findings, we speculated that activation of stemness- and stroma-related pathways caused an immunosuppressive TME, eventually resulting in the worst prognosis of CRP-B.

Furthermore, 341 DEGs were identified among three CRPs. The univariable Cox regression analysis revealed that 67 out of 341 DEGs were correlated with the OS and termed CSGs. GO analysis indicated that CSGs revolved around transcription and translation regulation, highlighting the central role of the circadian clock in mediating various biological activities. Based on CSGs, three CGCs were identified and found to be closely correlated to TME and OS. These results collectively unraveled that three circadian-related subtypes did exist in HNSCC, in which CRGs and CSGs were vital for the clinical outcomes of HNSCC by shaping TME.

The complexity of cancer has long constituted a main obstacle to precision medicine [46]. To overcome heterogenicity and provide accurate prediction, we innovatively established a six CSG-based circadian score system by performing LASSO and Cox regression analyses. Patients with higher risk scores tended to have shorter OS than those with lower risk scores in both the training cohort from the TCGA dataset and the validation cohort from the GEO dataset. Moreover, we constructed a quantitative nomogram by integrating the circadian score and other clinical features, which exhibited extraordinary reliability as examined by calibration plots. Notably, these six CSGs assume essential roles in numerous human cancers and have been reported as crucial regulators of various immune responses. For instance, MTHFD2 can induce immune escape of cancer cells by promoting the expression of basal and IFN-γ-stimulated PD-L1 [47]. Inhibition of POMP can result in the accumulation of cyclin-dependent kinase inhibitors and thus hamper cancer cell proliferation in vitro and in vivo [48].

However, our study still has some limitations. Since this is a retrospective analysis, there may be unavoidable bias in our results. Future prospective studies are warranted to validate the reliability of our conclusion. On the other hand, although this study has depicted the relationship between circadian rhythm and TME, the specific mechanisms remain to be further investigated.

Evidence has reported the critical role of circadian clock in mediating immunity, inflammation, and tumorigenesis [40, 41]. However, the influence of CRGs in HNSCC has not been comprehensively delineated.

Here, we systematically evaluated the molecular landscape and prognostic value of 13 CRGs and discovered that PER3 expression was significantly correlated to a better prognosis (Figure 7). Gene enrichment analysis indicated that proliferation-related pathways were enriched in the PER3-low expression samples, suggesting the potential of PER3 in inhibiting cancer cell proliferation. By combining cellular experiments, we found that PER3 depletion could promote HNSCC cell proliferation, consistent with its prognostic value.

5. Conclusions

In conclusion, this study investigated the molecular landscape and clinical significance of CRGs in HNSCC. By mediating various pathways and shaping distinct TME, CRGs play an indispensable role in the heterogeneity of HNSCC. A comprehensive investigation of circadian rhythm of individual patient would contribute to the understanding of cancer development and guide more personalized treatment strategies.

Abbreviations

AUC:	Area under curve
CCK-8:	Cell Counting Kit-8
CGC:	Circadian genomic cluster
CNV:	Copy number variation
CRG:	Circadian regulator genes
CRP:	Circadian regulator pattern
CRY:	Cryptochrome
CSG:	Circadian signature genes
CTL:	Cytotoxic T lymphocytes
DEG:	Differentially expressed gene
EMT:	Epithelial-mesenchymal transition
GEO:	Gene Expression Omnibus
GO:	Gene Ontology
GSEA:	Gene set enrichment analysis
GSVA:	Gene set variation analysis
HNSCC:	Head and neck squamous cell carcinoma
KEGG:	Kyoto Encyclopedia of Genes and Genomes
LASSO:	Least shrinkage and selection operator
NC:	Negative control
OS:	Overall survival
PER:	Period
real-time PCR:	Real-time polymerase chain reaction
ROC:	Receiver operating characteristic
shRNA:	Short-hairpin RNA
ssGSEA:	Single sample GSEA
TCGA:	The Cancer Genome Atlas
TME:	Tumor microenvironment
TPM:	Transcripts per kilobase million
TTFL:	Transcription-translation feedback loop
UCSC:	California Santa Cruz.

Data Availability

Genomic and transcriptomic data and patients’ clinical data are available at UCSC Xena (https://xena.ucsc.edu/) and GEO database (https://www.ncbi.nlm.nih.gov/geo/).

Ethical Approval

All animal experiments were approved by the Laboratory Animal Welfare and Ethics Committee of Shanghai Eye & ENT Hospital, Fudan University.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Authors’ Contributions

FW, LH, and DW contributed to conceptualization, supervision, and project administration. LA, ZW, FC, YX, and JZ were responsible for methodology and investigation. LA was responsible for the software, formal analysis, and data curation. ZW and FC validated the study. LA, ZW, FC, and YX wrote the original draft. LA, FW, LH, and DW wrote, reviewed, and edited the manuscript. LA and ZW visualized the study. FW and DW contributed to funding acquisition. LA, ZW, FC, and YX contributed equally to this work and share the first authorship. All authors read and approved the final manuscript.

Acknowledgments

This work was supported by the Medical Research Data Center of Fudan University and grants from the Shanghai Shenkang Hospital Development Center (SHDC12018118), Clinical Research Plan of Shanghai Shenkang Hospital Development Center (SHDC2020CR2005A), Shanghai Municipal Science and Technology Major Project (Grant No. 2017SHZDZX01), National Key Research and Development Program of China (2018YFC1005004), and National Natural Science Foundation of China (32070653).

Supplementary Materials

Figure 1: workflow of our study. The workflow of our study and analytical pipeline. Figure 2: comparison of PER3 expression in different subgroups. Box plots of PER3 expression in different T (A), N (B), M (C), stage (D), CRP subtypes (E), and CGC subtypes (F). Figure 3: CGC exhibited different TME and prognosis. (A) Consensus clustering matrix of HNSCC patients based on CSGs (). (B) Expression level of CSGs in the TCGA cohort. (C) A Sankey diagram exhibiting the relationship of CGCs, CRP and clinical features (; ; ; ; ; ; ; ; ; and ). (D) OS of different CGC subtypes. (E) A boxplot of CD8+ T cell abundance for three CGCs. (F) Difference in immune signatures among different CGCs. (G) GSVA analysis among different CGCs. Figure 4: validation of the performance of circadian score in different subgroups. Kaplan-Meier curves of OS based on circadian score in different T (A), N (B), M (C), stage (D), grade (E), and gender (F). (Supplementary Materials)

References

C. R. Leemans, B. J. Braakhuis, and R. H. Brakenhoff, “The molecular biology of head and neck cancer,” Nature Reviews. Cancer, vol. 11, no. 1, pp. 9–22, 2011.
View at: Publisher Site | Google Scholar
R. L. Siegel, K. D. Miller, H. E. Fuchs, and A. Jemal, “Cancer statistics, 2021,” CA: a Cancer Journal for Clinicians, vol. 71, no. 1, pp. 7–33, 2021.
View at: Publisher Site | Google Scholar
D. E. Johnson, B. Burtness, C. R. Leemans, V. W. Y. Lui, J. E. Bauman, and J. R. Grandis, “Head and neck squamous cell carcinoma,” Nature Reviews Disease Primers, vol. 6, no. 1, 2020.
View at: Publisher Site | Google Scholar
A. G. Sacco, R. Chen, F. P. Worden et al., “Pembrolizumab plus cetuximab in patients with recurrent or metastatic head and neck squamous cell carcinoma: an open-label, multi-arm, non-randomised, multicentre, phase 2 trial,” The Lancet Oncology, vol. 22, no. 6, pp. 883–892, 2021.
View at: Publisher Site | Google Scholar
T. Lu, Z. Zhang, J. Zhang et al., “CD73 in small extracellular vesicles derived from HNSCC defines tumour-associated immunosuppression mediated by macrophages in the microenvironment,” Journal of Extracellular Vesicles, vol. 11, no. 5, article e12218, 2022.
View at: Publisher Site | Google Scholar
J. L. Low, D. P. Lau, X. Zhang et al., “A chemical genetic screen identifies Aurora kinases as a therapeutic target in EGFR T790M negative, gefitinib-resistant head and neck squamous cell carcinoma (HNSCC),” eBioMedicine, vol. 64, article 103220, 2021.
View at: Publisher Site | Google Scholar
W. Ruan, X. Yuan, and H. K. Eltzschig, “Circadian rhythm as a therapeutic target,” Nature Reviews Drug Discovery, vol. 20, no. 4, pp. 287–307, 2021.
View at: Publisher Site | Google Scholar
N. W. Gentry, L. H. Ashbrook, Y. H. Fu, and L. J. Ptacek, “Human circadian variations,” Journal of Clinical Investigation, vol. 131, no. 16, 2021.
View at: Publisher Site | Google Scholar
J. Kovac, J. Husse, and H. Oster, “A time to fast, a time to feast: the crosstalk between metabolism and the circadian clock,” Molecules and Cells, vol. 28, no. 2, pp. 75–80, 2009.
View at: Publisher Site | Google Scholar
A. M. Zuber, G. Centeno, S. Pradervand et al., “Molecular clock is involved in predictive circadian adjustment of renal function,” Proceedings of the National Academy of Sciences, vol. 106, no. 38, pp. 16523–16528, 2009.
View at: Publisher Site | Google Scholar
K. Kinouchi and P. Sassone-Corsi, “Metabolic rivalry: circadian homeostasis and tumorigenesis,” Nature Reviews. Cancer, vol. 20, no. 11, pp. 645–661, 2020.
View at: Publisher Site | Google Scholar
C. L. Partch, C. B. Green, and J. S. Takahashi, “Molecular architecture of the mammalian circadian clock,” Trends in Cell Biology, vol. 24, no. 2, pp. 90–99, 2014.
View at: Publisher Site | Google Scholar
H. Cho, X. Zhao, M. Hatori et al., “Regulation of circadian behaviour and metabolism by REV-ERB-α and REV-ERB-β,” Nature, vol. 485, no. 7396, pp. 123–127, 2012.
View at: Publisher Site | Google Scholar
X. Cao, Y. Yang, C. P. Selby, Z. Liu, and A. Sancar, “Molecular mechanism of the repressive phase of the mammalian circadian clock,” Proceedings of the National Academy of Sciences of the United States of America, vol. 118, no. 2, 2021.
View at: Publisher Site | Google Scholar
D. C. Bersten, A. E. Sullivan, D. J. Peet, and M. L. Whitelaw, “bHLH-PAS proteins in cancer,” Nature Reviews. Cancer, vol. 13, no. 12, pp. 827–841, 2013.
View at: Publisher Site | Google Scholar
U. E. Nydegger, P. M. Escobar, L. Risch, M. Risch, and Z. Stanga, “Chronobiology and circadian rhythms establish a connection to diagnosis,” Diagnosis, vol. 1, no. 4, pp. 295–303, 2014.
View at: Publisher Site | Google Scholar
G. Sulli, A. Rommel, X. Wang et al., “Pharmacological activation of REV-ERBs is lethal in cancer and oncogene-induced senescence,” Nature, vol. 553, no. 7688, pp. 351–355, 2018.
View at: Publisher Site | Google Scholar
M. Baraldo, “The influence of circadian rhythms on the kinetics of drugs in humans,” Expert Opinion on Drug Metabolism & Toxicology, vol. 4, no. 2, pp. 175–192, 2008.
View at: Publisher Site | Google Scholar
I. C. Mason, J. Qian, G. K. Adler, and F. Scheer, “Impact of circadian disruption on glucose metabolism: implications for type 2 diabetes,” Diabetologia, vol. 63, no. 3, pp. 462–472, 2020.
View at: Publisher Site | Google Scholar
D. Montaigne, X. Marechal, T. Modine et al., “Daytime variation of perioperative myocardial injury in cardiac surgery and its prevention by Rev-Erbα antagonism: a single-centre propensity-matched cohort study and a randomised study,” The Lancet, vol. 391, no. 10115, pp. 59–69, 2018.
View at: Publisher Site | Google Scholar
L. Fu and C. C. Lee, “The circadian clock: pacemaker and tumour suppressor,” Nature Reviews. Cancer, vol. 3, no. 5, pp. 350–361, 2003.
View at: Publisher Site | Google Scholar
C. Yang, J. Wu, X. Liu et al., “Circadian rhythm is disrupted by ZNF704 in breast carcinogenesis,” Cancer Research, vol. 80, no. 19, pp. 4114–4128, 2020.
View at: Publisher Site | Google Scholar
T. Papagiannakopoulos, M. R. Bauer, S. M. Davidson et al., “Circadian rhythm disruption promotes lung tumorigenesis,” Cell Metabolism, vol. 24, no. 2, pp. 324–331, 2016.
View at: Publisher Site | Google Scholar
M. J. De La Cruz, D. D. Shuboni-Mulligan, N. Briceno et al., “Association of circadian clock gene expression with glioma tumor microenvironment and patient survival,” Cancers, vol. 13, no. 11, 2021.
View at: Publisher Site | Google Scholar
R. V. Puram, M. S. Kowalczyk, C. G. de Boer et al., “Core circadian clock genes regulate leukemia stem cells in AML,” Cell, vol. 165, no. 2, pp. 303–316, 2016.
View at: Publisher Site | Google Scholar
I. Aiello, M. M. Fedele, F. Román et al., “Circadian disruption promotes tumor-immune microenvironment remodeling favoring tumor cell proliferation,” Science Advances, vol. 6, no. 42, 2020.
View at: Publisher Site | Google Scholar
P. Chen, W. H. Hsu, A. Chang et al., “Circadian regulator CLOCK recruits immune-suppressive microglia into the GBM tumor microenvironment,” Cancer Discovery, vol. 10, no. 3, pp. 371–381, 2020.
View at: Publisher Site | Google Scholar
G. Wichmann, M. Rosolowski, K. Krohn et al., “The role of HPV RNA transcription, immune response-related gene expression and disruptive TP53 mutations in diagnostic and prognostic profiling of head and neck cancer,” International Journal of Cancer, vol. 137, no. 12, pp. 2846–2857, 2015.
View at: Publisher Site | Google Scholar
M. E. Ritchie, B. Phipson, D. Wu et al., “limma powers differential expression analyses for RNA-sequencing and microarray studies,” Nucleic Acids Research, vol. 43, no. 7, article e47, 2015.
View at: Publisher Site | Google Scholar
G. Yu, L.-G. Wang, Y. Han, and Q.-Y. He, “clusterProfiler: an R package for comparing biological themes among gene clusters,” Omics: a Journal of Integrative Biology, vol. 16, no. 5, pp. 284–287, 2012.
View at: Publisher Site | Google Scholar
M. Kanehisa, M. Furumichi, M. Tanabe, Y. Sato, and K. Morishima, “KEGG: new perspectives on genomes, pathways, diseases and drugs,” Nucleic Acids Research, vol. 45, no. D1, pp. D353–D361, 2017.
View at: Publisher Site | Google Scholar
A. Subramanian, P. Tamayo, V. K. Mootha et al., “Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles,” Proceedings of the National Academy of Sciences of the United States of America, vol. 102, no. 43, pp. 15545–15550, 2005.
View at: Publisher Site | Google Scholar
S. Hänzelmann, R. Castelo, and J. Guinney, “GSVA: gene set variation analysis for microarray and RNA-seq data,” BMC Bioinformatics, vol. 14, no. 1, pp. 1–15, 2013.
View at: Publisher Site | Google Scholar
S. Mariathasan, S. J. Turley, D. Nickles et al., “TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells,” Nature, vol. 554, no. 7693, pp. 544–548, 2018.
View at: Publisher Site | Google Scholar
K. Yoshihara, M. Shahmoradgoli, E. Martínez et al., “Inferring tumour purity and stromal and immune cell admixture from expression data,” Nature Communications, vol. 4, no. 1, p. 2612, 2013.
View at: Publisher Site | Google Scholar
A. M. Newman, C. L. Liu, M. R. Green et al., “Robust enumeration of cell subsets from tissue expression profiles,” Nature Methods, vol. 12, no. 5, pp. 453–457, 2015.
View at: Publisher Site | Google Scholar
P. Charoentong, F. Finotello, M. Angelova et al., “Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade,” Cell Reports, vol. 18, no. 1, pp. 248–262, 2017.
View at: Publisher Site | Google Scholar
L. Barrueto, F. Caminero, L. Cash, C. Makris, P. Lamichhane, and R. R. Deshmukh, “Resistance to checkpoint inhibition in cancer immunotherapy,” Translational Oncology, vol. 13, no. 3, article 100738, 2020.
View at: Publisher Site | Google Scholar
S. J. Turley, V. Cremasco, and J. L. Astarita, “Immunological hallmarks of stromal cells in the tumour microenvironment,” Nature Reviews Immunology, vol. 15, no. 11, pp. 669–682, 2015.
View at: Publisher Site | Google Scholar
S. Masri and P. Sassone-Corsi, “The emerging link between cancer, metabolism, and circadian rhythms,” Nature Medicine, vol. 24, no. 12, pp. 1795–1803, 2018.
View at: Publisher Site | Google Scholar
A. A. Shafi and K. E. Knudsen, “Cancer and the circadian clock,” Cancer Research, vol. 79, no. 15, pp. 3806–3814, 2019.
View at: Publisher Site | Google Scholar
R. L. Ferris, “Immunology and immunotherapy of head and neck cancer,” Journal of Clinical Oncology, vol. 33, no. 29, pp. 3293–3304, 2015.
View at: Publisher Site | Google Scholar
G. I. Evan and K. H. Vousden, “Proliferation, cell cycle and apoptosis in cancer,” Nature, vol. 411, no. 6835, pp. 342–348, 2001.
View at: Publisher Site | Google Scholar
T. Yaguchi, Y. Goto, K. Kido et al., “Immune suppression and resistance mediated by constitutive activation of Wnt/β-catenin signaling in human melanoma cells,” Journal of Immunology, vol. 189, no. 5, pp. 2110–2117, 2012.
View at: Publisher Site | Google Scholar
M. J. Parsons, T. Tammela, and L. E. Dow, “WNT as a driver and dependency in cancer,” Cancer Discovery, vol. 11, no. 10, pp. 2413–2429, 2021.
View at: Publisher Site | Google Scholar
N. McGranahan and C. Swanton, “Clonal heterogeneity and tumor evolution: past, present, and the future,” Cell, vol. 168, no. 4, pp. 613–628, 2017.
View at: Publisher Site | Google Scholar
M. Shang, H. Yang, R. Yang et al., “The folate cycle enzyme MTHFD2 induces cancer immune evasion through PD-L1 up-regulation,” Nature Communications, vol. 12, no. 1, p. 1940, 2021.
View at: Publisher Site | Google Scholar
X. Zhang, R. Schulz, S. Edmunds et al., “MicroRNA-101 suppresses tumor cell proliferation by acting as an endogenous proteasome inhibitor via targeting the proteasome assembly factor POMP,” Molecular Cell, vol. 59, no. 2, pp. 243–257, 2015.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2023 Ling Aye et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

362

Downloads

385

Citations