The circadian clock is associated with prognosis and immune infiltration in stomach adenocarcinoma

Data sets

We downloaded gene-expression data of STAD from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) on January 12, 2021, and the TCGA STAD dataset (n = 415) was isolated for analysis. Moreover, we also downloaded the clinical information of each STAD patient on January 12, 2021, such as sex, tumor grade and survival status. Gene expression between STAD and gastric tissue was analyzed with Student's t-test. A sanguini diagram was built based on the R software package ggalluvial. R software v4.0.3 was utilized to perform Statistical analyses with a p-value of 0.05 as the threshold value.

Genetic mutation and pathway activity analysis

The “maftools” package was utilized to analyze and visualize genetic mutation data obtained from the TCGA database. Number of mutated sample/all STAD Sample were defined as the gene SNV percentage. The horizontal histogram showed the genes with higher mutation frequency in STAD patients. Pathway activity module of GSCALite, a TCGA visualization web portal, was used to detect the role in the pathway activity analysis. And the method was described before [10].

Drug sensitivity analysis

To analyze the correlation between the gene expression of core circadian clock genes and drug sensitivity, we collected small molecules or drugs from the Therapeutics Response Portal (CTRP). We downloaded the area under the dose-response curve (AUC) values for drugs and gene expression profiles for all cancer cell lines. The correlation between gene expression with drug sensitivity (IC50) was explored with Spearman correlation analysis.

Functional enrichment analysis

The differential expression of mRNAs with |Log (Fold Change)| >3 thresholds was firstly identified by using “limma” package. After that, the "ggplot2" package in R software was used to perform Gene Ontology (GO) analysis, including the biological process (BP), cellular component (CC), and molecular function (MF) categories. Moreover, this package was also utilized to perform Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis.

Survival analysis

The median gene expression was set as cut-off value to divide STAD cohort into two groups (high vs low). The survival difference was detected with a Kaplan-Meier analysis using TCGA STAD and GSE62254 dataset. Log-rank test was conducted to calculate the p-values and hazard ratios (HRs) with 95% confidence intervals (CIs). All analytical methods above were conducted with R software. Moreover, univariate and multivariate Cox regression analyses were conducted to find out the appropriate factors to construct a nomogram, which could be utilized to compute the prognostic risk for STAD patients.

Validation of clinical tissues

After approved by the Ethics Committees of Second Affiliated Hospital of Nanchang University, STAD tissues (n = 46) and the corresponding gastric tissues (located 2 cm away from STAD lesions) were obtained from patients and each patient had received and told the informed consent. None of patients receive treatment before operation. Total RNA of tissues and qRT-PCR were conducted as described of our study before [11]. Supplementary Table 1 presented the primer sequences. The difference of PER1/NR1D1 expression and the prognosis of PER1/NR1D1 in STAD were evaluated with Student's t-test and Kaplan-Meier analysis.

Immune infiltration

