Disruptions of Circadian Genes in Cutaneous Melanoma—An In Silico Analysis of Transcriptome Databases

Circadian genes are a set of genes that regulate the body’s internal clock and influence various physiological processes, including sleep–wake cycles, metabolism and immune function. Skin cutaneous melanoma (SKCM) is a type of skin cancer that arises from the pigment-producing cells in the skin and is the most deadly form of skin cancer. This study has investigated the relevance of circadian gene expression and immune infiltrations in the outcomes of cutaneous melanoma patients. In the present study, in silico methods based on the GEPIa, TIMER 2.0 and cBioPortal databases were performed, so as to investigate the transcript level and prognostic value of 24 circadian genes in SKCM and their relationship with the immune infiltration level. The in silico analysis showed that significantly more than half of the investigated circadian genes have an altered transcript pattern in cutaneous melanoma compared to normal skin. The mRNA levels of TIMELES and BHLHE41 were upregulated, whereas those of NFIL3, BMAL1, HLF, TEF, RORA, RORC, NR1D1, PER1, PER2, PER3, CRY2 and BHLHE40 were downregulated. The presented research shows that SKCM patients with at least one alteration of their circadian genes have decreased overall survival. Additionally, majority of the circadian genes are significantly corelated with the immune cells’ infiltration level. The strongest correlation was found for neutrophils and was followed by circadian genes: NR1D2 r = 0.52 p < 0.0001, BMAL1 r = 0.509 p < 0.0001; CLOCK r = 0.45 p < 0.0001; CSNKA1A1 r = 0.45 p < 0.0001; RORA r = 0.44 p < 0.0001. The infiltration level of immune cells in skin tumors has been associated with patient prognosis and treatment response. Circadian regulation of immune cell infiltration may further contribute to these prognostic and predictive markers. Examining the correlation between circadian rhythm and immune cell infiltration can provide valuable insights into disease progression and guide personalized treatment decisions.


Introduction
Cutaneous melanoma (SKCM) is a type of skin cancer arising from the melanocytes, which has demonstrated increasing incidence rates over the past few decades. It is considered to be one of the most aggressive forms of skin cancer and is responsible for the majority of skin cancer-related deaths [1]. The highest incidence rates of cutaneous melanoma have been reported in Australia, followed by the southern states of the United States [2]. Study of the circadian rhythm in cutaneous melanoma is an area of research that has gained considerable attention in recent years [3,4]. This is because disruptions in the circadian rhythm have been linked to an increased risk of cancer and it is believed that understanding the circadian rhythm in cutaneous melanoma may lead to the new prevention and treatment strategies [5].
Every organism needs to maintain homeostasis to function properly. Circadian rhythm is the basic process that helps maintain it. Over the course of evolution and adaptation to the external environment, all organisms on Earth have developed a mechanism that regulates the rhythms of all physiological processes [6]. Circadian rhythm is a fundamental process that organizes the functioning of an organism in response to external stimuli. This adaptation process is observed in the entire organism at a physiological level, as well as peripherally-at the levels of tissues and cells-by influencing various molecular pathways regulated by circadian genes transcriptional-translational feedback loops, such as metabolism, cell proliferation, inflammation and DNA damage repair [6]. Molecular clocks are controlled by multiple genes that form basal and auxiliary transcriptional and translational feedback loops with both positive and negative characteristics. To date, approximately 20 candidate genes related to the circadian rhythm generation and maintenance in the hippocampus and suprachiasmatic nucleus (SCN), as well as in peripheral tissues including skin, have been identified (Supplementary Materials Table S1). The circadian-driven feedback loop consists of transcription-inducing complexes of the negative clock regulators PER and CRY, along with other positive factors BMAL1/CLOCK or BMAL1/NPAS2. In the cytoplasm, PER and CRY are phosphorylated by CSNK1A1. Accumulated clock proteins translocate to the nucleus and inhibit its BMAL1 expression in the form of protein complexes. Other bHLH-PAS transcription factor family members-RORs (RAR-related orphan receptor), NR1Ds or REV-ERBs (nuclear receptor subfamily 1, group D members)-constitute the stabilizing loop by repressing BMAL1 transcription via retinoic acid-related orphan receptor response element (RORE) in the promoter region of BMAL1 [7,8]. The regulatory role (particularly of the BMAL1/CLOCK heterodimer) in controlling expression of various genes involved in multiple biological processes in a cell, such as genetic alterations in metabolism, inflammatory pathways, cell cycle arrest, DNA damage repair, cell proliferation, maintenance of genomic stability, oxidative stress and apoptosis, is a unique feature of circadian genes. Moreover, circadian genes also affect the expression of genes involved in skin barrier function, inflammation and response to UV radiation [9]. Perturbations in these processes are hallmarks of cancer including skin cancer, while chronic circadian rhythm disruption predisposes tumor development [7,8].
Inflammation is another factor that has been implicated in the development and progression of cutaneous melanoma. Chronic inflammation has been linked to an increased risk of cancer. It is also believed that inflammation plays a role in the growth and spread of melanoma cells. In fact, studies have found that melanoma cells can produce inflammatory cytokines that promote tumor growth and suppress the immune system.
The present article summarizes circadian gene alterations at the transcriptomic level and the impact of circadian rhythm on the immune infiltration level of cutaneous melanoma using publicly available transcriptomic databases.

