Identifying and Validating Circadian Rhythm Biomarkers

Introduction Circadian rhythm refers to the physiological oscillation system present in living organisms, characterized by a cycle of approximately 24 hour.1 This rhythm plays a crucial role in regulating a wide array of physiological, biochemical, and behavioral processes across diverse organisms, ranging from unicellular entities to humans, thereby facilitating synchronization with the cyclical alternation of day and night on Earth.2,3 In mammals, the fundamental mechanism underlying circadian rhythm is governed by a network of clock genes that operate through transcriptional-translational feedback loops (TTFLs).4,5 The heterodimer composed of circadian locomotor output cycles kaput (CLOCK) and brain and muscle ARNT-Like 1 (BMAL1) functions as an activating factor, whereas the period (PER) and cyptochrome (CRY) protein families serve as inhibitory factors, thereby sustaining the stability of the circadian rhythm through their cyclical expression and degradation.3,6,7 These genes are essential for maintaining homeostasis within the body by coordinating critical physiological functions such as sleep-wake cycles, metabolism, immune responses, and DNA repair.3,8 In contemporary society, various factors—including shift work, transmeridian travel, and exposure to artificial lighting—can disrupt circadian rhythms, leading to circadian rhythm disruption (CRD).9 Such disruptions have been associated with an increased risk of metabolic syndrome, neurodegenerative diseases, and cancer.10,11 Disruption of circadian rhythms contributes to the onset and progression of various diseases through multiple mechanisms, including the disturbance of the balance between cell proliferation and apoptosis, impairment of immune surveillance, and interference with metabolic pathways. In the context of metabolic-related diseases, the association between circadian rhythm disruption and conditions such as obesity and diabetes is well-established.12,13 With regard to cardiovascular diseases, chronic disruptions of circadian rhythms increase the risk of several conditions, including coronary artery disease, myocardial infarction, heart failure, atrial fibrillation, and stroke.10,14 In the realm of neurological disorders, patients with Parkinson’s disease and Alzheimer’s disease often exhibit abnormal expression of clock genes.15 Moreover, the circadian rhythm system is intricately linked to the onset, progression, and therapeutic response of cancer.11 For instance, implementing a 6-hour time-restricted feeding (TRF) regimen has been shown to reprogram the circadian rhythm and glycolytic metabolism associated with clock genes, thereby inhibiting the progression of lung cancer in murine models.16 These findings underscore the pivotal role of clock genes in the regulation of malignant tumors. Within the TTFL of the circadian rhythm, the nuclear receptor family of retinoic acid-related orphan receptors (RORs) functions as a pivotal positive regulator. This family establishes an antagonistic feedback mechanism with the REV-ERB family to uphold circadian stability.17 The ROR family comprises three subtypes: RORα, RORβ, and RORγ, with RORγ further differentiated into the isoforms RORγ1 and RORγt, both encoded by the RORC gene.18,19 RORC specifically binds to RORE sequences within the PER1, PER2, CRY1, and BMAL1 genes, thereby enhancing their transcription to stabilize negative feedback and ensure circadian accuracy.20 Additionally, it regulates clock control genes (CCGs) to convey central clock signals to peripheral tissues, thereby synchronizing the rhythms of various organs. Beyond its role in circadian regulation, RORC is essential for the differentiation of immune cells, such as Th17 cells,21 and plays a significant role in glucose-lipid metabolism and reproductive development.22 Thus, RORC serves as a crucial nexus linking circadian rhythms to a wide array of physiological and pathological processes. Gastric cancer (GC), characterized as a malignant tumor with one of the highest mortality rates globally, continues to pose significant challenges in clinical diagnosis and treatment due to its high incidence and poor prognosis.23 While risk factors such as Helicobacter pylori infection and irregular dietary habits have been extensively investigated,24,25 the mechanistic link between circadian rhythm and gastric cancer remains in the exploratory phase.26 Recent studies suggest that individuals engaged in long-term night shift work are at an elevated risk of developing gastric cancer relative to the general population.27 Additionally, aberrant expression of clock genes has been linked to the malignant progression of gastric cancer cells. In cases of HER2-positive gastric cancer that exhibit resistance to trastuzumab, glycolysis has been observed to vary in accordance with circadian oscillations, which are regulated by the BMAL1-CLOCK-PER1- hexokinase 2 (HK2) axis.26 Nevertheless, these investigations primarily concentrate on individual genes or are limited by small clinical sample sizes, thereby lacking a comprehensive analysis of the overall expression patterns of circadian rhythm-related genes (CRGs) and their association with gastric cancer prognosis, particularly in the context of identifying potential prognostic markers. In light of these considerations, this study employed transcriptomic data of GC obtained from publicly available TCGA and GEO databases and integrated it with CRGs. We identified the principal genes associated with rhythm regulation in patients with gastric cancer and conducted comprehensive bioinformatics analyses to examine the key genes. Furthermore, we investigated their biological functions within the context of gastric cancer. This research not only advances the understanding of the molecular mechanisms underlying circadian rhythm disruptions in the progression of GC but also lays a theoretical foundation for the development of novel prognostic assessment tools and targeted therapeutic strategies. Materials and Methods Data Acquisition and Processing The gene expression and clinical data of gastric cancer patients were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/), comprising 36 normal tissue samples and 412 gastric cancer tissue samples. Additionally, data were sourced from the GSE79973 dataset within the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/), which includes 10 samples of gastric cancer tissues and 10 samples of adjacent normal tissues. Based on previous studies and reviews (doi:10.3389/fimmu.2022.777724; doi:10.1038/nrg2430), we identified circadian rhythm-related genes (CRGs) to elucidate the core regulatory role of circadian rhythms in various biological processes. Seventeen CRGs were selected for further analysis, namely ARNTL, ARNTL2, CLOCK, CRY1, CRY2, CSNK1D, CSNK1E, NPAS2, NR1D1, NR1D2, PER1, PER2, PER3, RORA, RORB, RORC, and TIMELESS. The selection of these 17 CRGs is grounded in their established roles in circadian rhythm regulation and cancer research. Differential Analysis The analysis of differentially expressed genes (DEGs) between cancerous and adjacent normal tissues within the TCGA gastric cancer database was conducted utilizing the “DESeq2” package. Similarly, the DEGs between gastric cancer tissues and adjacent normal tissues in the GSE79973 dataset were analyzed using the “limma” package. The threshold for significance was set at an adjusted P< 0.05 and |log2 FC| > 0.585. Visualization of these DEGs was performed through the generation of a volcano plot using the “ggplot2” package (version 3.4.4). Additionally, a heatmap illustrating the upregulated and downregulated genes was created using the same “ggplot2” package. Candidate Genes Acquisition and Enrichment Analysis Integration of the TCGA, GSE79973, and CRG datasets was achieved through the “ggvenn” software package to identify candidate genes. Subsequently, functional enrichment analysis of these candidate genes was performed using the “clusterProfiler” software package, encompassing both Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses. The results of these analyses were visualized using the “enrichplot” R package. Immune Cell Estimation and Gene Set Variation Analysis (GSVA) In this study, the initial step involved configuring the working directory and loading pertinent R packages, including GSVA. Subsequently, standardized expression data and immune-related gene set files were imported. The ssGSEA analysis was performed using the GSVA package to derive immune-related pathway scores within cancer samples. Following standardization, the results were preserved. The ssGSEA result data were then accessed, and expression data for the RORC gene were extracted and sorted based on expression values, categorizing them into high and low expression groups. Concurrently, data for other immune-related genes were extracted, sorted, and standardized in the same sequence. Ultimately, a bar chart depicting RORC expression (ordered by expression value) and the Z-score heat map of the other genes were generated, combined, and saved as a PDF. Additionally, for each immune-related gene and RORC, a composite graph was created, comprising a scatter plot, a top gene expression histogram, and a right RORC histogram. Correlations were calculated and annotated, and the final summary was saved as a PDF. Drug Sensitivity Analysis This study employed the “pRRophetic” package to predict the sensitivity of cancer samples exhibiting varying levels of RORC expression (categorized into high and low expression groups) to multiple drugs, as quantified by IC50 values. The differences in drug sensitivity between these groups were statistically analyzed using the Wilcoxon test. Genetic Alteration Analysis This study conducted an analysis of the gene mutation characteristics in tumor samples. The samples were stratified based on the expression levels of RORC genes, and the “maftools” package was employed to perform a visual analysis of the mutation data for each group. Somatic mutation data were extracted from the TCGA database, and mutation-related information was systematically screened and organized, encompassing mutant genes, mutation types (such as missense mutations and frameshift mutations), and mutation frequencies. The “plotmafSummary” function was utilized to generate a mutation summary chart, illustrating the mutation burden, distribution of mutation types, and high-frequency mutant genes within each group. Additionally, a waterfall plot was created using the “oncoplot” function to depict the mutation status of the top 20 high-frequency mutant genes in the samples, facilitating a comparison of mutation spectra differences across the various groups. RT-qPCR The expression levels of genes were assessed using reverse transcription quantitative polymerase chain reaction (RT-qPCR). Total RNA was extracted utilizing the TRIzol reagent (TaKaRa, Japan), and complementary DNA (cDNA) synthesis was conducted with the TaKaRa reverse transcription kit (TaKaRa, Japan). RT-qPCR was carried out employing the SYBR Premix Ex Taq II (TaKaRa, Japan), with β-actin serving as the reference gene for data normalization. Primer sequences for β-actin were as follows: 5′-ATAGCACAGCCTGGATAGCAACGTAC/CACCTTCTACAATGAGCTGCGTGTG-3′. Primer sequences for RORC were as follows: 5′- GAGGAAGTGACTGGCTACCAGA/GCACAATCTGGTCATTCTGGCAG-3′. Cell Lines and Cell Culture The human gastric cancer cell lines, AGS and HGC-27, were procured from Procell (Wuhan, China). AGS cells were maintained in F12 medium (Gibco, USA) supplemented with 10% fetal bovine serum (FBS; Gibco, USA) and 1% penicillin-streptomycin (P/S), under conditions of 5% CO2 and 37 °C in a humidified incubator. Similarly, HGC-27 cells were cultured in RPMI-1640 medium (Gibco, USA), also supplemented with 10% FBS and 1% P/S, under identical incubator conditions. Cell Transfection The cells were inoculated into 6-well plates and allowed to adhere to the surface. Cell transfection was conducted using Lipofectamine™ 2000 reagent (Invitrogen, USA), with pLVX-IRES-puro-RORC and empty vector control plasmids introduced into the cells, respectively. Following a 6-hour incubation period at 37°C, the complete culture medium was added to the cells. Cell Viability Assay To assess cell viability, the Cell Counting Kit-8 (CCK-8) assay was employed. Cells were seeded into a 96-well culture plate at a density of 1000 cells per 100 microliters per well. Subsequently, 10 microliters of CCK-8 reagent were added to each well, and the plate was incubated at 37°C for 2 hours. Following incubation, the optical density was measured at a wavelength of 450 nm. Colony Arrangement Assay To assess the proliferative capacity of cells, plate cloning experiments were conducted. A total of 2000 cells per well were inoculated into a 6-well plate and incubated at 37°C with 5% CO2 for a duration of two weeks. Subsequently, the cells were fixed and stained using 4% paraformaldehyde and crystal violet. Cell Migration and Invasion Analysis For migration assays, cells were resuspended in a serum-free medium and introduced into the upper chamber of a 24-well Transwell plate (pore size 8 μm, Corning, USA) at a density of 4 × 104 cells per well. The lower chamber was supplemented with culture medium containing 10% FBS. For invasion assays, the chamber was pre-coated with matrix gel (BD, USA). After a 24-hour incubation period, cells were fixed and stained with 4% paraformaldehyde and crystal violet. Statistical Analysis Version 4.2.2 of the R software was employed for bioinformatics analyses. Each experiment was independently conducted in triplicate. All experimental data are reported as mean ± standard deviation (SD). Data analysis was performed using Prism 10.0 (GraphPad Software). A P-value of less than 0.05 was considered indicative of a statistically significant difference. Result Differentially Expressed mRNA in Gastric Cancer To identify significant differentially expressed genes influencing the progression of gastric cancer, we performed an analysis utilizing multiple public databases. The data were standardized by eliminating batch effects from the gene chip data using “sva” R package. By applying the criteria of |logFC| > 0.585 and P < 0.05, we identified 3,507 differentially expressed mRNAs in the GSE79973 dataset, comprising 1,863 up-regulated and 1,644 down-regulated mRNAs. We employed R to generate volcano plots of the differentially expressed mRNAs, with significant genes highlighted in distinct colors: up-regulated genes in red and down-regulated genes in blue (Figure 1A). We also plotted the top 50 mRNAs with high significance (P < 0.01) and notable expression changes in GSE79973 (Figure 1B). In the TCGA database, we identified 7,745 differentially expressed genes, of which 4,261 were up-regulated and 3,484 were down-regulated (Figure 1C). A Venn diagram was utilized to illustrate the genetic intersections among the GSE79973, TCGA data, and CRG, revealing the differentially expressed gene RORC (Figure 1D). The expression differences of the RORC gene in gastric cancer and normal tissues were analyzed by Wilcoxon rank sum test, showing that the RORC gene was significantly lowly expressed in gastric cancer tissues, with p values all less than 0.05 (Figure 2A and B). The GO and KEGG Functions of Differentially Expressed Genes are Enriched The differentially expressed genes were subjected to GO and KEGG pathway analyses. The GO analysis indicated that these genes were significantly enriched in biological processes (BP) such as the positive regulation of circadian rhythm, cellular response to sterol, and T-helper 17 cell differentiation. Furthermore, they were associated with molecular functions (MF) including nuclear receptor activity, ligand-activated transcription factor activity, and sterol binding (Figure 3A). The KEGG pathway analysis further identified significant enrichment of these genes in pathways related to circadian rhythm, inflammatory bowel disease, and Th17 cell differentiation (Figure 3B). Immune Correlation Analysis We identified the expression of the RORC gene and its correlation with immune-related markers. The upper bar chart presents the RORC gene expression levels following Log2 transformation. Samples are arranged in ascending order of expression values and categorized into low (blue) and high (red) expression groups, with a dotted line demarcating the two groups. The heat map below illustrates the Z-score expression of various immune-related markers (eg, Mast cells, Type II IFN Response), maintaining the same sample order as the bar chart. The markers are annotated on the right side of the heat map, with some indicators denoted for statistical significance (*p<0.05, **p<0.01, ***p<0.001). Notably, several markers, such as Tfh, T cell co-stimulation, and aDCs, exhibit significant correlations with RORC expression. The Z-score ranges from −4 to 4, indicating the differential expression of immune markers between the high and low RORC expression groups (Figure 4). Analysis of RORC Gene-Related Mutations The mutation summary chart indicates variability in mutation burden across the different sample groups, with median values of 109 and 94, respectively (Figure 5A and B). Missense mutations constitute the predominant mutation type, exhibiting the highest proportion. Within the top 10 mutant genes, TTN, TP53, MUC16, and LRP1B demonstrate notably high mutation frequencies. Specifically, the mutation rate for TTN ranges from 49% to 54%, while TP53 exhibits a mutation rate between 43% and 49%. The waterfall plot further elucidates the distribution of the top 20 mutant genes across each sample group, with distinct colors denoting various mutation types (eg, Frame_Shift_Del, Nonsense_Mutation) (Figure 5C and D). Notably, key genes such as TP53 and ARID1A display relatively high mutation frequencies in both groups. Furthermore, disparities in the mutation distribution of certain genes, such as SYNE1 and CSMD3, between the groups suggest that the gene mutation pattern may be associated with the expression level of RORC genes. Analysis of Drug Sensitivity Related to RORC Genes Using a threshold of p<0.05 for screening, we identified drugs that exhibited a significant association between RORC expression levels and drug sensitivity. For these drugs, we illustrated the distribution of drug sensitivity in the RORC high-expression group (depicted in red) and the low-expression group (depicted in blue) using box plots. The statistical comparisons between the groups were annotated within the plots. The box plots for the significantly associated drugs were labeled as “urgsenstivity.” The analysis revealed that the IC50 values for VX-680, MG-132, Sunitinib, Paclitaxel, Imatinib, Midostaurin, BMS-754807, GNF-2, TAE684, S-Trityl-L-cysteine, XMD8-85, Parthenolide, Lisitinib, BMS-509744, 5-Fluorouracil, Rapamycin, Pyrimethamine, A-443654, and CGP-082996 were lower in the low-expression group compared to the high-expression group (Figure 6A–T). This suggests that these drugs may be more effective in treating patients with RORC gene deficiency. RORC Attenuates the Proliferation and Migration of GC Cells To elucidate the biological function of RORC in these cells, we performed experiments involving the transfection of GC cell lines to induce upregulation of RORC expression (Figure 7A). Compared to the control group, RORC overexpression significantly reduced the proliferation of GC cells (Figure 7B and C). Furthermore, invasion and migration assays demonstrated that RORC markedly diminished the invasive and migratory capacities of GC cells (Figure 7D). In conclusion, the findings suggest that RORC plays a role in inhibiting the progression of GC cells. Discussion The circadian rhythm is integral to the regulation of metabolic homeostasis, behavioral modulation, and immune function within the body. Central to the circadian mechanism are core clock genes, such as BMAL1 and CLOCK, which form heterodimers that interact with promoter regions of downstream genes, thereby facilitating the transcription of PER and CRY genes.6 PER and CRY proteins accumulate in the cytoplasm, subsequently translocate to the nucleus, and inhibit the activity of the BMAL1-CLOCK complex.3 This inhibition constitutes a negative feedback loop that regulates their own transcription, culminating in oscillatory expression patterns with a periodicity of approximately 24 hours.3,6 In recent years, with the advancement of research, the significance of the retinoic acid-related orphan receptor (ROR) family in the regulation of circadian rhythms has garnered increasing attention. The ROR family comprises three members: RORα, RORβ, and RORγ (encoded by RORC), which are integral components of the nuclear receptor superfamily and are involved in various physiological processes, including cellular differentiation, metabolism, and immune regulation.18,19 Within the circadian rhythm regulatory network, ROR family members modulate the expression of associated genes by binding to specific DNA response elements, thereby influencing the maintenance and regulation of circadian rhythms.19 As a pivotal member of the ROR family, RORC encodes two isoforms: RORγ1 (also referred to as RORc1 or RORγ), which is broadly expressed across various tissues, and the thymus-specific RORγ2 (also known as RORc2 or RORγt).28 The RORC protein comprises a DNA binding domain (DBD) and a ligand binding domain (LBD), enabling it to bind to the ROR response element (RORE) in a monomeric form and regulate downstream gene transcription. In the context of circadian rhythm regulation, RORγ directly influences the transcription of several clock factors (such as BMAL1, Npas2, and CRY1) by binding to ROREs,20 thereby playing a crucial role in modulating the expression of genes involved in lipid, glucose, and steroid metabolism.29 In triple-negative breast cancer, RORγ serves as a key driver of the cholesterol biosynthesis program. The application of RORγ-selective antagonists has been shown to induce significant tumor regression and inhibit metastasis in TNBC models.30 Furthermore, RORγ is overexpressed in osteosarcoma patients, where it promotes the mitochondrial oxidative phosphorylation program, positioning it as a viable therapeutic target for advanced OS.31 The primary role of thymus-specific RORγt is to function as a crucial transcription factor in the differentiation of Th17 cells, thereby contributing to autoimmune and inflammatory responses through the regulation of factors such as IL-17. Aberrant methylation patterns can lead to abnormal RORC expression levels in various cancers, and its dysregulation is significantly associated with cancer initiation, progression, and prognosis.32 A comprehensive pan-cancer analysis revealed that RORC is highly expressed in six types of tumors, including breast cancer and lung adenocarcinoma, while it is underexpressed in fourteen types, such as renal clear cell carcinoma and hepatocellular carcinoma. Studies indicate that RORC is overexpressed and amplified in metastatic castration-resistant prostate cancer, where it drives the expression of androgen receptors within the tumor.33 In bladder cancer, reduced RORC expression enhances glycolysis and chemoresistance by activating the PD-L1/ITGB6/STAT3 signaling axis.34 The mechanism of RORC action varies across different cancers, and its specific molecular regulatory network warrants further detailed investigation. In this study, we conducted a novel investigation into the role and potential mechanisms of RORC in gastric cancer. Our findings indicate that the expression level of RORC is significantly reduced in gastric cancer tissues compared to adjacent normal tissues. Functional assays demonstrated that RORC overexpression inhibits the proliferation, migration, and invasion of gastric cancer cells. These findings propose that RORC may function as a potential tumor suppressor gene in the regulation of gastric cancer initiation and progression. Notably, this level of inhibition from single-gene overexpression is relatively uncommon when compared to the moderate (30–50%) changes generally seen with the knockdown of classical oncogenes like MYC and RAS. This observation suggests the possibility that the inhibitory effects of RORC overexpression may involve non-specific cytotoxicity or the induction of apoptosis, rather than being solely mediated by a targeted anti-proliferative mechanism. Future research should aim to delineate the contributions of specific anti-tumor signaling pathways versus non-specific cytotoxic effects associated with RORC overexpression, thereby elucidating the precise molecular basis of its functional role. In conclusion, this study investigates the circadian rhythm-associated gene RORC and its involvement in gastric cancer. Our findings indicate that RORC acts as a tumor suppressor in gastric cancer, with its reduced expression correlating with disease progression and exerting inhibitory effects on the malignant characteristics of gastric cancer cells. These results suggest that RORC could serve as a circadian rhythm-related biomarker for the diagnosis and prognosis of gastric cancer, as well as a potential therapeutic target for the development of personalized treatment strategies. However, this study has certain limitations. Although we have validated the tumor-suppressive role of RORC in gastric cancer, further research is required to elucidate the precise molecular mechanisms underlying its regulatory function in cancer progression and to assess its clinical applicability in larger patient cohorts and prospective studies.