Genetic Landscapes and Prognostic Implications of Circadian Rhythm in Non-small Cell Lung Cancer

doi:10.21203/rs.3.rs-39450/v1

Download PDF

Research

Genetic Landscapes and Prognostic Implications of Circadian Rhythm in Non-small Cell Lung Cancer

https://doi.org/10.21203/rs.3.rs-39450/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted 07 Jul, 2020

You are reading this latest preprint version

Background: Circadian rhythm plays an influential role in the vast majority of behavior and physiology, even of tumorigenesis and development. In this study, we tried to explore the prognostic implication for non-small cell lung cancer (NSCLC) based on the expression profiles of circadian clock-related genes (CCRGs), and describe the changes of immune infiltration and cell functions of related to the circadian rhythm.

Methods: Univariate and multivariate Cox proportional hazard regression were performed to determine and construct a CCRGs risk-score significantly correlated with overall survival (OS) of the training set and validation set. GO, KEGG, and GSVA indicated discrepant changes in cellular processes and signaling pathways associated with these CCRGs. Further, immune cell infiltration, mutation rates, and related drugs and pathways analyses were investigated by the online analysis platform and the algorithm provided by works of literature.

Results: The risk-score based on ten-gene signatures could independently predict the OS both in TCGA lung adenocarcinoma (p < 0.0001, HR: 2.117, 95% CI: 1.546 to 2.900) and lung squamous cell carcinoma (p < 0.0001, HR: 2.066, 95% CI: 1.552 to 2.751), respectively. The risk-score also had superior accuracy and predictability (LUAD: AUC 0.788; LUSC: AUC 0.738). The prognostic performance of the ten-gene signature was finally validated in validation set based on GEO cohorts. The circadian oscillations driven by CCRGs could disturb the metabolism and cellular functions of cancer cells. We also found that the infiltration level of critical cells in specific anti-tumor immunity process was suppressed apparently. In contrast, the infiltrating of inflammatory cells and immune cells with negative regulatory effects were promoted in the high-risk group. CCRGs were evolutionarily conserved with low mutation rates, which brought difficulties to explore therapeutic targets.

Conclusions: We described the landscapes and prognostic implications of CCRGs. The risk-score based on CCRGs was an independent predictor of prognosis and could significantly stratify patient outcomes and immune cell infiltration levels. CCRGs were evolutionarily conserved with low mutation rates, which brought difficulties to explore therapeutic targets.Our study revealed that circadian rhythms might play an influential role in the NSCLC.

Translational Medicine

circadian rhythm

non-small cell lung cancer

risk-score

prognosis

cellular function

immune infiltration

Circadian rhythm exists endogenously in almost all organisms. The expression of circadian clock genes drives oscillatory changes in innumerable behavioral and physiological processes, including tumorigenesis and development(1, 2). According to the organism and cell type, the circadian clock promotes the rhythmic expression of 1% to over 60% of the genome, serving as the molecular basis for rhythmic control at the system's level(3). In recent years, it is becoming an increasing focus of the role of the circadian clock in tumorigenesis, cancer hallmarks, therapeutic options, and discussions of how circadian clock genes can lead a new dimension in future medicine(4, 5). Accumulating evidence identified that there was a tight association between cancer and disruption of circadian in curative effect and prognosis, and the core circadian transcripts are generally altered in many kinds of cancers(6–8). Nevertheless, the regulatory mechanism of circadian clock genes and their effects on clinical prognosis is not precise yet.

Nowadays, lung cancer is the most leading cause of cancer death (18.4% of the total cancer deaths) with the highest incidence (11.6% of the total cases) around the world(9). Non-small cell lung cancer (NSCLC) is the most common type and accounts for about 85% of total cases(10). Previous studies have shown that the dysregulation of the circadian rhythms can influence cancer development by regulating tumor cell apoptosis, autophagy, immune infiltration, and tumor cell-host interactions via related genes oscillatory and differential expression, which has been studied in many kinds of cancers including breast, colorectal cancer, and head and neck squamous cell carcinoma(11–13). However, landscapes and implications of differentially expressed circadian clock-related genes (CCRGs) in lung cancer have remained poorly defined. A potential role of these genes differential expression cannot be ignored between cancer cells and normal cells. In this present study, we purposed to explore a risk-score as a "Classifier" to predict the prognosis of patients based on the genomic expression profiles from public databases. Comprehensively, a total of 1,382 CCRGs with oscillatory transcripts with experimentally validated by techniques including RT-PCR, Northern blot, in situ hybridization, and Microarray or RNA-seq were analyzed in the present study(14). Among them, we identified 290 and 447 differentially expressed CCRGs in lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC) tissues, respectively based on the Cancer Genome of Atlas (TCGA) cohort. Then, a combination of univariate and Cox regression hazard regression analysis was used to screen the differential expression of CCRGs, which associated with overall survival. Subsequently, we established two optimal risk-score models to divide LUAD and LUSC patients into the high or low risk-score group, followed by verification in combined Gene Expression Omnibus (GEO) validation sets, respectively.

Further, we assessed the prognosis value and complementary value of molecular and clinical characteristics by survival, receiver operating characteristic (ROC) curve, and correlation analysis. We also identified the differences in the critical signaling pathways among these differential expression CCRGs using Gene Ontology (GO), the Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Variation Analysis (GSVA) methods. Finally, the correlative immune infiltrates in LUAD and LUSC, genetic alteration of CCRGs in risk-score models, and related drugs and signal pathways in cancer cells were also explored to provide novel ideas for the clinical translation of circadian genes.

