Introduction

In terms of prevalence, GC ranks fifth globally. GC also ranks fourth in terms of global mortalities that are related to cancer. It was projected that 1.08 million and 770,000 cases of GC new cases and deaths occurred, correspondingly, in 20201. Early-stage gastric cancer is mainly treated by surgical resection, while advanced-stage gastric cancer is primarily treated with a combination of chemotherapy and immunotherapy. Most GC patients are diagnosed with local or distant metastasis, leading to poor prognosis2. Stomach adenocarcinoma (STAD) represents the primary histological subtype of GC. The occurrence and progression of STAD are complex and heterogeneous processes influenced by genetic variations, the gastrointestinal microbiota, and external environmental factors. Due to the heterogeneity of STAD, recurrence and metastasis remain the major causes of treatment failure. Hence, identifying molecular biomarkers that are linked to not only the prognosis but also the therapeutic response of STAD patients is highly important.

The circadian rhythm is a conserved phenomenon in the evolution of eukaryotes that modulates the rhythmicity of physiological and behavioral processes in organisms3. In humans, over 50% of protein-coding genes are controlled by the circadian clock4. Previous studies have indicated that disruption of the circadian rhythm promotes cancer proliferation and metastasis through biological control of the behavior of tumor cells and the tumor microenvironment5,6. Some studies have also reported the influence of circadian rhythm genes on gastric cancer cell proliferation and metastasis7. Recently, although the circadian rhythm has emerged as a prominent subject of investigation in cancer research, the specific mechanisms of its involvement in cancer progression remain unclear. Currently, there is no available literature on whether circadian rhythm genes have a function as prognostic biomarkers for STAD patients.

Using STAD RNA-seq data and clinical information gleaned from the TCGA database, we developed a risk prediction signature for CRG. To validate the signature, we further assessed its performance using the Gene Expression Omnibus (GEO) dataset. In this study, we investigated the correlation between CRG and various factors, including the tumor microenvironment (TME), immune checkpoint inhibitors (ICIs), and drug sensitivity. These findings offer valuable technical and theoretical support for utilizing a CRG signature to evaluate STAD patients’ prognosis and inform clinical treatment decisions.

Materials and methods

Data collection and processing