The correlation between the gene expression of core circadian clock genes and the abundances of six immune infiltrates (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells) was analyzed with TIMER (https://cistrome.shinyapps.io/timer/), a web portal for immune infiltration. Moreover, a two-sided Wilcoxon rank-sum test was performed to detect the difference of tumor infiltration levels affecting by different somatic copy number alterations.

Defining core circadian clock genes in STAD

Gene expression profile was used to explore the mRNA level of core circadian genes in gastric tissues at different times, and it was found that these genes showed oscillation at different times (Figure 1A). Using the TCGA STAD dataset, we then computed the mRNA level of these genes in STAD tissues and gastric tissues. The result is shown in Figure 1B. The mRNA levels of CLOCK (p = 0.0173), CRY1 (p = 3.99E-5) and NR1D1 (p = 3.97E-9) were upregulated, while the mRNA levels of PER1 (p = 0.0007), RORA (p = 0.0173), CRY2 (p = 0.0045) and PER3 (p = 0.0044) and were decreased in STAD tissues versus gastric tissues (Figure 1B). This evidence demonstrates that core circadian clock genes are widely altered at the mRNA level in STAD.

Figure 1. The circadian rhythm and mRNA level of core circadian clock genes. (A) The circadian rhythm of core circadian genes, including CLOCK(A), ARNTL(B), CRY1(C), CRY2(D), NR1D1(E), PER2(F) and PER3(G). (B) Box plots derived from TCGA STAD dataset comparing the expression of CLOCK, ARNTL, CRY1, CRY2, NR1D1, PER1, PER2, PER3 and RORA in STAD tissue and normal tissues.

Gene mutation landscape, hallmark pathways and drug sensitivity analysis

The gene mutation landscape of core circadian clock genes is presented in Figure 2. The results showed that genetic mutations of core circadian clock genes included inframe_mutation, Missense_nutation, splice_mutation, truncating mutation, amplification, deep deletion, mRNA high and mRNA low (Figure 2A). Among these genes, NR1D1 was the most mutated core circadian clock gene, with a genetic mutation rate of 11% in the TCGA STAD cohort (Figure 2A). To elaborate underlying functions of disorder of the circadian clock, our study then analyzed association between core circadian clock genes and well-known cancer-related hallmark pathways in STAD. We found that core circadian genes exerted inhibitory effects on the apoptosis pathway, cell cycle pathway, and DNA damage response pathway (Figure 2A). Moreover, certain circadian genes exerted promoting effects on the RAS/MAPK pathway and the RTK pathway (Figure 2B). Correlation analysis suggested a positive correlation between each member of the core circadian genes (Figure 2C). To verify whether the core circadian genes could serve as therapeutic targets, it is important to evaluate their association with existing drug targets. Clock-related drug administration at specific times during circadian cycles may remarkably improve the effectiveness of the drug and reduce the toxicity. In our study, drug sensitivity analysis revealed that the expression of PER2, CRY2 and PER1 was linked to drug sensitivity (negative correlation), while the expression of NR1D1, CLOCK, PER3 and ARNTL was linked to drug resistance (positive correlation) (Figure 2D).

Figure 2. Genetic mutation landscape and drug sensitivity analysis of core circadian clock genes in STAD. (A) Oncoplot displaying genetic mutation landscape of core circadian clock genes in TCGA STAD cohort. (B) the role of circadian clock genes in the famous cancer related pathways. (C) A heat map of the correlation between each member of core circadian clock genes. (D) the correlation between core circadian clock genes and drug or small molecules. The Spearman correlation represent the core circadian clock genes expression correlates with the drug. The positive correlation means that the gene high expression is resistant to the drug, vise verse. ^*p < 0.05, ^**p < 0.01.

Enrichment analysis of differentially expressed genes (DEGs) in STAD

We then performed GO and KEGG analyses using DEGs and the heatmap of DEGs in STAD is shown in Figure 3A. GO analysis demonstrated the correlation of these DEGs in organelle fission, nuclear division, mitotic nuclear division, condensed chromosome, glycosaminoglycan binding, rhythmic process, and chemokine activity (Figure 3B). Moreover, KEGG pathway demonstrated that these DEGs were primarily enriched in the cell cycle, cytokine-cytokine receptor interaction, circadian rhythm, p53 signaling pathway, STAD and DNA replication (Figure 3C).

Figure 3. The enrichment analysis of different expression genes in STAD. (A) Heatmap of different expression genes in STAD. (B)The enriched items in GO analysis. (C) The enriched items in KEGG analysis.

Prognosis value of core circadian genes in STAD

The risk score of core circadian genes in the TCGA STAD cohort was detected with Cox regression analysis. Kaplan-Meier curves were drawn on the basis of the risk score of each patient. The data of prognostic analysis of core circadian genes in STAD was shown in Table 1. Among these core circadian genes, the expression of PER1 and NR1D1 was linked to overall survival (OS), progression-free survival (PFS), and disease-free survival (DFS) rates in STAD patients. For the overall survival rate, the risk score of each patient is shown in Figure 4A and Figure 5A. STAD patients with high expression of NR1D1 (Figure 4B, p = 0.002, HR (95% CI) = 1.69 (1.21–2.36)) and PER1 (Figure 5B, p = 0.007, HR (95% CI) = 1.59 (1.14–2.22)) had poor OS rates with 5-year AUCs of 0.611 (Figure 4C) and 0.58 (Figure 5C), respectively. For the progression-free survival rate, the risk score of each patient is shown in Supplementary Figure 1A and Supplementary Figure 2A. STAD patients with high expression of NR1D1 (Supplementary Figure 1B, p = 0.008, HR (95% CI) = 1.62 (1.13–2.31)) and PER1 (Supplementary Figure 2B, p = 0.007, HR (95% CI) = 1.64 (1.14–2.35)) had poor PFS rates with 5-year AUCs of 0.615 (Supplementary Figure 1C) and 0.562 (Supplementary Figure 2C), respectively. For the disease-free survival rate analysis, the risk score of each patient is shown in Supplementary Figure 3A and Supplementary Figure 4A. STAD patients with high expression of NR1D1 (Supplementary Figure 3B, p = 0.002, HR (95% CI) = 1.92 (1.27–3.02)) and PER1 (Supplementary Figure 4B, p = 0.010, HR (95% CI) = 1.69 (1.10–2–60)) had poor DFS rates with 5-year AUCs of 0.642 (Supplementary Figure 3C) and 0.622 (Supplementary Figure 4C), respectively. To further verify the prognostic value of NR1D1 and PER2 in STAD, the GSE62254 dataset was obtained from the GEO database. The results indicated poor OS and PFS rates in STAD patients with high NR1D1 (Figure 6A, all p < 0.05) and PER2 (Figure 6B, all p < 0.05) expression. These results indicated that NR1D1 and PER2 serve as potential prognostic biomarkers in STAD.

Figure 4. The overall survival analysis of NR1D1 in STAD. (A) The risk score, survival status and gene expression of each patient. (B) Kaplan-Meier overall survival curve of NR1D1 in STAD patients with high and low NR1D1 expression. (C) Time-dependent ROC of NR1D1 in predicting the prognosis of STAD patients. ROC receiver operating characteristic.

Figure 5. The overall survival analysis of PER1 in STAD. (A) The risk score, survival status and gene expression of each patient. (B) Kaplan-Meier overall survival curve of PER1 in STAD patients with high and low PER1 expression. (C) Time-dependent ROC of PER1 in predicting the prognosis of STAD patients. ROC receiver operating characteristic.

Figure 6. The prognosis value and expression of core circadian clock genes in the subtypes of STAD tissues. (A, B) The overall survival and progression-free survival curve of STAD patients with high/low NR1D1(A) and PER1(B) expression. (C, D) The relative expression of NR1D1 and PER1 in STAD tissues and normal tissues in clinical cohort. (E, F) The overall survival curve of STAD patients with high/low NR1D1 and PER1 expression in clinical cohort. (G) The expression of NR1D1 in the subtypes of STAD tissues. (H) The expression of PER1 in the subtypes of STAD tissues. Each column represents a characteristic variable, different colors represent different subtypes of STAD tissues, and lines represent the distribution of the same sample in different characteristic variables.

Validation of the expression and prognostic value

Clinical tissues were used to further detect the mRNA level and prognostic significance of PER1/NR1D1 in STAD. Not surprisingly, NR1D1 expression was increased in STAD (Figure 6C, p = 0.00001) while PER1 expression was decreased in STAD tissues (Figure 6D, p = 0.0069) compared with that in normal tissues. Moreover, As shown in Figure 6E and Figure 6F, prognosis analysis suggested a poor overall survival in STAD patients with high NR1D1 expression (p = 0.012) and PER1 expression (p = 0.0064), These data further verified above results.

The expression of NR1D1 and PER1 in STAD patient subtypes

The above results (result 3.4 and result 3.5) revealed that NR1D1 and PER1 may serve as prognostic biomarkers in STAD and are predictors of poor patient prognosis. Thus, NR1D1 and PER1 were selected for further analysis. We then analyzed the association between NR1D1/PER1 expression and clinicopathologic features. As shown in Figure 6G, STAD patients with high pTNM stage and high tumor grade had low NR1D1 expression compared with that of patients with low pTNM and low tumor grade. The PER1 results are presented in Figure 6H. Male STAD patients and patients with high tumor-grade had low PER1 expression.

NR1D1/PER1 were correlated with immune infiltration in STAD

Immune infiltration is a vital factor affecting sentinel lymph node status and prognosis in cancers [12–14]. Recent evidences have elaborated the close correlation between circadian rhythm and lung cancer [5]. To identify whether core circadian genes could act as immunotherapy targets for STAD, we first clarified the correlation between the mRNA level of NR1D1/PER1 and immune infiltration in STAD. We found that the mRNA level of NR1D1 was negatively associated with the abundance of CD8+ T cells (Cor = –0.164, p = 1.58E-3), macrophages (Cor = –0.11, p = 3.41E-2), neutrophils (Cor = –0.188, p = 2.75E-4), and dendritic cells (Cor = –0.172, p = 8.57E-4) (Figure 7A). PER1 expression was positively correlated with the abundance of B cells (Cor = 0.212, p = 3.39E-5), CD8+ T cells (Cor = 0.238, p = 3.69E-6), CD4+ T cells (Cor = 0.417, p = 7.79E-17), macrophages (Cor = 0.499, p = 1.15E-24), neutrophils (Cor = 0.233, p = 5.88E-6) and dendritic cells (Cor = 0269, p = 1.36E-7) (Figure 7B). Moreover, we found that somatic copy number alteration (SCNA) of NR1D1 (Figure 7C) and PER1 (Figure 7D) could partially suppress immune infiltration in STAD.

Figure 7. The correlation between core circadian clock genes and immune infiltration (TIMER). (A) The correlation between NR1D1 expression and the abundance of CD8+ T cells, CD4+ T cells, Macrophage, Neutrophils and Dendritic cells. (B) The correlation between PER1 expression and the abundance of CD8+ T cells, CD4+ T cells, Macrophage, Neutrophils and Dendritic cells. (C) the correlation between SCNA of NR1D1 and immune cell infiltration. (D) the correlation between SCNA of PER1 and immune cell infiltration. (E) the correlation between core circadian clock genes and the expression of chemokines in STAD. (F) the correlation between core circadian clock genes and the expression of chemokine receptors in STAD. Red color represents positive correlation, blue color represents negative correlation. SCNA, somatic copy number alterations; ^*P < 0.05, ^**P < 0.01, ^***P < 0.001.

Additionally, we sought to calculate the association between NR1D1/PER1 expression and immune cell markers. As shown in Table 2, NR1D1 expression was correlated with certain biomarkers of these immune cells, while PER1 expression showed a positive correlation with most of the biomarkers of these immune cells (Table 2). Chemokines and their receptors modulate immune surveillance in oncogenesis, metastasis, and response to immunotherapy [15]. In our study, the correlation analysis demonstrated a strong correlation between PER1/NR1D1 and chemokines as well as chemokine receptors (Figure 7E–7F). These data demonstrated the possible correlation between the circadian clock and immune infiltrates in STAD.

Consent for publication

Our research was approved and supported by an ethics committee of the Second Affiliated Hospital of Nanchang University. All patients provided informed consent.

Availability of data and materials

The analyzed data sets generated during the study are available from the corresponding author on reasonable request.