Survival Analysis According to Gene Expression Level of the Core Circadian Genes in Skin Cutaneous Melanoma
The search for overall survival or disease-free survival biomarkers in melanoma remains desirable. Therefore, we performed an analysis of the core circadian genes in the GEPIA database. We found that higher expression levels of BMAL1, NR1D2 and NR1D2 were associated with disease free survival (DFS). Opposite results were found for CSNKE1, where low expression at marginal significance was associated with DFS ( Figure  3). Similar associations were found for the overall survival (OS), where significantly higher expression levels of NR1D2 and BMAL1 were associated with longer survival among skin cutaneous melanoma patients. Otherwise, lower expression levels of HLF and CSNK1E were significantly associated with OS among SKCM patients ( Figure 4). The further Multivariate Cox regression analysis performed by the TIMER database showed that the origin of SKCM played a role. High hazard ratio of overall survival was linked to Analysis performed by the GEPIA database using log 2 (TPM + 1) for log-scale.

Survival Analysis According to Gene Expression Level of the Core Circadian Genes in Skin Cutaneous Melanoma
The search for overall survival or disease-free survival biomarkers in melanoma remains desirable. Therefore, we performed an analysis of the core circadian genes in the GEPIA database. We found that higher expression levels of BMAL1, NR1D2 and NR1D2 were associated with disease free survival (DFS). Opposite results were found for CSNKE1, where low expression at marginal significance was associated with DFS ( Figure 3). Similar associations were found for the overall survival (OS), where significantly higher expression levels of NR1D2 and BMAL1 were associated with longer survival among skin cutaneous melanoma patients. Otherwise, lower expression levels of HLF and CSNK1E were significantly associated with OS among SKCM patients ( Figure 4). The further Multivariate Cox regression analysis performed by the TIMER database showed that the origin of SKCM played a role. High hazard ratio of overall survival was linked to TIMELESS and RORA expression in primary cutaneous melanoma patients, with values of 3.05 (1.36-6.846) and 2.29 (1.195-4.37), respectively (Table 1).       The groups were created based on the median values. Analysis was performed by the GEPIA database. Data for CSNK1A1L is not available.

Infiltration Level
Our study presents a correlation between circadian genes expression and immune infiltration in cutaneous melanoma. According to the in silico analysis performed using the TIMER database, there are many significant correlations between immune cell infiltration, including B cells, CD4+ T cells, CD8+ T cells, dendritic cells, macrophages and neutrophils, level and circadian gene expression levels in skin cutaneous melanoma patients. All the results are presented in Tables 2 and 3. The strongest correlation between circadian gene expression and the level of immune cell infiltration was observed among neutrophils and NR1D2, r = 0.52, p < 0.0001; BMAL1, r = 0.509, p < 0.0001; CLOCK, r = 0.45, p < 0.0001; CSNKA1A, r = 0.45, p < 0.0001; RORA, r = 0.44, p < 0.0001. Additionally, we found an association between mutation of selected circadian genes and expression of CTLA4, CD86, CD80 and CD28. We observed that patients with mutated PER2 presented higher expression levels of crucial CTLA4, CD86, CD80 and CD28 ( Figure 5). Moreover, PER2 expression level was correlated with CD86 and CD28 expression levels ( Figure 6B,C). Other interesting and significant correlations between PDCD1 expression and CRY2, BMAL1 and CSNK1E expression levels were r = 0.13, p = 0.04; r = 0.13, p = 0.04; r = −0.20, p = 0.0003, respectively ( Figure 6A).

