This article is mentioned in:

Introduction

According to the World Cancer Report 2020, ovarian cancer (OV) represents the eighth most common cancer affecting females worldwide (1). Due to difficult detection at the early stage, OV has a high fatality rate, i.e., a 5-year survival rate of ~30%, and poor outcomes in patients. Despite years of great breakthroughs in diagnosis and treatment, OV remains the deadliest malignancy in females globally, and its incidence and death rates generally increase with age (1,2). OV presents a major therapeutic challenge due to its high heterogeneity, aggressive nature and lack of efficient targeted therapies (3). Furthermore, the pathogenesis of OV is complex and requires a deeper understanding of genetics as well as potentially modifiable risk factors such as circadian rhythms.

A large number of organisms possess a biological clock organizing oscillations in physiological and behavioral events, which is referred to as circadian rhythms (4–6). Circadian rhythms are generated by an intracellular clock mechanism that is involved in long-term biological evolution, representing an essential feature of life activity (7–9). Emerging evidence suggests that circadian clock genes (CCGs) and downstream effectors control multiple cancer-associated biological events, including metabolism, inflammatory responses, DNA damage repair and cell cycle (10). Therefore, studying the association of CCG dysregulation with tumor progression is critical for developing effective clinical strategies. Of note, CCG dysregulation was indicated to be tightly associated with carcinogenesis in multiple cancer types. In females, disordered circadian rhythms resulting from unfavorable work shifts or other stressors were observed to induce menstrual disorders, as well as breast cancer and OV (11,12). Hence, an overall understanding of CCGs may provide important insight into tumor biology and help develop effective models for OV prognosis.

The present study was the first to utilize OV cohorts from both the Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases to explore the association of circadian clock aberrations resulting from imbalanced circadian rhythm with OV to identify reliable circadian biomarkers and novel predictive and therapeutic targets. First, CCGs that were able to predict poor outcome in patients with OV were identified. Subsequently, reticulon 2 (RTN2) gene expression was assessed in both normal and OV tissues. A flow chart of the present study is provided in Fig. 1.

Figure 1. Flow chart of the hub circadian genes screening process. OV, ovarian cancer; CCG, circadian clock gene; TCGA, The Cancer Genome Atlas.

Materials and methods

Data retrieval and curation