Patient information and databases in training and validation sets

Thoroughly, a list of 1,382 homo sapiens (human) CCRGs validated by experiments including RT-PCR, Northern blot, and in situ hybridization were obtained from the Circadian Gene Database (The CGDB: http://cgdb.biocuckoo.org/index.php/)(14). In the training set, both the gene expression profiles (HTSeq - FPKM) and patient clinical information of 594 LUAD samples (normal count: 59; tumor count: 535) and 551 LUSC samples (normal count: 49; tumor count: 502) were downloaded from the TCGA database (https://portal.gdc.cancer.gov/). Patients who lacked follow-up information were exclude in the survival analysis. In the testing set, Microarray expression profiles and clinical information were obtained from the GEO database (ncbi.nlm.nih.gov/geo/) using the accession number GSE30219, GSE31210, GES3141, GSE37745, GSE50081, GSE68465, which contained more than a thousand samples of patients with lung cancer. Then, all samples were classified into the LUAD or LUSC type due to histological criteria to verify the risk-score model, respectively.

Data processing

To avoid the heterogeneity among different datasets and ensure a unified standard, the RNA-seq profiles were transformed using the formula log2(x + 1) and normalized. R version 3.6.1 (https://www.r-project.org/) software was used to normalize and process the data. Also, all data processing, analysis, and mapping were done using the R version 3.6.1 software and the Perl Programming Language version 5.28.1 (https://www.perl.org/) in the present study.

Functional enrichment analysis

To explore the pathways and interactions that are affected among these differential expression CCRGs, the statistical and visualize analysis of functional annotation (GO), including biological process, cellular component, and molecular function, and the KEGG pathway enrichment analysis and visualization were performed by using the R package "clusterProfiler" (http://www.bioconductor.org/). Finally, we also introduce the GSVA, a gene set enrichment method, to estimate the variation of pathway activity between the high risk-score and low risk-score groups in an unsupervised manner(15).

Risk-score model construction

We performed the univariate analysis and Cox proportional hazard regression to conducted to screen the differential expression CCRGs significantly associated with prognosis in the training cohort. Then, a risk score for each patient of prognostic risk was calculated respectively in LUAD and LUSC sets, according to the regression coefficients of the individual CCRGs screened from the multivariate Cox regression model and the expression value of each of the selected CCRGs. The computational formula used for this analysis was risk-score = h₀(t)*exp(β₁X₁ + β₁X₁+…+β_nX_n). β refers to the regression coefficient, and a hazard ratio (HR) value can be obtained after taking the natural logarithm exp (β). H₀(t) is the function of a benchmark risk, and h (t, X) is a risk function associated with X (covariant quantity) at time t. X (covariant quantity) represented the relative expression profiles of every CCRGs, which standardized by z-score. After modeling by multivariate Cox regression, the value of the risk score calculated by function "predict ()" is h (t, X). Then, all patients were into a high or low-risk group according to the median value as a cutoff to separately dichotomize the training sets in LUAD and LUSC, and a low-risk score indicates a superior prognosis for the patients. Finally, the risk-score signature of prognosis was verified in the testing sets, which integrated from the GEO cohorts (GSE30219, GSE31210, GES3141, GSE37745, GSE50081, GSE68465).

Immune infiltration analysis

We analyzed the correlation of these CCRGs included in the risk-score model with the abundance of immune infiltrates, including B cells, CD4 + T cells, CD8 + T cells, neutrophils, macrophages, and dendritic cells, via the TIMER (https://cistrome.shinyapps.io/timer/), in LUAD and LUSC. TIMER, including 10,897 samples across 32 cancer types from TCGA, is a comprehensive resource for systematic analysis of immune infiltrates across diverse cancer types, which depends on a deconvolution method to infer the abundance of tumor-infiltrating immune cells from RNA-Seq expression profiles(16). Gene expression levels against tumor purity are displayed on the left-most panel. All the genes expression levels were displayed with log2 RSEM.

Further, we also compared the difference of immune cell infiltration between the high-risk group and low-risk group by using the CIBERSORTx (https://cibersortx.stanford.edu/). CIBERSORTx is an extension of the CIBERSORT, which provides an analytical method to infer cell-type-specific gene expression profiles without digital cytometry(17). By the transcriptome profiling of single cells or sorted cell subpopulations based on a machine learning method, CIBERSORTx provides new possibilities for applying the signature matrix to bulk tissue expression profiles to infer cell-type proportions and represent cell type expression signatures. Based on this, we provided an analytical estimation of the abundances and distribution difference of 22 immune cell types in a mixed cell population of high-risk and low-risk group samples, using CCGRs expression data in LUAD and LUSC respectively. During the analysis, gene expression was corrected by the normalization using the R package "limma," and samples after the abundance estimation were filtered with the p-valve (< 0.05). When within-group data were merged, we combined the data by global averages.

Mutation analysis of CCRGs in risk-score model

The cBioPortal for Cancer Genomics (http://cbioportal.org) provides a visualized Web resource for exploring and analyzing multidimensional cancer genomics data, which is now developed and maintained by a multi-institutional team, such as the Memorial Sloan Kettering Cancer Center, the Dana Farber Cancer Institute, Princess Margaret Cancer Centre of Toronto, Children's Hospital of Philadelphia. In this portal, molecular profiling data of cancer tissues and cell lines were reduced into readily visual and understandable genetic, epigenetic, gene expression, proteomic events, and clinically relevant events. According to the cBioPortal, we investigate mutations and expressions of CCRGs in the risk-score model in lung adenocarcinoma (TCGA PanCancer Atlas) (566 samples) and lung squamous Cell Carcinoma (TCGA PanCancer Atlas) (487 samples) respectively.

Related drugs and pathways analysis

We also analyzed related anti-tumor drugs and related pathways of the CCRGs included in the risk-score model. All relevant information was retrieved by searching the GeneCards (The Human Gene Database) (https://cistrome.shinyapps.io/timer/) and supplemented by retrieving relevant literature in the PubMed (https://pubmed.ncbi.nlm.nih.gov/). The results are presented in the form of the Sankey diagram.

Statistical analysis

All statistical analyzes in the present study were performed using the R version 3.6.1 software (https://www.r-project.org/), and p-value < 0.05 was regarded as statistically significant for all the analyses. The Kruskal-Wallis test and one-way ANOVA were used to check the expression differences among the genes and the association of the risk-score with clinical signatures. OS was analyzed by the Kaplan-Meier survival curve and the log-rank test to check the significant difference between the high-risk and low-risk groups. The univariate analysis and multivariate Cox proportional hazards regression model were used to analyze the key CCRGs that affect the prognosis of NSCLC patients. The ROC curve analysis was used to evaluate the sensitivity and specificity of prognostic prediction of the CCRGs signature risk-score model. The prognostic accuracy was presented by the area under the ROC curve (AUC). All tests were two-sided.

Identification and screening of differential expression genes

We analyzed the differential expression profile of a total of 1,382 homo sapiens (human) CCRGs in 594 LUAD samples (normal count: 59; tumor count: 535) and 551 LUSC samples (normal count: 49; tumor count: 502) by the Wilcoxon signed-rank test respectively. Finally, we identified 290 CCRGs (downregulated CCRGs: 126; upregulated CCRGs: 164), and 447 CCRGs (downregulated CCRGs: 226; upregulated CCRGs: 221) differentially expressed in LUAD and LUSC samples respectively. Results of expression analyses are illustrated as heatmaps, volcano plots, and boxplots (Fig. 1).

Functional enrichment analysis of differential expression genes

Firstly, we analyzed the association of these differentially expressed CCRGs with the GO terms of the biological process (BP) and cellular component (CC) categories. For LUAD samples, the top five enriched BP terms were 'regulation of hemopoiesis,’ 'neutrophil degranulation,’ 'neutrophil activation involved in immune response,’ 'neutrophil-mediated immunity,' and 'neutrophil activation (Fig. 2A-D).' The top two enriched CC terms were 'chromatin' and 'secretory granule membrane (Fig. 2A-D).' For LUSC samples, the top five enriched BP terms were 'neutrophil activation involved in immune response,’ 'neutrophil activation,’ 'neutrophil degranulation,’ 'neutrophil-mediated immunity,' and 'negative regulation of immune system process,' which was similar to the results of the former, while more inclined towards the immunomodulation (Fig. 2E-H). For the results of CC, the top two enriched terms were precisely the same as the former (Fig. 2E-H). Then, results of KEGG analysis indicated that altered CCRGs were mainly involved in the systemic lupus erythematosus and osteoclast differentiation of LUAD samples, meanwhile in the systemic lupus erythematosus, apoptosis, and B cell receptor signaling pathway of LUSC samples (Fig. 3).

Finally, we further investigated the differential distribution of signal pathway enrichment of differentially expressed CCRGs set between the high risk-score and low risk-score groups using the GSVA method. In LUAD samples, compared with the low-risk group, CCRGs of the high-risk group mainly enriched in the G2M_CHECKPOINT, COAGULATION, MYC_TARGETS_V2, and MITOTIC_SPINDLE. Meanwhile, CCRGs of the low-risk group mainly enriched in the MYOGENESIS and NOTCH_SIGNALING. However, the results of the LUSC samples are quite different. CCRGs of the high-risk group mainly enriched in the UV_RESPONSE_DN, KRAS_SIGNALING_UP, TNFA_SIGNALING_VIA_NFKB, INFLAMMATORY_RESPONSE, and INTERFERON_GAMMA_RESPONSE, while CCRGs were mainly involved in the E2F_TARGETS, UNFOLDED_PROTEIN_RESPONSE, MYC_TARGETS_V1, G2M_CHECKPOINT, and DNA_REPAIR in the low-risk group (Fig. 4). It indicated the differential expression and regulation of CCRGs in LUAD and LUSC. This difference may lead to different biological processes related to the circadian rhythm in cancer cells.

Construction of the prognostic risk-score model in the training set

To explore the prognostic signatures of these differentially expressed CCRGs, we preliminarily performed a univariate analysis of the standardized expression of the 290 CCRGs of LUAD and 447 CCRGs of LUSC in the training sets to identify the prognostic CCRGs respectively. The results showed that the expression of 31 CCRGs and 70 CCRGs each significantly correlated with the OS (p < 0.05) of LUAD patients and LUSC, respectively (Fig. 5A-B). Further, we performed Cox proportional hazard regression to screen the ultimate CCRGs of risk-score models. Subsequently, results indicated that the expression of 10 CCRGs correlated with the OS of LUAD and LUSC patients, respectively (Table 1). According to this, we constructed a risk-score model for predicting the prognosis of LUAD and LUSC patients using the calculation formula mentioned in the method part, respectively. Finally, CDA, POU2AF1, TUBB6, SPAG8, NT5E, ARRB1, DDIT4, HAL, PHLDB2, and AGMAT as risk genes in the risk-score model of LUAD and ALOX5AP, RALGAPA2, TIGD3, PNPLA6, ALPL, TREM1, VSIG4, CD300C, HIST1H2BH, and WNT10A as risk genes in the LUSC risk-score model (Fig. 5C-F). Survival analysis revealed that there was a significant difference between the high-risk group and the low-risk group in OS, and patients in the high-risk group significantly correlated with an inferior prognosis (LUAD: p < 0.0001, HR: 2.117, 95% CI: 1.546 to 2.900; LUSC: p < 0.0001, HR: 2.066, 95% CI: 1.552 to 2.751) (Fig. 6A-B). Finally, we also ranked the risk scores of LUAD and LUSC patients for OS and explored the distribution features (Fig. 6C-D). The dot plots revealed the status of each patient in the training sets (Fig. 6E-F). The heat maps showed the differential expression of the feature CCRGs in the high-risk and low-risk groups (Fig. 6G-H). As results show, there was an upregulation of HAL, PHLDB2, AGMAT, DDIT4, CDA, NT5E, and TUBB6 as high-risk genes and downregulation of POU2AF1, SPAG8, and ARRB1 as protective genes of LUAD patients in the high-risk score group comparing the low-risk score group. Moreover, samples with high-risk scores of LUSC patients suggested upregulation of TREM1, WNT10A, CD300C, ALPL, ALOX5AP, VSIG4, RALGAPA2 and PNPLA6 as high-risk genes and downregulation of TIGD3 and HIST1H2BH as protective genes.

Table 1

Descriptions of CCRGs of the risk-score model in Cox proportional hazard regression analysis.
Type	Gene ID	Full name or Description	Risk coefficient	HR	HR.95L	HR.95H	p-value
LUAD	CDA	Cytidine Deaminase	-0.089059134	0.914791	0.822109	1.017923	0.102253
	POU2AF1	POU Class 2 Associating Factor 1	-0.295832916	0.743912	0.631686	0.876075	0.000392
	TUBB6	Tubulin Beta 6 Class V	0.253501154	1.288529	1.031295	1.609925	0.025669
	SPAG8	Sperm Associated Antigen 8	-0.215285427	0.806311	0.602534	1.079006	0.14751
	NT5E	5'-Nucleotidase Ecto	0.165120683	1.179535	1.040611	1.337006	0.009806
	ARRB1	Arrestin Beta 1	-0.229636783	0.794822	0.638692	0.989119	0.039588
	DDIT4	DNA Damage Inducible Transcript 4	0.128789666	1.137451	0.974568	1.327557	0.102412
	HAL	Histidine Ammonia-Lyase	0.254600207	1.289946	1.101163	1.511093	0.001613
	PHLDB2	Pleckstrin Homology Like Domain Family B Member 2	0.232732475	1.262044	1.006844	1.581927	0.043472
	AGMAT	Agmatinase	0.276379132	1.318348	1.00177	1.73497	0.048543
LUSC	ALOX5AP	Arachidonate 5-Lipoxygenase Activating Protein	-0.204079367	0.815398	0.638768	1.040868	0.10134
	RALGAPA2	Ral GTPase Activating Protein Catalytic Subunit Alpha 2	0.209187368	1.232676	0.963034	1.577816	0.096734
	TIGD3	Tigger Transposable Element Derived 3	-0.338879739	0.712568	0.445497	1.139747	0.157327
	PNPLA6	Patatin Like Phospholipase Domain Containing 6	0.262619794	1.300332	0.94122	1.79646	0.11125
	ALPL	Alkaline Phosphatase, Biomineralization Associated	0.088447658	1.092477	0.977686	1.220746	0.118395
	TREM1	Triggering Receptor Expressed on Myeloid Cells 1	0.205116552	1.227668	1.028907	1.464826	0.022834
	VSIG4	V-Set and Immunoglobulin Domain Containing 4	0.21911239	1.244971	1.000762	1.548772	0.049206
	CD300C	CD300 Antigen-Like Family Member C	-0.355141335	0.701074	0.440876	1.114837	0.133452
	HIST1H2BH	histone cluster 1, H2bh	-0.161357251	0.850988	0.749454	0.966278	0.012805
	WNT10A	Wnt Family Member 10A	0.127271005	1.135725	1.008189	1.279394	0.036246

Validation of the prognostic risk-score model in the testing set

We next validated the stability and accuracy of the prognostic risk-score model in the testing sets, which included LUAD and LUSC cohorts from the GEO database. The OS was selected as the key indicator to compared the groups and samples were divided into low and high risk-score groups based on the calculated risk score. The formula is as mentioned before. For the testing set of LUAD type, 519 and 544 samples were separated into low and high risk-score groups, respectively. Survival analysis showed that there was a significant difference between the high and low risk-score groups (p < 0.0001, HR: 1.493, 95% CI: 1.248 to 1.787) (Fig. 7A). Similarly, 88 samples of the low-risk group and 89 samples of the high-risk group were included in the survival analysis of the LUSC testing set. Results also suggested a significant difference between the high and low risk-score groups (p = 0.0486, HR: 1.453, 95% CI: 1.002 to 2.105) (Fig. 7B). To summarize, our results confirmed that these two risk-score models based on CCRGs signatures were all stable and accurate in predicting the prognosis of patients.

Clinical characteristics correlation analysis

In this section, we further explored the stability and reliability of the risk score as a clinical indicator. Seven variables, including age, gender, T stage, N stage, M stage, pathologic stage, and risk-score, were analyzed using the univariate analysis and Cox proportional hazard regression in the training sets of LUAD and LUSC respectively (LUAD: p < 0.001, HR: 1.229, 95% CI: 1.157 to 1.305; LUSC: p < 0.001, HR: 2.609, 95% CI: 2.073 to 3.284) (Fig. 8A-D). Our analysis showed that risk-score was found to be an independent prognostic indicator both in LUAD patients and LUSC patients. Then, we constructed ROC curves for different variables to evaluate the risk-score as classifiers, and the AUC was calculated and considered as the basis for evaluation (LUAD: AUC 0.788; LUSC: AUC 0.738) (Fig. 8E-F). Our results indicated that the risk-score had superior accuracy and predictability comparing other clinical characteristics both in LUAD and LUSC samples.

Immune correlation and infiltration analysis

Firstly, correlation analysis between immune cell infiltration levels and CCRGs of the risk-score model was performed using the TIMER. We initially assessed the correlations of characteristic CCRGs expression with immune infiltration levels in LUAD and LUSC samples using the TIMER one by one. We found that almost all CCRGs expression in the risk-score model has significant correlations with infiltrating levels of B cells, CD8 + T cells, CD4 + T cells, macrophages, neutrophils, and dendritic cells both in LUAD and LUSC (Supplementary Fig. 1–2). Further, we explored the immune infiltration level of the diverse immune infiltrating cells between the high-risk group and the low-risk group by the CIBERSORTx. We first calculated the abundances of 22 immune cell types in each sample using standardized CCGRs expression profiles in LUAD and LUSC, respectively. Then we filtered each sample according to the p-value to eliminate the bias caused by inaccurate estimation and grouped the samples based on the risk-score. Finally, 23 samples of the low-risk group and ten samples of the high-risk group were selected in LUAD. Meanwhile, of LUSC samples, 19 samples of the low-risk group, and 26 samples of the high-risk group were included for analysis. Our results suggested that immune infiltration levels in the high-risk group of naive B cells, plasma cells, naive CD4 + T cells, CD8 + T cells, dendritic cells, M2 macrophages, and mast cells decreased significantly in the LUAD compared with the low-risk group (Fig. 9A-B). Furthermore, there were significantly increased levels of infiltrates of activated memory CD4 + T cells, T follicular helper cells, regulatory T cells (Tregs), natural killer (NK) cells, monocytes, neutrophils, and M0/1 macrophage in the high-risk group. Moreover, for LUAC samples, immune infiltration levels in the high-risk group of naive B cells, resting memory CD4 + T cells, naive CD4 + T cells, CD8 + T cells, T cells follicular helper, NK cells, dendritic cells, and M0/1 macrophage decreased significantly comparing the low-risk group (Fig. 9C-D). Meanwhile, there were significantly increased levels of infiltrates of plasma cells, activated memory CD4 + T cells, Tregs, mast cells, monocytes, M2 macrophages, and neutrophils in the high-risk group. In summary, we found that there were similar commonalities in the level of immune cell infiltrations between the high-risk group and the low-risk group in LUAD and LUSC. For instance, our results showed that immune cells, which play a vital role in tumor immunity, were significantly decreased in the high-risk group, such as CD8 + T cells, dendritic cells, B cells, and naive CD4 + T cells. Meanwhile, negative regulation and inflammation-related cells, such as Tregs, mast cells, monocytes, and neutrophils, were significantly increased. These results suggested that compared with the low-risk group, anti-tumor immune responses were inhibited while the inflammatory process and hyperresponsiveness may be significantly enhanced in the high-risk group. However, there were still some differences and changes in the level of infiltration of some immune cells between the LUAD and LUSC, which might illustrate the difference between LUAD and LUSC at the level of circadian-rhythm-related immune infiltration.

Genetic alteration analysis

We investigated the genetic alteration of these CCRGs in the risk-score model to further understand their contributions to carcinogenesis by the cbioportal. We found that these risk-associated CCRGs were relatively conservative, and the mutation rates of these genes were all lower than 3% (most the percent of which were around 1%-2%) both in LUAD and LUSC (Fig. 10). The low frequent genetic alterations suggested the high stability and conservatism of these CCRGs and the crucial roles of these genes in the genetic, epigenetic, and development of lung cancer.

Related drugs and related pathways

We analyzed related drugs and signaling pathways related to these CCRGs of the risk-score model in tumor cell lines. All relevant information was retrieved by searching the GeneCards and supplemented by retrieving relevant literature in PubMed. We found that although some anti-tumor drugs targeting these CCRGs have been applied in clinical practice or were in clinical trials, there were still many blank spots in relevant research, especially for LUSC patients (Fig. 11). As the results showed, most of the existing drugs played a corresponding role in cancer through classic signaling pathways, such as the EGFR signaling pathway, NF-kB signaling pathway, and PI3K-Akt signaling pathway.

Circadian biology is being considered as the fourth dimension of modern medicine(18). There is clear evidence for the differential expression of CCRGs in a variety of diseases, and cancers are no exception(19). At the same time, these self-perpetuating differential expressions were also synchronized by the external environment (light, temperature, ingestion, and activities) that were known as the zeitgebers. Also, the process of differential expression was 24-hour periodicity and might also be affected by the seasons. Universally recognized, implications of regulating CCRGs expression in epigenetic control mechanisms have been described during the tumor initiation and progression, which included circadian metabolic changes and tumor-derived macroenvironment, which has been reported in studies of breast cancer and LUAD in mouse models(20). The implication also has been indicated by epidemiological studies(21). Changes and disruptions of circadian rhythms in humans significantly impinged on the increasing risk of tumorigenesis(22, 23). Evidently, circadian biology is becoming a critical involvement in improving the understanding of molecular mechanisms involved in cancer cells. Nevertheless, its importance has sparsely been well recognized in clinical studies and practice, and even more when translating to the bedside. Based on this, we attempted to represents a substantial step toward that direction, which aims first to describe landscapes and implications of these differentially expressed CCRGs and investigate the connection between impingement of circadian rhythms and prognostic significance, in the most common and malignant tumor. We finally integrated and analyzed the expression profiles of 1,382 human CCRGs in NSCLC wholly and systematically via the CGDB, TCGA, and GEO database. We pioneering proposed a CCGRs-based risk-score model better to assess the effects of circadian rhythm on prognosis accordingly. In addition, according to the score, we further focused on the infiltration changes of immune cells, genetic alteration, and the possibility of being a pharmacological target in these samples.

Recently, a study based on integration and analysis of data from the TCGA database has investigated the association between 14 kinds of clock genes and prognostic signatures in NSCLC patients, which also showed that differentially expressed clock genes constitute their characteristic asynchronous circadian rhythms(24). To date, thousands of genes and proteins are considered to be related to the circadian rhythms' oscillation. Given the significance of circadian rhythms in lung cancer, it is reasonable to speculate that CCRGs hold excellent promises in prognostic prediction and that a risk score based on multiple-gene signatures derived from dependable algorithms would be more reliable and superior to any single molecules in predicting prognosis of NSCLC. We, therefore, put forward a risk-score model, in which ten-gene signatures were selected and calculated for evaluating the prognostic risk of LUAD and LUSC training sets, and the predictive validity of the risk-score model was validated in several GEO NSCLC cohorts, respectively. Fortunately, the risk scores significantly stratified patient outcomes and immune cell infiltration levels between the high-risk and low-risk groups. Further, the risk-score also showed its excellent stability and accuracy as a classifier in the Cox proportional hazard regression, including risk-score and other clinical variables. It was evident that the high-risk group has an inferior prognosis and a more reduced anti-tumor immune response in our analysis. In the risk-score model containing ten genes of LUAD and LUSC respectively, we found that these genes are involved in many biological processes, such as cell cycle control, metabolism, immune-modulating, inflammatory reaction, cytoskeletal reorganization, chromatin remodeling, apoptosis in response to DNA damage repair, and protein synthesis and transportation, through systemic functional analysis. We concluded that these above processes in tumor cells might be affected by the circadian rhythm. Several studies have demonstrated that a wide range of core circadian clock components is epigenetically altered, and this perturbation could promote tumorigenesis, progression, and decreased survival in lung cancer, which also suggested an essential position of circadian homeostasis in the tumor-suppressive role(25, 26).

Interestingly, we found that the infiltration level of critical cells in specific anti-tumor immunity process, such as CD4 + T cells, CD8 + T cells, and dendritic cells, were suppressed apparently, while the activity and infiltrating of inflammatory cells and Tregs with negative regulatory ability were promoted in the high-risk group. It established that circadian rhythms and related genes played a vital role in the tumor immune and tumor-associated inflammatory response. The latest studies have confirmed our results. To date, current notion suggests that CCRGs express in most immune cells universally and present a circadian oscillation with a fixed rhythm, which performs essential roles in a wide range of immunomodulation process, including the phagocytosis, apoptosis, the synthesis, and release of cytokines, chemokines, and cytolytic factors, the response occurring through pattern recognition receptors(27). Differential expression of CCRGs also plays a vital role in the development and specification of immune cell lineages(28). This view also reflected in our analysis. For instance, immune infiltration level of resting memory CD4 + T cells and naive CD4 + T cells were decreased, while the level of activated memory CD4 + T cells was increased in the high-risk group. Consequently, it is evidence that alterations in circadian rhythms due to differential expression of genes in cancer cells may lead to disturbed the immune responses, and these changes may be caused by clock gene mutation, environmental disruption, or the age and tumor itself. A study of circadian rhythm reprogramming during the lung inflammation suggested that the early events in lung injury may produce a complex reorganization of cellular and molecular circadian rhythms and further regulate immune responses of the host(29). It will be essential to determine the mechanism and causality of oscillations driven by CCRGs in cellular function, metabolism and immunity, and whether the critical drivers for oscillations are the time of day/season-dependence. If so, it might strengthen our fundamental understanding of how the circadian rhythm disturbs metabolism and immune functions to anticipate changes in the environment, and provide a bridge between the circadian rhythms and novel insights to facilitate the development of chronotherapies for fighting cancer and other diseases.

Besides, our genetic alteration analysis also suggested the low mutation rates of these CCRGs in the risk-score model, which was also in line with the current view that CCRGs were evolutionarily conserved in eukaryotes(18). Moreover, this conservatism would affect plenty of critical cell functions, such as immunomodulatory. Studies over the last decade indicate that immune responses related to the circadian oscillators are a consequence of this Darwinian selection process, and the circadian rhythm could minimize costs and maximizes benefits of immunity to optimize organismal fitness in a given environment(30). Thus, the disruption of the normal circadian rhythmic may result in the appearance of CCRGs differential expression and metabolic rhythms, which might function to support host immunity but also increase the probability of tissue damage and a catastrophic vulnerability(31). Furthermore, because of this, in the analysis of related anti-tumor drugs of these CCRGs in the risk-score model, there are few direct targets of these genes. Meanwhile, we also found the subtle difference between CCRGs differential expression and immune cell infiltration in LUAD and LUSC, which might result from the specific contexts of different types of cancers. These are still urgent questions needed to be studied and solved today.

All organisms on Earth are exposed to regular environmental cycles generated by the rotation and revolution of the Earth. This, in turn, has led to the evolution of circadian rhythms driven by CCRGs, which facilitate lives to anticipate and adapt to the internal and external changes during their environment. We preliminary explored a risk-score based on ten CCRGs signatures based on TCGA and GEO database in LUAD and LUSC, respectively. This risk-score was an independent predictor of prognosis. Further analysis of cell functions and immune infiltration between the high-risk and low-risk group and genetic alteration of these GGRGs also investigated in our study. Differential expression of CCRGs also regulated the immune cell infiltration level in NSCLC. These CCRGs were evolutionarily conserved with low mutation rates. Further studies and experimental confirmations are needed to explore potential drug targets for therapy.

NSCLC

non-small cell lung cancer; CCRGs:circadian clock-related genes; LUAD:lung adenocarcinoma; LUSC:lung squamous cell carcinoma; TCGA:The Cancer Genome of Atlas; GEO:Gene Expression Omnibus; GO:Gene Ontology; KEGG:Kyoto Encyclopedia of Genes and Genomes; GSVA:Gene Set Variation Analysis; ROC:receiver operating characteristic; CGDB:Circadian Gene Database; HR:hazard ratio; AUC:area under the ROC curve; BP:biological process; CC:cellular component; Tregs:regulatory T cells; NK cells:natural killer (NK) cells;

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable

Availability of data and materials

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

Competing interests

There are no conflicts of interest to declare.

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 81672268) (www.nsfc.gov.cn/).

Authors' contributions

(I) Conception and design: Ming Li, MD; Tian Jiang, Ph.D.; Cheng Zhan, Ph.D.

(II) Administrative support: Junjie Xi, MD; Cheng Zhan, Ph.D.; Mingxiang Feng, MD

(III) Provision of study materials or patients: Ming Li, MD; Zhencong Chen, MD; Tian Jiang, Ph.D.

(IV) Collection and assembly of data: Ming Li, MD; Zhencong Chen, MD;Xiaodong Yang, MD; Yajing Du, MS; Jiaqi Liang, MS, Lin Wang, Ph.D.

(V) Data analysis and interpretation: Ming Li, MD, Tian Jiang, Ph.D.; Zhencong Chen, MD

(VI) Manuscript writing: Ming Li, MD.

(VII) Final approval of manuscript: All authors.

Acknowledgments

None.

Sulli G, Lam M, Panda S. Interplay between Circadian Clock and Cancer: New Frontiers for Cancer Treatment. Trends Cancer. [Journal Article; Review]. 2019 2019-08-01;5(8):475–94.
Xu C, Ochi H, Fukuda T, Sato S, Sunamura S, Takarada T, et al. Circadian Clock Regulates Bone Resorption in Mice. J BONE MINER RES. [Journal Article; Research Support, Non-U.S. Gov't]. 2016 2016-07-01;31(7):1344–55.
Li S, Zhang L. Circadian Control of Global Transcription. BIOMED RES INT. [Journal Article; Review]. 2015 2015-01-20;2015:187809.
Sulli G, Manoogian E, Taub PR, Panda S. Training the Circadian Clock, Clocking the Drugs, and Drugging the Clock to Prevent, Manage, and Treat Chronic Diseases. TRENDS PHARMACOL SCI. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't; Research Support, U.S. Gov't, Non-P.H.S.; Review]. 2018 2018-09-01;39(9):812 – 27.
Ercolani L, Ferrari A, De Mei C, Parodi C, Wade M, Grimaldi B. Circadian clock: Time for novel anticancer strategies? PHARMACOL RES. [Journal Article; Review]. 2015 2015-10-01;100:288 – 95.
Shostak A. Circadian, Clock. Cell Division, and Cancer: From Molecules to Organism. INT J MOL SCI. [Journal Article; Review]. 2017 2017-04-20;18(4).
Kelleher FC, Rao A, Maguire A. Circadian molecular clocks and cancer. CANCER LETT. [Journal Article; Review]. 2014 2014-01-01;342(1):9–18.
J CLIN ONCOL. [Journal Article

1993 1993-07-01

Hrushesky WJ, Bjarnason GA. Circadian cancer therapy. J CLIN ONCOL. [Journal Article; Research Support US, Gov't Non-PHS, Research Support US, Gov't PHS, Review]. 1993 1993-07-01;11(7):1403–17.
Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin [Journal Article]. 2018;68(6):394–424. 2018-11-01.
Zakowski MF. Cytology nomenclature and 2015 World Health Organization classification of lung cancer. CANCER CYTOPATHOL. [Journal Article; Review]. 2016 2016-02-01;124(2):81–8.
Alhopuro P, Bjorklund M, Sammalkorpi H, Turunen M, Tuupanen S, Bistrom M, et al. Mutations in the circadian gene CLOCK in colorectal cancer. MOL CANCER RES. [Journal Article; Research Support, Non-U.S. Gov't]. 2010 2010-07-01;8(7):952–60.
Blakeman V, Williams JL, Meng QJ, Streuli CH. Circadian clocks and breast cancer. BREAST CANCER RES. [Journal Article; Review; Research Support, Non-U.S. Gov't]. 2016 2016-09-02;18(1):89.
Rahman S, Kraljevic PS, Markova-Car E. Circadian (De)regulation in Head and Neck Squamous Cell Carcinoma. INT J MOL SCI. [Journal Article; Review]. 2019 2019-05-30;20(11).
Li S, Shui K, Zhang Y, Lv Y, Deng W, Ullah S, et al. CGDB: a database of circadian genes in eukaryotes. NUCLEIC ACIDS RES. [Journal Article; Research Support, Non-U.S. Gov't]. 2017 2017-01-04;45(D1):D397-403.
GSVA: gene set variation analysis for microarray and RNA-seq data

[Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't]. 2013 2013-01-16;14:7.

Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC BIOINFORMATICS. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't]. 2013 2013-01-16;14:7.
Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. CANCER RES. [Journal Article]. 2017 2017-11-01;77(21):e108-10.
Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. NAT BIOTECHNOL. [Journal Article; Research Support, N.I.H. Extramural; Research Support, Non-U.S. Gov't; Research Support, U.S. Gov't, Non-P.H.S.]. 2019 2019-07-01;37(7):773–82.
Cederroth CR, Albrecht U, Bass J, Brown SA, Dyhrfjeld-Johnsen J, Gachon F, et al. Medicine in the Fourth Dimension. CELL METAB. [Journal Article; Review]. 2019 2019-08-06;30(2):238–50.
Kelleher FC, Rao A, Maguire A. Circadian molecular clocks and cancer. CANCER LETT. [Journal Article; Review]. 2014 2014-01-01;342(1):9–18.
Masri S, Sassone-Corsi P. The emerging link between cancer, metabolism, and circadian rhythms. NAT MED. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't; Review]. 2018 2018-12-01;24(12):1795–803.
Masri S, Kinouchi K, Sassone-Corsi P. Circadian clocks, epigenetics, and cancer. CURR OPIN ONCOL. [Journal Article; Research Support, N.I.H. Extramural; Research Support, Non-U.S. Gov't; Review]. 2015 2015-01-01;27(1):50–6.
Kloog I, Haim A, Stevens RG, Portnov BA. Global co-distribution of light at night (LAN) and cancers of prostate, colon, and lung in men. CHRONOBIOL INT. [Journal Article]. 2009 2009-01-01;26(1):108–25.
Schernhammer ES, Kroenke CH, Laden F, Hankinson SE. Night work and risk of breast cancer. EPIDEMIOLOGY. [Journal Article]. 2006 2006-01-01;17(1):108–11.
Qiu M, Chen YB, Jin S, Fang XF, He XX, Xiong ZF, et al. Research on circadian clock genes in non-small-cell lung carcinoma. CHRONOBIOL INT. [Journal Article; Research Support, Non-U.S. Gov't]. 2019 2019-06-01;36(6):739–50.
Papagiannakopoulos T, Bauer MR, Davidson SM, Heimann M, Subbaraj L, Bhutkar A, et al. Circadian Rhythm Disruption Promotes Lung Tumorigenesis. CELL METAB. [Journal Article]. 2016 2016-08-09;24(2):324–31.
Kettner NM, Katchy CA, Fu L. Circadian gene variants in cancer. ANN MED. [Journal Article; Research Support, N.I.H., Extramural; Research Support, U.S. Gov't, Review] Non-PHS. 2014 2014-06-01;46(4):208–20.
Man K, Loudon A, Chawla A. Immunity around the clock. SCIENCE. [Journal Article; Research Support, N.I.H. Extramural; Research Support, Non-U.S. Gov't; Review]. 2016 2016-11-25;354(6315):999–1003.
Labrecque N, Cermakian N. Circadian Clocks in the Immune System. J Biol Rhythms. [Journal Article; Research Support, Non-U.S. Gov't; Review]. 2015 2015-08-01;30(4):277–90.
Haspel JA, Chettimada S, Shaik RS, Chu JH, Raby BA, Cernadas M, et al. Circadian rhythm reprogramming during lung inflammation. NAT COMMUN. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't]. 2014 2014-09-11;5:4753.
Stearns SC. Evolutionary medicine: its scope, interest and potential. Proc Biol Sci. [Journal Article; Review]. 2012 2012-11-07;279(1746):4305–21.
Haspel JA, Chettimada S, Shaik RS, Chu JH, Raby BA, Cernadas M, et al. Circadian rhythm reprogramming during lung inflammation. NAT COMMUN. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't]. 2014 2014-09-11;5:4753.

Supplementaryfigure2.tif
Supplementaryfigure2. Correlation of expression of CCRGs in the risk-score model with immune infiltration level in LUSC.
Supplementaryfigure1.tif
Supplementaryfigure 1.Correlation of expression of CCRGs in risk-score model with immune infiltration level in LUAD.

Download PDF

Version 1

posted 07 Jul, 2020

You are reading this latest preprint version

Genetic Landscapes and Prognostic Implications of Circadian Rhythm in Non-small Cell Lung Cancer

Status:

Version 1

Abstract

Figures

Background

Methods

Patient information and databases in training and validation sets

Data processing

Functional enrichment analysis

Risk-score model construction

Immune infiltration analysis

Mutation analysis of CCRGs in risk-score model

Related drugs and pathways analysis

Statistical analysis

Results

Identification and screening of differential expression genes

Functional enrichment analysis of differential expression genes

Construction of the prognostic risk-score model in the training set

Validation of the prognostic risk-score model in the testing set

Clinical characteristics correlation analysis

Immune correlation and infiltration analysis

Genetic alteration analysis

Related drugs and related pathways

Discussion

Conclusions

Abbreviations

Declarations

References

Supplementary Files

Status:

Version 1