Infiltration Level
Our study presents a correlation between circadian genes expression and immune infiltration in cutaneous melanoma. According to the in silico analysis performed using the TIMER database, there are many significant correlations between immune cell infiltration, including B cells, CD4+ T cells, CD8+ T cells, dendritic cells, macrophages and neutrophils, level and circadian gene expression levels in skin cutaneous melanoma patients. All the results are presented in Tables 2 and 3. The strongest correlation between circadian gene expression and the level of immune cell infiltration was observed among neutrophils and NR1D2, r = 0.52, p < 0.0001; BMAL1, r = 0.509, p < 0.0001; CLOCK, r = 0.45, p < 0.0001; CSNKA1A, r = 0.45, p < 0.0001; RORA, r = 0.44, p < 0.0001. Additionally, we found an association between mutation of selected circadian genes and expression of CTLA4, CD86, CD80 and CD28. We observed that patients with mutated PER2 presented higher expression levels of crucial CTLA4, CD86, CD80 and CD28 ( Figure 5). Moreover, PER2 expression level was correlated with CD86 and CD28 expression levels ( Figure 6B,C). Other interesting and significant correlations between PDCD1 expression and CRY2, BMAL1 and CSNK1E expression levels were r = 0.13, p = 0.04; r = 0.13, p = 0.04; r = −0.20, p = 0.0003, respectively ( Figure 6A).

Mutations of Circadian Genes in Cutaneous Melanoma
According to the cBioPortal and TIMER 2.0 database, the most mutated circadian genes among cutaneous melanoma patients were PER1, PER2, RORB and RORC (Table 4). Additionally, these mutations were most frequent among skin cutaneous melanoma compared to other cancer patients. Types of mutations that occurred among circadian genes included inframe mutations, missense mutations, splice mutations, truncating mutation, structural mutation and amplification, as well as deep deletion (Figure 7). Some of these mutations had reflections in mRNA expression of circadian genes in cutaneous melanoma patients. Mutations of BMAL1 and NR1D2 were associated with the increased expression of their mRNA levels ( Figure 8). The analysis of overall survival among cutaneous melanoma patients showed that the patients with at least one alteration within the selected circadian genes (24 genes) had a worse prognosis compared to the patients without any alterations in these genes (Figure 9). The unaltered group of patients lived 52 months longer than the altered group of patients (logrank test p-value = 0.03).

Methylation of Circadian Genes
Correlation between mRNA level of the selected circadian genes and methylation level was analyzed using thecBioPortal database. The strongest negative correlation between the transcription level and methylation was observed for PER3, RORC, CSNK1E and PER2. Otherwise, the strongest positive correlation was observed for RORA and HLF (Table 5). Moreover, we found that methylation of selected CpG sites of circadian genes was associated with overall survival of skin cutaneous melanoma patients. The in silico analysis showed that higher methylation levels of BHLHE41 (cg03046445), CLOCK (cg05960024), CSNK1E (cg01346718), CSNK1E (cg01441777), PER1 (cg16545079), PER2 (cg05664072), PER3 (cg25514503), RORC (cg18149207) and RORA (cg27167601) were correlated with longer overall survival ( Figure 10). Moreover, significant relationships between methylation and OS were reported for CpG islands located in the transcription start site (TSS) ( Table 5).