A total of 1,409 CCGs were retrieved from the Circadian Gene DataBase (CGDB), which indicates alterations in mRNA amounts of these CCGs confirmed by previously published reports by reverse transcription-quantitative (RT-q)PCR, northern blot and in situ hybridization (13). Subsequently, gene expression profiles for 53 OV and 10 normal ovary tissue specimens were retrieved from GEO (https://www.ncbi.nlm.nih.gov/geo; accession no. GSE18520). Gene expression analysis was performed with the Affymetrix Human Genome U133 Plus 2.0 Array (Affymetrix; Thermo Fisher Scientific, Inc.). Raw data processing utilized the Robust Multi-Array Average (RMA) method and the ‘Oligo’ package from Bioconductor (http://www.bioconductor.org) for data normalization and probe annotation. Another public OV dataset with 57 tumor and 12 normal tissue samples, GSE66957, obtained with the Rosetta/Merck Human RSTA Custom Affymetrix 2.0 microarray, was downloaded. The matrix data and GPL10379 files were utilized for generating normalized data and annotated probes.

Analysis of CCG expression

The R package limma was utilized for differential expression analysis of CCGs between OV and noncancerous tissue specimens (14), and genes with the P-value of the false discovery rate (FDR) <0.05 were defined as differentially expressed. A Venn diagram containing four lists of differentially expressed genes was drawn online (https://bioinfogp.cnb.csic.es/tools/venny/index.html) to identify those genes that were overexpressed in these datasets.

Survival analysis

The effects of CCGs on survival of patients with OV were assessed in the GSE49997 dataset and validated in the TCGA OV cohort. The whole gene expression profile data and associated OV patient features were obtained from GEO and TCGA (https://cancergenome.nih.gov/). Kaplan-Meier curves were established for assessing the associations of overall survival (OS) with CCG expression levels using the log-rank test. The optimal cutoff of the gene expression value was determined with the survminer package (v0.4.6,https://github.com/kassambara/survminer/) based on expression levels, survival time and survival status (15).

Protein-protein interaction (PPI) network construction

In total, 217,249 PPI pairs were retrieved from Reactome (v.2014; http://www.reactome.org) (16), based on BioGrid, the Database of Interacting Proteins (17), Human Protein Reference Database (18), I2D (19), IntACT (20) and MINT (21), in addition to gene co-expression data generated by high-throughput techniques such as yeast two-hybrid, mass spectrometry pull-down and DNA microarray assays (22). A PPI network was generated with Cytoscape (v.3.2.1; http://www.cytoscape.org) (23).

Pathway enrichment analysis for the functional interaction network

Reactome FIViz was utilized in Cytoscape for pathway analysis (24). Cell Map (http://www.pathwaycommons.org/pc/dbSnapshot.do?snapshot_id=8), Reactome (16), Kyoto Encyclopedia of Genes and Genomes (25), Panther Pathways (26), NCI-Pathway Interaction Database (NCI-PID) (27) and BioCarta (http://www.biocarta.com/genes/index.asp) were used as pathway annotation sources, with an FDR of 0.05 as the cut-off criterion.

Construction of the prognostic signature and calculation of the risk score

For the 9 validated prognostic CCGs, the Cox proportional hazard regression model was used to construct prognostic models using the GSE49997 cohort. An equation was established for calculating risk score as follows: Risk score=βULK1 × expressionULK1 + βCRY2 × expressionCRY2 + βGAS7 × expressionGAS7 + βNPTXR × expressionNPTXR + βRTN2 × expressionRTN2 + βATF3 × expressionATF3 + βCSF3R × expressionCSF3R + βDAAM2 × expression DAAM2 + βPPP1R15A × expression PPP1R15A. β was an equation coefficient. Risk scores were then determined for all patients in the TCGA OV cohort. With the median risk score as the cutoff, the cases were assigned to the high- and low-risk groups. OS times in both risk groups were analyzed by the Kaplan-Meier method using the log-rank test for comparison. The area under the curve for the survival receiver operating characteristic (ROC) curve was determined with the survival ROC function in R software for validating the performance of the prognostic signature (28).

Patients and specimens

A total of 42 OV and 20 noncancerous tissue specimens (from different subjects than the OV group; 8 paracancerous tissues and 12 ‘healthy’ control tissues) were collected from patients enrolled at Ningbo First Hospital, Zhejiang University School of Medicine (Ningbo, China) between January 2015 and December 2021. The demographic and clinicopathological characteristics of the patients with OV are provided in Table SI. All subjects (patients with and without OV) who had undergone surgery provided written informed consent. OV samples were obtained from patients with no previous chemotherapy or radiotherapy prior to surgery. After surgical removal, tissue samples were frozen immediately in liquid nitrogen and stored at −80°C. The present study was approved by the Research Ethics Committee of Ningbo First Hospital (Ningbo, China; no. 2021-R210).

RT-qPCR

Freshly collected and then frozen tissue specimens were used for total RNA extraction with TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc.). Samples of total RNA with a 260/280 nm absorption ratio between 1.8 and 2.0 were further analyzed. cDNA was synthesized using a reverse transcription kit (Takara Bio Inc.) and amplified using a SYBR qRT-PCR Kit (Takara Bio, Inc.), as directed by the manufacturer. qPCR was performed in a Cobas z480 real-time PCR system (Roche Diagnostics), the qPCR reaction was performed at 95°C for 5 min, 60°C for 30 sec, followed by 40 cycles at 95°C for 30 sec and 58°C for 30 sec (29). The 2−∆∆Cq method (30) was employed for data analysis. GAPDH was used as an internal control. The following primers were used: RTN2 forward, 5′-GACCTGCTGTACTGGAAGGAC-3′ and reverse, 5′-ACGGACACGATGCTAAAGTGC-3′; GAPDH forward, 5′-AGGTCGGTGTGAACGGATTTG-3′ and reverse, 5′-TGTAGACCATGTAGTTGAGGTCA-3′.

Immunohistochemistry (IHC)

IHC staining was carried out according to standard procedures described in a previous report. Paraffin-embedded normal ovarian tissue and OV tissue samples were routinely processed and then incubated with rabbit primary antibodies against RTN2 overnight at 4°C (cat. no. 11168-1-AP; 1:200 dilution; Proteintech Group, Inc.). After a 1-h incubation with HRP-linked anti-rabbit secondary antibodies (cat. no. A0545; 1:4,000 dilution; MilliporeSigma) at room temperature, the DAB reagent was utilized for development, followed by hematoxylin counterstaining. A Leica DM 2000 microscope (Leica Microsystems) was used for analysis. Positive staining of RTN2 in tumor cells was assessed by IHC signal intensity. Scoring was conducted according to the ratio and intensity of positively stained cells: 0–5% scored 0; 6–35% scored 1; 36–70% scored 2; and >70% scored 3. According to the final score of RTN2 expression, a sample was designated as having low or high expression as follows: Low expression: Score 0/1; high expression: Score 2/3.

Statistical analysis

R v3.6.1 was employed for data analysis. Clinicopathological parameters in the high- and low-expression groups were compared by the Student's t-, χ2- and Fisher's exact tests, as appropriate. Univariate and multivariate Cox regression analyses were performed for determining factors independently predicting survival. P<0.05 was considered to indicate statistical significance.

Results

Expression of CCGs in OV

To assess the role of CCGs in OV, the gene expression profiles of 53 OV tumor samples and 10 normal ovary tissues from the GSE18520 dataset were compared and another public OV dataset with 57 tumor and 12 normal tissue samples (GSE66957) was examined. Heat maps (Fig. 2A and B) revealed gene expression differences between the tumor and non-tumor expression groups. Genes with the P-value of FDR <0.05 were considered to have differential expression. A total of 730 and 136 genes were upregulated and downregulated, respectively, in OV specimens compared with normal specimens in GSE66957. Similarly, 417 and 144 genes were upregulated and downregulated, respectively, in OV specimens compared with normal specimens in the GSE18520 dataset. Furthermore, Venn diagrams (Fig. 2C) revealed that in the two datasets, 264 and 15 genes were commonly upregulated and downregulated DEGs.

Figure 2. Cancer driver function of 279 circadian clock genes. (A) Heatmap of the DEGs in tumor vs. non-tumor tissues in the GSE18520 dataset with P(FDR)<0.05. (B) Heatmap of the DEGs in tumor vs. non-tumor tissues in the GSE66957 dataset with P(FDR)<0.05. (C) Venn diagrams presenting the number of upregulated or downregulated DEGs in tumor and normal groups (GSE18520 and GSE66957), generated by the Venn Diagram function. FDR, false discovery rate; DEG, differentially expressed gene.

Hub CCGs

In order to identify CCGs with active participation in OV carcinogenesis and progression, differentially expressed CCGs with a significant association with patient outcomes (P<0.05) were determined. PPI network analysis revealed hub genes in this dataset (Fig. 3A). Functional enrichment analysis indicated that the above genes were actively involved in the pathways such as circadian rhythm, CXCR4-mediated signaling events and rRNA processing (Fig. 3B). The hub genes were STAT3, CSNKIE, ARRB1, DVL3, MDM2, SOCS3, CDK4, NOP58, CBL, PSMD8 and TUBA4A (Fig. 3C). The above major nodes were associated with other genes, suggesting they may be able to affect the prognosis of OV.

Figure 3. PPI network of CCGs module. (A) PPI network analysis result of CCGs. (B) Pathway Enrichment Analysis result for the Functional Interaction Network. (C) The 11 hub genes in the PPI network. PPI, protein-protein interaction; CCG, circadian clock gene.

Association of CCGs with survival

In order to explore the prognostic value of individual DEGs, the associations of CCGs with OS in OV cases in the GSE49997 and TCGA databases were analyzed. A total of 9 validated prognostic CCGs were found in the GSE49997 OV discovery cohort and TCGA OV validation cohort, all of which were upregulated, with significant associations with reduced OS (months) (log-rank P<0.05), including ULK1, ATF3, CRY2, CSF3R, DAAM2, GAS7, NPTXR, PPPIR15A and RTN2 (Fig. 4A and B). In the GSE49997 OV discovery cohort, the associations of various DEGs in OS based on high (yellow line) and low (blue line) expression levels were confirmed (Fig. 4A). Furthermore, in the TCGA OV validation cohort, the associations of various DEGs in OS based on high (yellow line) and low (blue line) gene expression levels were also confirmed (Fig. 4B).

Figure 4. Prognostic value of 9 circadian clock genes in OV in the Gene Expression Omnibus and TCGA cohorts. (A) Kaplan-Meier curves for OS according to the expression of individual DEGs extracted from the comparison of groups in the GSE49997 OV discovery cohort. (B) Kaplan-Meier curves for OS according to the expression of individual DEGs extracted from the comparison of groups in the TCGA OV validation cohort. Yellow line, high expression; blue line, low gene expression group. OS, overall survival (months); DEG, differentially expressed gene; OV, ovarian cancer; TCGA, The Cancer Genome Atlas.

Validation of differentially expressed CCGs in the GEO cohort

In order to further confirm these data, public datasets were retrieved from GEO, including GSE18520 and GSE66957. First, the above 9 genes were examined in 53 OV and 10 noncancerous ovary tissue specimens in GSE18520. As indicated in Fig. 5A, the levels of ULK1, ATF3, CRY2, CSF3R, DAAM2, GAS7, NPTXR, PPPIR15A and RTN2 were elevated in OV samples compared with those in non-malignant tissues (P<0.05). In addition, the expression of the 9 genes in 57 OV and 12 noncancerous ovary tissue samples in GSE66957 was detected. As displayed in Fig. 5B, the levels of ULK1, ATF3, CRY2, CSF3R, DAAM2, GAS7, NPTXR, PPPIR15A and RTN2 were elevated in OV compared with those in the normal group (P<0.05).

Figure 5. Validation of different expression levels of 9 CCGs in Gene Expression Omnibus cohorts. (A) Expression levels of 9 validated prognostic genes in 53 tumorous samples and 10 normal samples in the GSE18520 dataset. (B) Expression levels of 9 validated prognostic genes in 57 tumorous samples and 12 normal samples in the GSE66957 dataset. The values in the box plots were presented as median and interquartile range. The relative gene expression normalized by the Robust Multi-Array Average method is presented on the y-axis. CCG, circadian clock gene.

Construction of a prognostic signature based on the 9 CCGs

From the 9 genes, ULK1, ATF3, CRY2, CSF3R, DAAM2, GAS7, NPTXR, PPPIR15A and RTN2, a prognostic signature model was generated, allowing for the calculation of a risk score for each patient according to the expression levels in the GSE49997 cohort (Fig. 6). The risk score may be an important tool for distinguishing among patients with OV based on potential discrete clinical outcomes (Fig. 6A and B). In Fig, 6C, the survival status of the patients is presented. The area under curve of the receiver operating characteristic (ROC) curve at 3-year was 0.66, suggesting moderate potential for the prognostic signature based on CCGs in survival monitoring (Fig. 6D). Kaplan-Meier curve revealed that patients in high-risk group had a shorter OS (P=0.006, Fig. 6E).

Figure 6. 9-gene signature model in the GSE49997 OV discovery cohort. (A) Rank of risk score of model and distribution of groups. (B) Heatmap of expression profiles of included genes. (C) Survival status of patients in different groups. (D) Survival-dependent ROC curve of prognostic value of risk score of model. (E) Patients in high-risk group suffered shorter OS. AUC, area under curve; OS, overall survival; ROC, receiver operating characteristic; OV, ovarian cancer; TP, true positive; FP, false positive.

Validation of the model based on 9 circadian clock-related genes in the TCGA dataset

The TCGA cohort with 353 OV patients was used for validation of the model (Fig. 7). The risk score also distinguished in TCGA OV cohort based on potential discrete clinical outcomes (Fig. 7A and B). OS of patients with alive and dead status is presented in Fig 7C. The area under the ROC curve at 3 years in the TCGA validation cohort was 0.56 and thus lower than that for the discovery cohort GSE49997, suggesting the model's external predictive power was limited (Fig. 7D). However, the Kaplan-Meier curve verified that high-risk cases among patients with OV had decreased OS (P=0.037, Fig. 7E).

Figure 7. 9-gene signature model in TCGA OV validation cohort. (A) Rank of risk score of model and distribution of groups. (B) Heatmap of expression profiles of included genes. (C) Survival status of patients in different groups. (D) Survival-dependent ROC curve of prognostic value of risk score of model. (E) Patients in high-risk group suffered shorter OS in validation cohort. TCGA, The Cancer Genome Atlas; AUC, area under curve; OS, overall survival; ROC, receiver operating characteristic; OV, ovarian cancer; TP, true positive; FP, false positive.

Confirmation by RT-qPCR and IHC analysis

The differential expression of the key gene RTN2 between OV and normal ovarian tissues was then experimentally validated in an internal cohort, there was no difference in age or BMI between the control and OV groups (Table SII). RTN2 mRNA levels were elevated in 42 OV specimens compared with 20 noncancerous tissue samples (P<0.05, Fig. 8). More importantly, compared with their normal counterparts, RTN2 expression determined by IHC was markedly increased in tumor tissues (P<0.001; Table I; Fig. 9). These data demonstrated that RTN2 was upregulated in OV.

Figure 8. RTN2 mRNA expression level in ovarian cancer tissues and normal ovarian tissues determined by reverse transcription-quantitative PCR. Data were presented as 2−ΔΔCq. *P<0.05. RTN2, reticulon 2.

Figure 9. Representative immunohistochemistry images for reticulon 2 protein expression (high, low, medium) in ovarian cancer tissues and normal ovarian tissues (magnification, ×200 in the left and ×400 in the right-hand panel; scale bar, 100 µm).

Table I. Association of RTN2 expression with OV determined by immunohistochemistry in the internal cohort.

Table I.

Association of RTN2 expression with OV determined by immunohistochemistry in the internal cohort.

	Expression of RTN2
Group	Low (n=23)	High (n=39)	Total	χ2	P-value
OV group	6	36	42	62	<0.001
Control group	17	3	20

[i] OV, ovarian cancer; RTN2, reticulon 2.

Discussion

OV is associated with poor prognosis, high malignancy and a rapid fatality potential. The reduced survival results from late presentation, early lymph node metastasis, common invasion of adjacent organs and poor chemotherapeutic response. The pathogenesis of OV remains to be fully elucidated and mostly involves malignant events, including apoptosis blockade, deregulated proliferation and gene mutation-induced differentiation, migration, adhesion, invasion and angiogenesis (2). For prognostic improvement in OV, effective therapeutic targets and predictive molecular markers should be developed.

The ovary is an important organ in females, which is both affected by and produces hormones. Hormone secretion is known to be strongly associated with OV occurrence and progression (31,32). Furthermore, female hormones are controlled by circadian rhythms, whose disruption increases the cancer risk. With the rapid pace of modern lifestyles, an increasing number of females, particularly professionals, ignore their normal biological rhythms in order to satisfy work demands (33). The ‘biological clock’ controls circadian rhythms in humans and is associated with multiple events in normal physiology as well as pathology. Epidemiological analyses suggested that unfavorable work shifts increase the risk of fatal OV (34,35). Most living organisms are active during the daytime, while resting at night, due to the biological 24 h rhythms controlled by a regular oscillation of the amounts of involved genes. Genes involved in circadian rhythms are highly expressed in the ovary for regulating ovulation, and disrupted circadian rhythms are associated with multiple risk factors for OV (36,37). Genomic data from ovarian tumors indicate that aberrant rhythmic changes may be involved in malignant biological behaviors (38). Hence, the present study aimed to examine CCGs and explore the mechanisms underlying tumorigenesis, and to further explore the associations of these genes for the development of OV therapies.

Exploring the mechanisms that regulate circadian rhythms helped identify CCGs, leading to a Nobel Prize award in 2017. In addition to controlling circadian rhythms, CCGs are also involved in multiple physiological and behavioral events, including sleep, feeding pattern, body temperature, hormone release and blood pressure (39,40). Interactions between tumor cells and disrupted circadian clock significantly contribute to cancer development (41). CCGs, particularly master tumor-related genes, are involved in the initiation, progression and evolution of OV (42,43). Multiple studies have reported that abnormalities in circadian rhythms are involved in various malignancies, including prostate (44), breast (42), endometrial (43), colorectal (45), liver (46) and lung (47) cancers, as well as leukemia (48). Therefore, identifying the abnormal expression of CCGs in OV is important as a biomarker for disease management. The analyses of the present study were performed based on public datasets, with the aim of expanding the current knowledge on how CCGs have roles in OV. Future studies by our group will further verify or validate these findings using in vitro or in vivo experiments.

In the present study, CCGs in OV cases were systematically assessed in GEO and TCGA datasets. Cancer drivers of 279 CCGs and the built PPI network further demonstrated that the assessed genes mostly contributed to cancer pathways, circadian rhythm pathways, CXCR4-mediated signaling events and rRNA processing. Furthermore, ATF3, CRY2, CSF3R, DAAM2, GAS7, NPTXR, PPPIR15A, RTN2 and ULK1 were identified as potential CCGs, which may be biomarkers for OV. Certain genes included in the present screening results have been previously reported (49–51) and are consistent with the present results, thus confirming the findings of the present study. For instance, the impact of ATF3 in tumorigenesis and immune cell infiltration of ovarian tumors was assessed in a previous bioinformatics study (49), ULK1 associated with progression-free survival in ovarian cancer, decreases autophagy and cell viability in high-grade serous ovarian cancer spheroids (50,51). However, various genes have not been in detail validated in vivo and in vitro. In the present study, a small sample validation in internal OV samples and controls was performed; based on the bioinformatics analyses and previous literature, RTN2 was expressed highly significantly in colon adenocarcinoma and gastric cancer (52,53), circadian clock genes were dysregulated and the expression levels have crucial roles in cancer (54). RTN2 was selected from the candidate genes for further experimental verification by RT-qPCR and IHC in OV samples; the results indicated that the differential expression was confirmed.

Given that the present study was based on bioinformatics, it had several limitations. First, it focused on CCG gene expression levels and their clinical significance. However, transcriptomics only reflects certain aspects but not all global alterations. Furthermore, data in various patients may have been obtained at distinct times, which may confound the analysis. Hence, whether these novel genes in combination would predict survival with higher potential in comparison with individual genes should be examined. Taken together, differentially expressed CCGs characterize asynchronous circadian rhythms in cancer and may represent a theoretical basis for chronotherapy.

In conclusion, the current study attempted to identify CCGs that contribute to OS in OV using the GEO and TCGA databases. Two independent cohorts were used to identify nine genes with potential utility in OV prognosis, including an independent OV discovery cohort and an OV validation cohort. Further investigation of CCGs controlling ovarian cell functions may help develop novel therapeutic targets for improving OV prognosis. These findings provide insight into the expression of the CCG RTN2 in clinical OV samples, which may provide information for further investigations of RTN2-associated mechanisms and drug development in OV.

Supplementary Material

Supporting Data

Acknowledgements

Not applicable.

Funding

This work was supported by the Natural Science Foundation of Zhejiang Province (grant no. LQ21H160011), the Medical and Health Plan of Zhejiang (grant no. 2021KY990, 2022KY310 and 2022KY1115), Ningbo Science and Technology Project (grant no. 2019F1003) and the Natural Science Foundation of Ningbo (grant no. 2019A610306 and 2019A610260).

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Author's contributions

XJZ, XL, LZ and JC designed the study. XJZ, XL and YH performed the experiments. XNZ, XL, YH and JC wrote and revised the manuscript. XL, LZ, JC and YH collected and analyzed the data. XJZ, LZ and JC confirm the authenticity of all the raw data. XJZ acquired funding. XNZ conceptualized and supervised the study. All authors have read and approved the final manuscript.

Ethics approval and consent to participate

This study was conducted in concordance with the tenets of the declaration of Helsinki and was approved by the Research Ethics Committee of Ningbo First Hospital (Ningbo, China; no. 2021-R210). Prior to enrollment, all patients provided written informed consent.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References