Using the TCGA database, we gathered gene expression data and clinical information, encompassing age, gender, stage, and survival time, from STAD patients (https://portal.gdc.cancer.gov/). The TCGA gene expression dataset encompassed 412 and tumor samples and 36 normal control samples, correspondingly. After removing the samples from patients with missing information, 407 tumor samples were left for further analysis. Pathological characteristics and clinical of the patients included in the study is displayed in Table 1. Moreover, we acquired the GSE84437 dataset from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), which contains gene expression data and survival information from 483 STAD patients. To ensure uniformity, we normalized the gene expression data from both the TCGA and GSE84437 datasets to transcripts per million. We converted using the log2 (x + 1) method. We acquired the CRG dataset from the Molecular Signatures Database (https://www.gsea-msigdb.org/), which includes 201 circadian rhythm genes. We deployed differential expression analysis utilizing the “limma” R package to ascertain CRG that manifested remarkable differential expression between tumor and normal samples. The criteria for significance were |log2FC|≥0.585 and a false discovery rate < 0.058.

Table 1 Clinical and pathological characteristics of STAD patients in the training and testing cohort.
Full size table

Construction of the gene signature

A total of 407 TCGA-STAD patients were segregated into two distinct groups via randomization in our study. These two groups encompassed the training cohort comprising 265 patients and a testing cohort consisting of 142 patients. The 483 STAD patients from the GSE84437 dataset served as the validation set. To further screen the CRG and identify candidate genes significantly linked to prognosis, univariate Cox regression analysis was performed (p < 0.05). Next, the LASSO regression was conducted to choose the optimal set of genes associated with overall survival9. Finally, the multivariate Cox regression analysis was performed to compute the coefficients of the identified candidate genes.

Assessment of the predictive performance of the gene signature

By utilizing the median risk score, the STAD samples were classified into high- and low-risk groups. The generated K‒M survival curves for the training cohort, testing cohort, and validation cohort were instrumental in the examination of variations in patient survival rates. Subsequently, ROC curves were generated for the gene signatures to further evaluate their effectiveness and accuracy. Additionally, progression-free survival (PFS), principal component analysis (PCA), concordance index (C-index), and clinical correlation analysis were utilized ascertaining the predictive ability of the gene signature. The “survival, " “survminer, " and “timeROC” R packages were utilized in this analysis.

Establishment of a nomogram

To improve the ability of the prognosis for STAD, we established a nomogram integrating the risk score as well as the clinical characteristics. Subsequently, we validated its predictive ability through survival curves together with an area under the curve (AUC). In this analysis, the “regplot” and “rms” R packages were utilized.

Assessment of mutations

The MAF files containing mutation data of STAD patients were obtained from the TCGA10. Utilizing the “tcgaCompare” function of the “maftools” R package, we compared the mutation burden in the MAF file with the mutation burden in the TCGA MC3 project across 33 cancer types and calculated the tumor mutation burden (TMB)11. K‒M survival curves evaluated the link between the TMB and survival time.

Assessment of the tumor immune signature

Using estimation of stromal and Immune cells in malignant Tumor tissues using Expression data (ESTIMATE) algorithm, gene expression data can predict stromal and immune cell infiltration levels in tumors12. This approach yields four scores: stromal score (infiltration of stromal cells), immune score (infiltration of immune cells), ESTIMATE score (sum of the stromal and immune scores), and tumor purity13.The CIBERSORT algorithm was used to evaluate the relative abundance of 22 distinct subtypes of immune cells, and 1000 simulations were conducted to ensure the robustness of the computations14. We conducted this analysis using the “limma” R package.

Immunotherapy prediction

The tumor immune dysfunction and exclusion (TIDE) score was executed by the TIDE online platform to anticipate the probability of tumor immune evasion by analyzing the gene expression profile of tumor samples15. The immune phenotype score (IPS) for STAD patients was downloaded from The Cancer Immunome Atlas (TCIA) database to anticipate the therapeutic response to ICI treatment16.

Drug sensitivity

Genomics of Drug Sensitivity in Cancer database provided the IC50 values of the therapeutic agents utilized in chemotherapy17. Drug sensitivity prediction was performed utilizing the “oncoPredict” R package18. We compared the response to 198 drugs used for chemotherapy drugs among various risk groups and identified drugs that exhibited significant differences in sensitivity (p < 0.001).

Enrichment analysis

Gene set enrichment analysis (GSEA) were performed to identify enriched functions and pathways among distinct risk groups19. We extracted the DEGs among the different risk groups (|log2FC|>1, FDR < 0.05). Subsequently, Gene Ontology (GO) functional enrichment analysis and KEGG pathway enrichment analysis were conducted for DEGs by utilizing the “clusterProfiler” R package20.

Protein–protein interaction networks

To ascertain the protein–protein interaction (PPI) among the DEGs, the String online database was utilized to establish a PPI network21. Significant interactions were determined based on an interaction that was larger than > 0.7. Next, the Cytoscape software (version 3.9.1, https://cytoscape.org/download_old_versions.html) was utilized to visualize the network. Finally, the CytoHubba was used to calculate the protein interaction degrees and determine the hub genes in the network22.

Assessment of hub genes

To ascertain the functions of the hub genes, STAD patients were assigned into high- and low-expression groups utilizing the median expression values of the hub genes in the TCGA dataset. Further investigations involving K‒M survival and immune cell infiltration analyses were executed with the intention of comparing the high- and low-expression groups.

Real-time quantitative polymerase chain reaction (RT-qPCR)

Eight patients with gastric tumors who underwent their first endoscopic examination at the Endoscopy Center of the Affiliated Hospital of Southwest Medical University from March 2024 to April 2024 were selected for this study. Tumor tissue samples were collected along with adjacent non-tumor tissue located 2 cm from the tumor margin. Among these patients, six were confirmed to have gastric adenocarcinoma postoperatively, and the adjacent tissues showed no evidence of cancer cell infiltration. These patients were selected for further PCR experiments. The Steady Pure Quick RNA Extraction Kit (Accurate Biotechnology, Hunan, China) was utilized to isolate the total RNA. The isolated RNA was then utilized as a template for reverse transcription reactions with TB Green® Premix Ex Taq™ II (Takara, Beijing, China) were utilized to conduct the RT-qPCR. The SYBR® Green Real-time PCR Master Mix (TOYOBO, Osaka, Japan) and a SLAN 96S Real-Time PCR System (Hongshi Medical Treatment Technology, Shanghai, China). Supplementary Table 1 outlines the primer nucleotide sequences used. GAPDH served as the housekeeping gene.

Statistics

We performed statistical analysis utilizing R software (version 4.3.1, https://cran.r-project.org/src/base/R-4/). Depending on the homoscedasticity and the Gaussian distribution of the data, either t tests or Mann‒Whitney U tests were employed. Pearson correlation investigated the link between two continuous variables. A distinctive mark was applied to indicate the significance level: not-significant (ns) for p ≥ 0.05; * denotes p < 0.05; ** denotes p < 0.01; *** denotes p < 0.001.

Results

Establishment of a circadian rhythm gene signature

The flow diagram of our study is presented in Fig. 1. To identify CRG linked to STAD prognosis, we contrasted the expression levels of CRG between tumor and normal tissue samples in the TCGA dataset. We identified 109 differentially expressed CRG and generated a gene expression heatmap (Fig. 2A). During the randomization of TCGA-STAD patients, two groups were established, a training cohort (n = 265) and a testing cohort (n = 142). To further select candidate CRG that met the criteria, we executed univariate Cox regression analysis on the 109 genes that exhibited differential expression. This analysis identified 8 candidate genes that met the criteria (ADIPOQ, CARTPT, HTR7, KCND2, OPN3, PROKR2, SERPINE1, and TP53; p < 0.05; Table 2). LASSO regression analysis was performed on the 8 genes to mitigate multicollinearity (Fig. 2B-C). Next, multivariate Cox regression analysis was applied to detect and estimate the coefficients of the candidate genes involved in the final modeling process. We constructed a 4-CRG signature (CARTPT, OPN3, SERPINE1, and TP53). The risk score was finally computed-utilizing the following formula:

$${\text{risk score}} = \left( {CARTPT \times 0.{\text{1485476}}0{\text{5}}} \right) + \left( {OPN3 \times 0.{\text{32}}0{\text{794967}}} \right) + \left( {SERPINE1 \times 0.{\text{19297654}}} \right) - \left( {TP53 \times 0.{\text{234329235}}} \right).$$
Fig. 1
figure 1

The flow diagram of this study.

Full size image
Fig. 2
figure 2

Establishment and evaluation of the circadian rhythm gene signature. (A) Heatmap of differentially expressed CRGs. Blue indicates low expression and red indicates high expression. (B) Adjusted parameter selection in the Lasso Cox model, each curve represents a gene, and the proximity of the y-axis values indicates the corresponding gene coefficients (C) Cross-validation adjusted parameter selection for LASSO regression. (D–F) K-M survival curves of the training cohort, testing cohort and validation cohort. p < 0.05 indicates that the difference in survival times between the two groups is statistically significant. (G–I) The ROC curves of the training cohort, testing cohort and validation cohort. The larger the AUC, the better the specificity and sensitivity of the signature.

Full size image
Table 2 Candidate genes of the CRG gene signature.
Full size table

Evaluation of the predicting ability of the gene signature

In the training as well as the testing cohorts, K‒M survival analysis affirmed that the high-risk group manifested a lower survival rate when contrasted with the low-risk group (p < 0.001; Fig. 2D-E), and consistent findings were affirmed in the validation cohort (p = 0.006; Fig. 2F). The AUC values for anticipating 1-, 3-, and 5-year OS were as follows: training cohort (0.668, 0.654, 0.733), testing cohort (0.683, 0.665, 0.659), and validation cohort (0.652, 0.621, 0.596) (Fig. 2G-I). These results suggested that CRG signature could predict prognosis of patient accurately.

The PCA results demonstrated that the CRG involved in the signature could significantly distinguish between high- and low-risk STAD patients, while all genes and all CRG could not (Fig. 3A-C). PFS refers to the time interval from the start of tumor observation to disease progression or patient death. The PFS K‒M curve results affirmed that the high-risk group had a worse survival time when a comparison is executed with the low-risk group (p < 0.001; Fig. 3D). Both univariate and multivariate Cox regression analyses suggested that the risk score could serve as an independent risk factor, with a hazard ratio higher than that of patient characteristics (Fig. 3E-F). Moreover, the C-index and AUC values pertinent to the risk score were greater than the other variables (Fig. 3G-H). These results indicate that CRG gene signature possess good prognostic predictive ability and are superior to traditional characteristics.

Fig. 3
figure 3

Assessment of predictive performance of the CRG signature. (A-C) Principal component analysis of all genes, CRGs and model CRGs. The dot in different colors represent patients with varying risk levels. (D) The K-M survival curve of progression free survival. (E-F) Univariate regression and multivariate regression analyses. Hazard Rate (HR) > 1 denotes a risk factor, HR < 1 denotes a protective factor. (G) The concordance index (CI) of risk score and clinical characteristics. The higher the value of the CI, the better the predictive ability of the model. (H) The ROC curve of the risk score and clinical characteristics.

Full size image

Clinical evaluation of the gene signature

We evaluated the risk score in distinct subgroups utilizing the clinical characteristics. Pearson correlation analysis showed that the risk score was significantly correlated only with the T stage, while it had no significant correlations with the rest characteristics (Fig. 4A). To estimate the reliability of the CRG signature in different clinical characteristic subgroups, the STAD samples were categorized into 4 subgroups. These subgroups encompassed age subgroups (<= 65 years and > = 65 years), grade subgroups (G1-2 and G3-4), gender subgroups (female and male), and stage subgroups (I-II and III-IV). K‒M survival analysis affirmed that in the age, gender, and stage subgroups, the high-risk group was affirmed to exhibit a lower survival rate than the low-risk group (p < 0.001, Fig. 4B-G). There is a similar result in both the G1-2 subgroup (p = 0.01, Fig. 4H) and the G3-4 subgroup (p < 0.001, Fig. 4I). These results showed that the CRG signature could predict the survival probabilities of patients across different clinical subgroups.

Fig. 4
figure 4

Clinical correlation and survival analysis of patients stratified by the CRG gene signature. (A) Risk score correlation with clinical characteristics. The colors represent different clinical characteristics. p-value < 0.05 indicates a significant correlation. Survival analysis according to (B–C) age, (D–E) gender, (F–G) tumor stage, and (H–I) tumor grade.

Full size image

Establishment and evaluation of a nomogram

To provide a more accurate prediction of the patient’s risk probability, the risk score was integrated with clinical characteristics to develop a prognostic nomogram (Fig. 5A). By calculating the scores of the variables and summing them, the total score could be used to predict patient’s survival time. According to the univariate Cox regression analysis, gender, grade, stage, and the nomogram were linked to prognosis, while age was not (Fig. 5B). Multivariate Cox regression implied that only the nomogram could anticipate the survival in of STAD patients as an independent risk factor (Fig. 5C). Subsequently, calibration curve of nomogram was plotted. The calibration curves for 1-year, 3-year, and 5-year were all close to the ideal curve (Fig. 5D), indicating that the nomogram has strong predictive capability. Next, we plotted ROC curves for clinical characteristics, the risk score, and the nomogram. The AUC value for the established nomogram was 0.758, and this value denoted the highest among those of the other variables (Fig. 5E). This result indicated that the nomogram could be an instrumental tool in anticipating the prognosis of STAD patients.

Fig. 5
figure 5

Establishment and evaluation of a nomogram. (A) The nomogram of the CRG signature. The total point is the sum of the scores of all variables. (B, C) Univariate regression and multivariate regression analyses of the nomogram. HR > 1 denotes a risk factor, HR < 1 denotes a protective factor. (D) The calibration curve of the nomogram. The gray straight line represents the ideal curve; the closer the model curve is to the ideal curve, the better the calibration of the model. (E) The ROC curve of the nomogram and clinical characteristics.

Full size image

Analysis of mutations

The onset and the progression of tumors can be attributed to somatic gene mutations23. Therefore, we assessed the relationship between the CRG signature and mutations in the TCGA-STAD cohort. Within the high-risk group, the following indicated the 5 genes that exhibit remarkable mutation frequencies: TTN (45%), TP53 (45%), MUC16 (26%), LRP1B (23%), and CSMD3 (21%), with missense mutations standing out as prevalent mutation type (Fig. 6A). In contrast, in the low-risk group, the most frequent mutations were in TTN (58%), TP53 (40%), MUC16 (36%), ARID1A (36%), and LRP1B (31%), with ARID1A having frameshift mutations as the dominant mutation type, while the others had missense mutations as the dominant mutation type (Fig. 6B). These results indicated that there are no significant differences in the mutated genes, mutation types, and mutation frequencies among different risk groups. Next, we compared the risk scores of the most frequently mutated genes (TTN, TP53, and MUC16) between the mutant-type and wild-type samples to determine whether there were any differences. The results affirmed that the risk score was greater for the TP53 mutant-type samples than for the wild-type samples (p = 0.012), while no notable changes were detected in the rest of the genes (Fig. 6C). Additionally, TMB analysis implied that patients in the high-risk group had a diminished TMB than those in the low-risk group (p < 0.001; Fig. 6D). Finally, K‒M curves revealed that high TMB group exhibited a greater survival rate in comparison to the low TMB group (p = 0.006; Fig. 6E). This results indicated that TMB is associated with the risk score, and lower mutational burden is generally linked to poorer prognosis, which is consistent with previous findings regarding “cold tumors” in immunotherapy24.

Fig. 6
figure 6

Analysis of mutations. (A-B) The gene mutation frequency and type in the high- and low-risk group. The colors represent different types of mutations, and the vertical lines indicate the patients with mutations. (C) The risk scores of patients with wild-type and mutant TP53. (D) The tumor mutation burden. (E) K-M survival curve according to the TMB.

Full size image

Assessment of the tumor microenvironment (TME)

The TME is the local environment whereby tumor cells exist and exert a substantial influence on tumor growth and metastasis25. Since the CRG signature could anticipate the survival probability of STAD patients, we investigated its correlation with the TME. The ESTIMATE algorithm was employed to assess TME in STAD samples. The findings revealed that compared to the high-risk group, the stromal score and ESTIMATE score were lower in the low-risk group (p < 0.001; Fig. 7A-B), while the tumor purity score was elevated (p < 0.001; Fig. 7C), and no notable changes in immune score between the two groups (Fig. 7D).Increase in stromal scores, along with a decrease in tumor purity, fostered the development of an immunosuppressive microenvironment, which subsequently contributes to resistance to immunotherapy26.To further ascertain the correlation between the CRG signature and the TME in STAD, the CIBERSORT algorithm estimated the relative proportions of immune cells in the STAD samples. The analysis of immune cells revealed that in comparison to the high-risk group, the low-risk group had greater proportions of memory B cells, CD8 + T cells, and regulatory T cells, but fewer neutrophils, resting mast cells, eosinophils, and M2 macrophages (p < 0.05; Fig. 7E). The analysis of immune checkpoints revealed that in comparison to the low-risk group, the high-risk group had greater expression of NRP1, CD276, CD86, CD44, CD200, PDCD1LG2, and TNFSF4, but lower expression of TNFRSF18, LGALS9, and TNFRSF14 (p < 0.05; Fig. 7F). These results indicated that the high-risk group is more prone to form an immunosuppressive microenvironment, leading to poorer prognosis.

Fig. 7
figure 7

Analysis of the tumor microenvironment. (A-D) Stromal score, ESTIMATE score, tumor purity score, and immune score of STAD in the high- and low-risk groups. (E) The relative abundance of 22 immune cell infiltrates between the high- and low-risk groups. (F) Immune checkpoint gene expression between the high- and low-risk groups.

Full size image

Response to ICIs and drug sensitivity

Immunotherapy and chemotherapy represent the primary therapeutic approaches for advanced gastric cancer. Both the TIDE and IPS have demonstrated their ability to anticipate the reaction of patients with cancer to ICIs27. Hence, we explored the potential link between the CRG signature and TIDE and IPS in STAD. TIDE analysis revealed that, compared to the low-risk group, the high-risk group exhibited higher TIDE scores (p < 0.05; Fig. 8A). IPS analysis revealed that among the four types—ctla4_neg_pd1_neg, ctla4_neg_pd1_pos, ctla4_pos_pd1_neg, and ctla4_pos_pd1_pos—the low-risk group displayed higher IPS scores than did the high-risk group (p < 0.001; Fig. 8B-E). These results indicated that the CRG signature can effectively predict the risk of immune evasion and the response to immunotherapy in STAD individuals.

Furthermore, we evaluated the performance of the CRG signature for predicting sensitivity to chemotherapy drugs. In the 198 anticancer drugs screened, we identified 63 drugs with differential sensitivity, of which 60 drugs manifested remarkable sensitivity in the high-risk group than in the low-risk group; these drugs included frequently utilized chemotherapeutic drugs for gastric cancer such as 5-fluorouracil, oxaliplatin, sorafenib, and afatinib (p < 0.001; Fig. 8F-I). These results indicate that the CRG signature can provide guidance for the selection of chemotherapy drugs for STAD individuals.

Fig. 8
figure 8

Response to ICIs and drug sensitivity. (A) Differences in TIDE scores between the high- and low-risk groups. The higher the TIDE value, the greater the likelihood of immune evasion. (B–E) Differences in the IPS between the high-risk and low-risk groups Ips_ctla4_neg_pd1_neg: Indicates no response to both anti-CTLA-4 and anti-PD-1 antibodies; ips_ctla4_neg_pd1_pos: Indicates a lack of response to anti-CTLA-4 with a positive response to PD-1.ips_ctla4_pos_pd1_neg: Indicates a response to anti-CTLA-4 with no response to PD-1. ips_ctla4_pos_pd1_pos: Indicates a positive response to both anti-CTLA-4 and anti-PD-1 antibodies. Higher IPS are associated with higher response rates. (G-I) Differences in the sensitivity to 5-fluorouracil, oxaliplatin, sorafenib, and afatinib between the high- and low-risk groups.

Full size image

Enrichment analysis and the development of protein–protein interaction networks

As per the GSEA enrichment analysis, the high-risk group manifested enrichment of pathways encompassing calcium signaling pathway, and focal adhesion, while the low-risk group exhibited enrichment of pathways including allograft rejection, oxidative phosphorylation, primary immunodeficiency, and ribosome (Fig. 9A-B). Among these pathways, the calcium signaling pathway and oxidative phosphorylation were closely linked to cancer progression.

We further extracted the DEGs between the high- and low-risk groups whereby we identified 526 genes. Then, GO and KEGG analyses were performed. GO analysis indicated that the DEGs were dominantly enriched in muscle contraction, neuronal cell body, and extracellular matrix structural constituent, none of which showed a significant association with cancer cells (Fig. 9C). KEGG enrichment analysis revealed pathways like the PI3K-Akt signaling pathway, cAMP signaling pathway, and calcium signaling pathway, all of which are closely associated with cancer cell proliferation, apoptosis, and metastasis (Fig. 9D).

To study the protein-protein interactions encoded by the DEGs, we constructed PPI networks. First, the String online database was utilized establish the PPI network. Then, we visualized the network using Cytoscape (Fig. 9E). Finally, using CytoHubba, 10 hub genes, namely, FN1, IL6, IGF1, FGF2, FLNC, NPY, NGF, FLNC, SST, and GCG, were identified in the network (Fig. 9F). These hub genes among DEGs may influence the prognosis of STAD.

Fig. 9
figure 9

Enrichment analysis and protein–protein interaction networks. (A-B) GSEA enriched pathways in the high- and low-risk group. The X-axis represents genes sorted in order, while the Y-axis indicates the corresponding Running Enrichment Score. The peaks in the line graph correspond to the enrichment scores for this gene set. Different colored curves represent various pathways, and the gray area plot displays the signal noise of the genes. (C-D) GO and KEGG enrichment analyses of DEGs between the between high- and low-risk groups. The size of the bubbles or the length of the bars represents the quantity of enriched genes. (E) Protein–protein interaction networks of DEGs between the high- and low-risk groups. The ellipses represent genes, while the lines indicate the presence of interaction relationships. (F) Hub genes of protein–protein interaction networks. The rectangles represent genes, while the lines indicate the presence of interaction relationships. Blue denotes ordinary genes, while colors represent hub genes.

Full size image

Estimation of hub genes

Hub genes are highly important in gene co-expression networks and can serve as potential therapeutic targets or biomarkers28. FN1 has a high degree of connectivity in the gene expression network. Previous literature has propounded that FN1 is a soluble factor secreted by cancer-associated fibroblasts, that can enhance tumor cell proliferation and metastasis29. To assess the significance of FN1 in STAD, we performed K-M analysis, and immune cell correlation analysis. First, we analyzed FN1 expression levels in both normal and tumor samples. The findings implied a substantial upregulation of FN1 expression in tumor samples (p < 0.001; Fig. 10A). Subsequently, with the median expression value of FN1 serving as a cutoff, we segregated the TCGA-STAD samples into high- and the low-expression groups. The survival analysis findings affirmed that the group with low expression levels exhibited a greater survival rate than the other expression group, indicating that FN1 may act as a pro-cancer factor (p < 0.001; Fig. 10B). Clinical correlation analysis manifested that the FN1 expression levels were greater in the T2 stage, T3 stage, and T4 stage subgroups than in the T1 stage subgroup (p < 0.001; Fig.  10C), while there were no substantial variations in age, gender, or other stage subgroups (Supplemental Fig. 1A-D). These results indicated that the expression of FN1 does not correlate with clinical characteristics. Immune cell infiltration analyses denoted that the low-expression group affirmed markedly greater proportions of plasma cells and regulatory T cells (Tregs) than did the high-expression group, with a significant increase in M2 macrophages (p < 0.001; Fig. 10D). Pearson correlation analysis of immune cells indicated a negative interrelationship between the expression of FN1 and plasma cells as well as activated dendritic cells, revealing a positive correlation with M2 macrophages (Fig. 10E-G). The immunological analysis results suggested that FN1 may be negatively correlated with anti-tumor immunity.

Fig. 10
figure 10

Estimation of Hub genes. (A) FN1 gene expression between normal and tumor samples. (B) Differences in survival probability between the FN1 high- and low-expression groups. (C) Differences of FN1 gene expression among T stages. (D) The relative abundance of 22 immune cell infiltrates between the FN1 high- and low-expression groups. (E-G) Pearson correlation analysis of FN1 expression and immune cells. The black dots represent patients, the yellow area indicates the number of patients with a certain level of FN1 expression, and the blue area represents the number of patients with a specific relative abundance of immune cells.

Full size image

Expression levels of these signature genes and hub gene

To examine the discrepancy in the expression levels of CRG participating in the signature and FN1 between normal and tumor tissues, we gathered 6 pairs tumor and adjacent tissues from STAD patients for RT‒qPCR analysis. The results revealed elevated expression of the OPN3, and FN1 genes in tumor tissues in comparison to adjacent tissues, while TP53 expression was reduced (p < 0.05, Fig. 11A-C). CARTPT and SERPINE1 expression were furthermore heightened in tumor tissues than adjacent tissues, even though the variation was not found to manifest a statistical significance (Fig. 11D-E).

Fig. 11
figure 11

Relative expression of 4-CRGs and FN1 in 6 pairs adjacent and cancer tissue samples determined by RT‒qPCR. (A–E) OPN3, FN1, TP53, CARTPT, and SERPINE1.

Full size image

Discussion

Disruption of the normal circadian rhythm is linked to dysregulation of circadian rhythm genes and various chronic diseases, including cancer30. Recent pieces of evidence propound a strong link between disturbances in the circadian cycle and the progression of STAD31,32. However, whether circadian rhythm genes can function as prognostic biomarkers for STAD patients is unclear. Here, we used the TCGA-STAD dataset to develop a multigene prognostic signature (CARTPT, OPN3, SERPINE1, TP53) associated with gastric cancer. Univariate and multivariate Cox regression analyses unequivocally manifested that the 4-CRG signature served as an independent prognostic factor for patients with STAD. K‒M and ROC curve analyses validated the model’s ability to predict OS and PFS among STAD patients. Furthermore, we confirmed the robust predictive capability of this model in the independent GSE84437 cohort. Additionally, we formulated a nomogram that integrates risk score with clinical characteristics to increase the prediction accuracy of the prognosis of STAD patients.

Combining our results with those of existing studies, the 4-CRG signature could function as a possible biomarker for the development and recurrence risk of STAD. SERPINE1 serves as a primary suppressor of urokinase plasminogen activator and tissue plasminogen activator and downregulating SERPINE1 can hamper the proliferation and migration of GC cells33. In our study, SERPINE1 was identified as a hazard factor for the prognosis of STAD patients (HR > 1), consistent with previous research showing that SERPINE1 is a procancer factor. Cocaine and amphetamine-regulated transcript (CART), encoded by the CARTPT gene, is considered a ligand for G protein-coupled receptor 16034. Recent studies have shown that in prostate cancer and glioblastoma cells, GPR160 overexpression promotes tumor cell proliferation, while GPR160 knockdown has the opposite effect35,36. OPN3 belongs to the opsin subfamily within the G protein-coupled receptor superfamily, and upregulation of OPN3 expression facilitates epithelial-mesenchymal transition and tumor migration in lung cancer37. The above studies suggest that the CARTPT and OPN3 genes act as procancer factors in several cancers, but their role in STAD is not yet clear. TP53 is a crucial gene for suppressing tumor formation and is the most frequently mutated gene in most cancer types38. In our study, TP53 was identified as a protective factor (HR < 1), consistent with previous research39. Additionally, the RT-qPCR results indicated that the expression levels of OPN3 was elevated in tumor tissues compared to adjacent tissues, while TP53 expression was found to be decreased. However, it is important to note that CARTPT and SERPINE1 did not show a statistically significant difference in expression levels, which constitutes a limitation of this study. Previous research had identified increased expression of CARTPT and SERPINE1 in STAD40,41. Therefore, we propose that the following factors may contribute to this discrepancy: (1) The relatively small sample size; (2) The presence of various confounding factors among different patient samples, which may lead to considerable individual variability. In future research, we aim to further investigate the expression and functional roles of CARTPT and SERPINE1 in the context of STAD.

Tumors with a high mutation burden are generally regarded to be linked to a good prognosis. This can be predominantly attributed to heightened, tumor-infiltrating immune cells and tumor-specific antigens42. According to our model, the low-risk group manifested a remarkable TMB than did the high-risk group, and the low TMB group exhibited a lower survival rate than did the high TMB group. However, specific genes may show opposite results. Interestingly, in our study, patients harboring TP53 mutations had a greater risk score than those with wild-type, possibly because TP53 mutation leads to decreased apoptosis, increased DNA repair, and increased tumor cell growth and migration, leading to poor prognosis in STAD patients43.

Solid tumors constitute predominantly of the extracellular matrix, tumor cells, immune cells, and stromal cells. Previous literature has implied that decreased tumor purity is markedly linked to a mesenchymal, invasive, and metastatic phenotype and is associated with a shorter survival time44. The low-risk group affirmed a greater tumor purity score, with a corresponding increase in the ESTIMATE score as per our findings. Immune cells have a crucial function in the TME. Disruption of the circadian rhythm alters the cellular composition of TME, elevating the ratio of immunosuppressive immune cells encompassing M2 macrophages and MDSCs, while diminishing the proportion of antitumor cells like CD4 + T cells, CD8 + T cells, and M1 macrophages, leading to an immunosuppressive microenvironment that favors tumor immune escape and proliferation6. According to our CRG signature, the abundance of M2 macrophages in the high-risk group was greater than that in the low-risk group, whereas the presence of CD8 + T cells was reduced. We established that patients belonging to the high-risk group manifested a shorter survival duration, indicating that immune infiltration is a key factor influencing the prognosis of STAD patients as per the K‒M curve analysis. To elucidate the function of CRG in STAD immunity, we further analyzed immune checkpoints. Notably, in our study, PDCD1LG2 expressions were elevated levels in the high-risk group. PD-L2, which is encoded by the PDCD1LG2 gene, is the second most significant ligand discovered to interact with PD-1 following PD-L1. The PD-L2 and PD-1 interactions suppresses the PD-1 signaling pathway, promotes the apoptosis of cytotoxic T cells, facilitate tumor cells to escape immune surveillance, and thereby facilitates rapid proliferation45.

Despite the significant benefits of ICIs in several advanced cancers, their efficacy in STAD patients is concerning46. It is a significant challenge to identify STAD patients who will positively respond to ICI therapy. Prior research has indicated that a higher TIDE score is linked to an increased likelihood of immune evasion by cancer cells47. The IPS integrates antigen presentation, effector cell, suppressive cell, and checkpoint scores, and a higher IPS signifies heightened immunogenicity and a better response to immunotherapy48. According to our study, compared to the high-risk group, the low-risk group manifested lower TIDE scores and higher IPS scores, indicating a more favorable response to ICIs in low-risk patients. Interestingly, in predicting the sensitivity of STAD to chemotherapeutic drugs, the high-risk group affirmed remarkably greater sensitivity than the low-risk group. The results further infer that our CRG signature has ability to anticipate the response of STAD patients to both ICIs and chemotherapy, which is extremely important for guiding the treatment of patients with advanced STAD.

Our study also identified several important biological processes associated with CRG and high recurrence risk. The involvement of the PI3K-Akt signaling pathway in cancer initiation, progression, and metastasis is widely recognized. Consistent with our study, targeting the PI3K-Akt signaling pathway can impede the progression of STAD49,50. cAMP and Ca2+ are important cellular second messengers that have long been recognized as key regulatory factors in cell growth and invasion. Previous literature implies that modulating the cAMP signaling pathway as well as the calcium signaling pathway can affect the proliferation and migration of STAD cells51,52.

Through PPI network analysis of DEGs, we discovered 10 hub genes. FN1 exhibited the highest level of connectivity in the PPI network. In our study, by qPCR, we confirmed that FN1 expression is elevated in cancer tissues compared to that in adjacent tissues. The findings of survival analysis affirmed that the group with lower FN1 expression exhibited a comparatively greater rate of survival that the group with higher FN1 expression. Immune cell correlation analysis revealed a negative correlation between the expression of FN1 and the activation of dendritic cells and plasma cells, while a positive correlation was found between the expression of FN1 and M2 macrophages. Activated dendritic and plasma cells play pivotal roles in antigen presentation and B cell activation, whereas M2 macrophages are primarily associated with immune suppression. Utilizing these findings, we can conclude that FN1 can influence the TME and is closely linked to poor patient prognosis, consistent with previous research53. Here, we discovered that the CRG signature is linked to procancer signaling pathways and the FN1 gene, suggesting that disruption of the circadian rhythm plays a significant function in STAD progression.

Although the constructed signature of CRG has demonstrated promising utility in STAD patients’ prognoses and aiding in treatment decisions, it is paramount to take into consideration the delimitations that affected our study. First, our study relies on existing public databases. Although the constructed signature exhibited satisfactory performance in the TCGA and GEO datasets, its predictive capacity requires further validation through randomized controlled clinical trials. Second, the hub genes as well as signaling pathways identified in this study using bioinformatics approaches need further experimental validation to elucidate their underlying mechanisms. Third, the RT-qPCR results showed no significant difference in the expression of CARTPT and SERPINE1 between adjacent tissues and tumor tissues, which may be related to the relatively small number of clinical samples we collected. Finally, further clinical validation is vital to verify the accuracy and reliability of our prediction regarding the response to immunotherapy in STAD patients, particularly considering that it is based on the correlation between the risk score and IPS.

Conclusion

In summary, we have achieved success in developing a prognostic signature for CRG, which exhibits a high level of effectiveness in accurately predicting both the prognosis as well as the immunotherapeutic response for STAD patients. These provide insights into the development of personal-centered treatment interventions for STAD patients.