Discussion
The presented in silico research comprehensively analyzed the mRNA expression levels of 24 circadian genes and correlations with the immune infiltration level in skin cutaneous melanoma by the use of multipleomics databases.
Cutaneous melanoma is the leading cause of skin cancer deaths and is typically diagnosed at an advanced stage, which results in a poor prognosis. Therefore, new genetic targets for therapy and valuable prognostic biomarkers are desirable [10]. Our study showed disruptions of circadian genes in broad aspects. The presented in silico analysis showed that more than half of the investigated circadian genes represented significantly different expression patterns compared to normal skin. We found that BHLHE41 and TIMELESS were upregulated in skin cutaneous melanoma patients, while BMAL1, BHLHE40, CRY2, HLF, NFIL3, NR1D1, PER1, PER2, PER3, RORA, RORC and TEF were downregulated. Moreover, a similar tendency for these genes was associated with more advanced stages (stages III and VI) of SCKM. Lengyel et al. have shown that mRNA levels of core circadian genes, including PER1, PER2, CLOCK and CRY1, were reduced in melanoma skin biopsies compared to adjacent normal skin, but in a limited number of patients [11]. Moreover, a later study conducted by de Assis et al. showed that patients with more advanced tumor stages (III and IV) had a lower PER3 expression [12]. The presented in silico analysis of a larger group of patients partially proved that observation. We did not observe significant differences for mRNA levels of CLOCK and CRY1 genes. An animal study conducted by de Assis et al. demonstrated that circadian gene expression was altered in melanoma compared to normal skin, pointing to the reduction in Per1 and Bmal1 in tumor melanomabearing mice [13]. Additionally, the same authors have shown that circadian genes were disrupted in metastatic cutaneous melanoma based on an in silico analysis [12]. We found a negative correlation between the transcript level of circadian genes and their methylation status, mostly in the case of CpG islands located in Transcription Start Site. An altered epigenetic pattern of circadian genes has been observed in cancer studies [14]. Methylation status of cfDNA cutaneous melanoma is considered as an early cancer detection biomarker. Our analysis showed that OS of SKCM patients was associated with methylation status of CpG site of circadian genes.
Additionally, mutations of circadian genes are also present in SKCM; the most mutated genes are PER2 and PER3, which are classified as tumor suppressor genes [15][16][17]. Moreover, patients with at least one alteration among circadian genes have decreased overall survival. These observations may indicate that disruptions of circadian genes may play an essential role in the induction of cancer formation [18]. The link between circadian genes and cutaneous melanoma may be explained by the role of circadian genes in DNA damage repair, cellular proliferation and apoptosis [9]. In 2022, for example, a study by de Assis et al. found that increased Bmal1 gene expression was associated with reduced melanoma cell proliferation [19]. Therefore, in silico analysis may be a powerful tool for investigating circadian genes dependent strategy for anticancer therapy of skin.
An observation suggesting that circadian genes may be used as potential prognostic biomarkers in SCFCM and are associated with poor prognosis in cutaneous melanoma patients (including disease-free survival (DFS) and overall survival (OS) rates) is another interesting result of the present study. A similar association has been found by de Assis et al., where BMAL1 was positively correlated with overall survival in metastatic melanoma patients. High expression of BMAL1 was associated with longer survival [12]. This observation was proven for a larger group of patients in our study. Other studies have already identified circadian genes as prognostic markers for other tumors. Triple-negative breast cancer patients with a high expression level of NR1D1 treated with chemotherapy had longer OS [20]. Our analysis showed that an increased NR1D2 expression level was correlated with positive results, including OS and DFS, in SKCM patients. A study on stomach adenocarcinoma has revealed opposite results, where high NR1D1 and PER1 levels were correlated with poor OS [21]. Another study has also shown that decreased levels of PERs genes were associated with poor OS rates in non-small cell lung cancer patients [22]. Additionally, CLOCK has been linked with a worse prognosis for hepatocellular carcinoma [23].
Another interesting aspect of our study is related to the infiltration level of many immune cells, which is correlated with the majority of the investigated circadian genes. The tumor microenvironment (TME) plays a significant role in tumorigenesis and SKCM progression, but the dynamic regulation of the immune and stromal components is not yet fully understood [24]. The tumor microenvironment (TME) involves many types of cells, including tumor cells, fibroblasts, stromal cells, immune and inflammatory cells including natural killer cells (NK), tumor-associated macrophages (TAM), natural killer neutrophil T cells (NKT), dendritic cells (DC), T lymphocytes, B lymphocytes and the extracellular matrix (ECM). Immune cells in the TME are key factors affecting tumor growth and treatment response. The TME influences clinical outcomes and contains potential targets for therapeutic modulation [25]. In the present study, we quantified the correlation between circadian genes' expression and the proportion of tumor-infiltrating immune cells, based on the TIMER 2.00 database [26,27].
Recently, increasing evidence has indicated that the tumor microenvironment (TME) is involved in tumor development. Interactions between cancer cells, stromal cells and tumorinfiltrating immune cells (TICs) are critical for malignant cancer progression, including promoting replicative immortality, invasion, metastasis and evasion of immune surveillance [28].
Multiple studies have reported that TIC represents a promising TME marker for assessing treatment effects. TIC components and their activation statuses are important parameters affecting patient prognosis and tumor characteristics. Anti-cytotoxic T lymphocyte antigen 4 (CTLA-4) therapy activates T cells and induces programmed deathligand 1 (PD-L1) expression in tumor cells and TICs [29]. Activation of CD8+ T cells can prolong patients' survival in many cancers, including SKCM [30]. A study has shown that increased CD8+ T cell trafficking contributes to the efficacy of anti-programmed death 1 (PD-1)/CTLA-4 therapy against melanoma metastasis and may represent an effective immunotherapeutic strategy [31]. Neutrophils also play a background role in melanoma and can actively switch to an antitumor mode [32]. Our analysis showed that the strongest correlation between circadian genes' expression and immune infiltration was observed for neutrophils. The abovementioned studies demonstrate that crosstalk between cancer cells and the TME plays an integral role in the CM development, making it challenging to accurately characterize the dynamic regulation of the TME by immune and stromal components [29]. Furthermore, the circadian clock regulates differentiation of CD4(+) T helper [33]. We found a significant correlation between the CD4+ T cell infiltration level and mRNA expression of CRY2, DBP, PER2 and TEF ( Table 4). The presented results also revealed that crucial cytokines, chemokine and other crucial immune biomarkers are significantly associated with circadian gene expression in both primary and metastatic melanoma (Table 3).
Aberrant inflammation, immunosuppression and evasion of immune surveillance promote tumorigenesis by increasing cell proliferation, cell survival, angiogenesis and invasiveness. Both the immune system and inflammatory processes are regulated by the circadian clock and its components [34,35]. NF-κB is one of the key regulators of inflammation, regulated by circadian clock components BMAL1, CLOCK and CRY [36]. BMAL1 is physiologically downregulated during activation of NF-κB-related inflammatory responses through p65 subunit phosphorylation, but our analysis showed a positive correlation between the expression of NF-κB and BMAL1 in both primary and metastasis cutaneous melanoma. On the other hand, CLOCK binds to p65 and increases NF-κB-induced transcription of inflammatory genes. CRY regulates NF-κB activation by modulating adenylate cyclase, an enzyme that regulates cAMP levels and increases protein kinase A (PKA). This, in turn, promotes NF-κB through p65 phosphorylation -κB activation, leading to a lack of CRY and reduced expression of NF-κB targets such as TNF, IL6 and CXCL1 in fibroblasts and macrophages. A similar observation has been indicated in the study where CRY1 had a strong positive correlation with the NF-κB expression level [37].
Increasing evidence indicates that the assessment of immune infiltration may predict the outcome of cancer patients. Recent studies on metastatic melanoma have shown that a high level of BMAL1 is manifested by an increased T cell infiltration level [12]. This correlation is also shown in the current analysis. Interestingly, de Asiss et al. have shown that patients with a decreased BMAL1 mRNA level had a worse response to anti-PD-1 therapy compared to those with the elevated BMAL1 expression [12]. The present analysis showed a positive correlation between expression of BMAL1 and PDCD1 in SKCM patients. The immune system is controlled by circadian rhythm and plays a significant role in the immune cells development. REV-ERB receptors regulate adhesion and motility of macrophages [38]. Rodent studies have shown that lack of Clock gene is manifested by a reduced number of Th1 cells. BMAL1 and CRY are involved in B-cell development [39].
In silico genomic analysis, which involves the use of computational models and simulations to study genomics and gene expression, has some limitations. The existing data may be incomplete or biased, leading to potential inaccuracies in the modeling and analysis. Therefore, in silico models and simulations need to be validated and confirmed through experimental studies. However, conducting experiments to validate every aspect of the model can be time-consuming, expensive or even ethically challenging. Consequently, the validation of in silico models may be limited and the predictions or findings may need to be interpreted with caution. Biological systems are dynamic and gene expression can be influenced by various factors, including environmental stimuli and cellular context. In silico models may struggle to accurately capture the dynamic nature of biological processes, especially when the underlying regulatory mechanisms are not well understood or easily modeled.
It is important to note that the roles of circadian genes in SKCM are complex and can vary depending on various factors. The research in this area is ongoing and further studies are needed to fully elucidate their precise roles and potential therapeutic implications in melanoma.
In conclusion, the circadian rhythm study in cutaneous melanoma is an area of research that has the potential to result in new strategies for disease prevention and treatment. Disruptions in the circadian rhythm have been linked to an increased risk of cancer. It is also believed that understanding of the circadian rhythm in cutaneous melanoma could lead to new insights into the development and progression of the disease. By targeting the circadian rhythm and reducing inflammation, it might be possible to slow down the progression of cutaneous melanoma and improve patient outcomes. Moreover, the infiltration level of immune cells and their relationship with circadian rhythm in SKCM provides insights into the immune response dynamics, potential prognostic and predictive markers, optimization of treatment timing and the impact of circadian disruption on immune suppression. These findings have clinical implications for the development of targeted immunotherapies, chronoimmunotherapy strategies and the identification of biomarkers for treatment response.

GEPIA-Gene Expression Profiling Interactive Analysis Database
In our study, the in silico analysis was performed in relation to 24 core circadian genes (Table 1) in the samples of tissues derived from skin cancer cutaneous melanoma (SKCM) and normal skin. For the analysis, an interactive web server GEPIA was used. The GEPIA website was set up to analyze the RNA sequencing expression data of 9736 tumors and 8587 healthy tissue samples from the TCGA and the GTEx projects, using a standard processing pipeline. The GEPIA website provides customizable functions such as tumor/normal differential expression analysis, profiling according to cancer types or pathological stages, patient survival analysis, similar gene detection, correlation analysis and dimensionality reduction analysis [40,41].
In the present study, we used a tab multiple gene expression analysis. In this section, an analysis of gene expression core circadian genes SKCM using the box plot, stage and survival analysis was performed.

TIMER 2.0-Database
TIMER 2.0 (Tumor IMmune Estimation Resource 2.0) is a publicly available database that provides comprehensive information on the immune cell infiltration in various types of human cancers. The database integrates gene expression data from over 10.000 tumor samples from The Cancer Genome Atlas (TCGA) and the Genotype-Tissue Expression (GTEx) project, with established computational algorithms to estimate immune cell infiltration levels in the samples. This version of the webserver features immune infiltrates' abundances estimated by multiple deconvolution methods. The platform allows comprehensive exploration of immunological, clinical and genomic features in tumors. In our analysis, we used the modules that allow to check the correlation between circadian genes' expression and immune infiltration level of six major immune cell types: B cells, CD4+ T cells, CD8+ T cells, dendritic cells, macrophages and neutrophils. It also enables users to explore the relationship between immune cell infiltration and gene expression, as well as clinical outcomes such as survival rates and response to immunotherapy. The correlation value between gene expression and immune infiltration level was assessed by the Spearman correlation and adjusted by tumor purity. Distributions in gene expression levels are displayed using box plots, with statistical significance of differential expression evaluated using the Wilcoxon test. TIMER 2.0 is accessible free of charge at http://timer.cistrome.org/ (accessed on 21 March 2023) [26,27].

cBioPortal
cBioPortal (https://www.cbioportal.org/ accessed on 21 March 2023) is a web-based platform for exploring and analyzing large-scale cancer genomics data sets. It provides access to molecular data, such as DNA mutations, mRNA and microRNA expression levels, protein expression and DNA copy number alterations from over 200 cancer studies. The platform enables researchers to visualize and analyze data using a variety of tools, including interactive heatmaps, scatterplots and survival curves. The cBioPortal was developed by the Memorial Sloan Kettering Cancer Center and the Dana-Farber Cancer Institute and is freely available to researchers and clinicians around the world [42,43]. The data analysis was performed on a dataset obtained from the TCGA PanCancer Atlas, including 448 samples of Skin Cutaneous Melanoma. Figure 11 presents basic features of SCCM patients included in the analysis.