Role of Melatonin-Related Signatures in the Prognosis of Glioblastoma

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

Download PDF

Research Article

Role of Melatonin-Related Signatures in the Prognosis of Glioblastoma

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

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

Melatonin (MT) can have a direct or indirect impact on the expression of numerous crucial genes in cancer cells. Consequently, this paper strived to investigate MT-associated prognostic genes in Glioblastoma (GBM) and their mechanisms of operation.

Methods

Based on the survival of samples in training set and the expression of 34 MT-related genes (MT-RGs), univariate Cox regression analysis and consistent cluster analysis were used to identify differential prognostic MT-RGs for follow-up analysis. Then, prognostic genes were obtained by least absolute shrinkage and selection operator (LASSO) regression analysis to construct risk model. On the basis of risk model, GBM patients in training set were sorted to high- and low-risk cohorts. Further, to predict patient survival, we combined risk score with other clinical characteristics to access independent prognostic factors for the construction of prognostic model. Moreover, immune cell infiltration analysis and drug sensitization were performed to probe the potential mechanism of prognostic genes. Additionally, single-cell analysis was utilized to recognize cell types in GBM patients, and pseudotemporal trajectories of predictive genes were predicted in key cell types.

Results

The GBM patients in the training set were divided into two clusters (Cluster 1 and Cluster 2). Based on two clusters and the expression of MT-RGs, 8 differential prognostic MT-RGs were selected by wilcoxon test. Ultimately, ECE1, CYP1B1, IRAK1, SIRT1, and ACHE were defined as prognostic genes via LASSO. Importantly, risk score and age associated with GBM prognosis were identified as independent prognostic factor to establish the prognostic model, which might be an effective tool for predicting survival in GBM patients. For immune analysis, neutrophils were markedly more abundant in high-risk cohort. And 198 drugs (e.g. Temozolomide, Cispiatin, Etoposide) were identified between two risk cohorts. Moreover, 7 cell types (Myeloid, Neoplastic, Oligodendrocytes, Astrocytes, Vcascular, OPCs, and Neurons, and Myeloid) were annotated in single-cell dataset. Pseudotemporal analysis revealed that ECE1, CYP1B1 and ACHE expression steadily decreased on Myeloid, Vascular and Neoplastic.

Conclusion

A new prognostic model for GBM associated with MT was created and verified, which might furnish new opinions for the study of MT-RGs in GBM.

Glioblastoma

Melatonin

prognosis

immune

single-cell

Glioblastoma (GBM) is one of the most malignancies of the central nervous system (CNS) and is associated with high rates of recurrence, morbidity, and mortality¹. Despite comprehensive intervention involving surgical resection, radio, and chemotherapies, patients often have a very poor outcome². Recently, a myriad of biological indicators were identified to facilitate the accurate diagnosis and prognosis of patients with multiple cancers. Globally, in 2021, the World Health Organization (WHO) revised their stratification of CNS tumors based on morphology and molecular variables, which indicates that molecular testing is crucial in GBM diagnosis³. However, despite much improvement in the molecular profile-based diagnosis and prognosis, patient outcomes are still unsatisfactory⁴. Therefore, novel and curative treatment strategies are urgently needed.

Melatonin (MT), or N-acetyl-5-methoxytryptamine, is a major secretory product of the pineal gland, but is also produced in extrapineal organs including the retina, skin, and ovary.^5,6 It is not only a well-known modulator of reproductive physiology in mammals, but also has a variety of biological activities, such as anti-inflammation, anti-oxidation, anti-convulsion, anti-anxiety, and application of oncostatic and tumorinhibitory effects^7,8. Melatonin reduces glioma cell proliferation both in vitro and in vivo, according to previous studies.^9,10 Notably, melatonin decreases GBM stem-like cells proliferation, self-renewal and clonogenic capabilities through the decrease of stem cell markers and the increase of differentiation markers^10,11. There is evidence that melatonin plays a role in reprogramming events involving transcriptional competency acquisition and gene activation in these cellular changes.

Melatonergic system axis signaling has a strong correlation with various genes. Colombo et al. reported that melatonin functioned as a suppressor of hepatocarcinoma angiogenesis by downregulating the expression of HIF-1 and VEGF¹². A recent study suggested that melatonin could play a pivotal role in tumor suppression via miR-424-5p/VEGFA axis in osteosarcoma¹³. Melatonin exerts an anti-angiogenic activity in HepG2 cells by interfering with the transcriptional activation of VEGF, via Hif1α and STAT3¹⁴. Thus, targeting the melatonergic system axis may have a synergistic anti-tumor effect with therapy or immunotherapy.

In this study, MT-RGs were selected through univariate Cox regression analysis. Based on that, we stratified GBM into two subgroups based on the expression patterns of MT-RGs. In addition, the correlation of MT-RGs with immune infiltration and drug sensitivity was evaluated. This study highlights a functional role for the MT-RGs signature and uncovers a potential prognostic biomarker for GBM, providing novel insights into potential therapeutic targets and strategies for the treatment of GBM.

2.1 Data source

The Cancer Genome Atlas (TCGA)-GBM cohort, consisting of 168 patient samples and 5 healthy samples, was retrieved from TCGA database (https://portal.gdc.cancer.gov/). The external validation set was downloaded from the Chinese Glioma Genome Atlas (CGGA) database (http://www.cgga.org.cn/index.jsp), which contained 85 GBM samples. The relevant gene expression matrices and the corresponding clinical data in the validation set were utilized to validate the risk model. Single-cell dataset GSE84465 (GPL18573) with GBM tumor samples was accessed from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). A total of 34 melatonin-related genes (MT-RGs) were mined from previously published literature ¹⁵.

2.2 Acquisition of differential prognostic MT-RGs

According to samples survival in training set and expression of 34 MT-RGs, univariate Cox regression analysis was first performed to identify MT-RGs associated with GBM prognosis (prognosis MT-RGs) using survival package (v 3.3-1)¹⁶, with HR ≠ 1 and P < 0.2 as screening conditions. Immediately after, correlations between prognosis MT-RGs were calculated by Spearman using corrplot package, and a correlation network map was drawn to demonstrate the complex relationships between prognosis MT-RGs. On the basis of the expression of prognosis MT-RGs, the GBM samples in the training set were consistently clustered via ConsensusClusterPlus package (v 1.62.0)¹⁷, and the optimal number of clusters was determined by the Cumulative Distribution Function (CDF). The survminer package (v 0.4.9) was then utilized to compare the differences in survivability between clusters. Then, prognosis MT-RGs were tested for differences between clusters using Wilcoxon method to get differential prognostic MT-RGs. Meanwhile, the association of different clusters and clinical characteristics of differential prognostic MT-RGs was presented in a heat map by pheatmap package (v 1.0.12)¹⁸. For further investigation of the functions of differential prognostic MT-RGs, enrichment analysis depending on Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) databases was performed through clusterProfiler package (v 4.2.2)¹⁹ and human gene annotation package (org.Hs.eg.db, v 3.14.0) with P.adj < 0.05. In addition, we used the STRING website (https://string-db.org) to construct a protein-protein interaction (PPI) network for differential prognostic MT-RGs to investigate whether there were interactions between these genes. The PPI network was then beautified using Cytoscape software. On the other hand, differentially expressed genes (DEGs) inter-clusters were selected by DESeq2 package (v 1.31.0)²⁰ with P.adj < 0.05 and |log₂FC|>1.

2.3 Construction and validation of risk model

According to differential prognostic MT-RGs, least absolute shrinkage and selection operator (LASSO) regression analysis via glmnet package (v 4.1-3)²¹ was utilized to select out prognostic genes. And risk score was computed utilizing a formula below.

$$risk score=\sum _{i=1}^{n}\left(\beta i*xi\right)$$

β and x indicated coefficients and gene expression, respectively. Meanwhile, a risk model was built and evaluated in training set. First, the GBM samples of training set were sorted into high- and low-risk cohorts by median of risk score, and risk curves for two cohorts as well as distribution of patients survival were plotted with the help of survival package (v 3.3-1)¹⁶. Furthermore, Kaplan-Meier (K-M) curves were plotted using survminer package (v 0.4.9) to analyze differences in survival between two risk cohorts. Meanwhile, principal component analysis (PCA) was executed for high- and low-risk cohorts. Based on the risk scores previously computed for samples in training set, the true positives, false positives, and corresponding area under the curve (AUC) values for each sample were calculated using survivalROC package (v 1.0.3)²², and receiver operating characteristic (ROC) curves were plotted at 3-, 5-, and 7-year to test the predictive ability of risk model for a given survival time. Similarly, risk model was verified by GSE84465 with risk, K-M survival and ROC curves.

2.4 Independent prognostic analysis

Five indicators: risk score and 4 clinical characteristics (age, gender, MGMT status, and subtype) were sequentially analyzed by univariate Cox regression, PH hypothesis test, and multivariate Cox regression analysis. Then, the independent prognostic factors derived from the above analyses were applied to plot in a nomogram as a prognostic model for predicting patient survival in clinical. Additionally, the predictive ability of nomogram was established by calibration curves. Further, we also analyzed the relationship between risk score and 4 clinical characteristics.

2.5 Enrichment analysis of risk cohorts

Based on KEGG gene set, we performed gene set enrichment analysis (GSEA) depending on log₂FC with threshold settings |NES| ≥ 1 and P.adj < 0.05 for all genes in high- and low-risk cohorts.

2.6 Analysis of immune cell infiltrates and drug sensitivity

Immune cell infiltration in GBM samples from training set was evaluated with CIBERSORT algorithm. Then, immune cells with an abundance of 0 in more than 75% samples were filtered out, and remained immune cells were utilized for further variance analysis in high- and low-risk cohorts by Wilcoxon test. And correlations between immune cells were analyzed through Spearman. Additionally, we analyzed the correlation between prognostic genes and immune infiltrating cells. To investigate potential drugs for GBM, the 50% inhibitory concentration (IC₅₀) of 198 drugs was extrapolated from Cancer Genome Project (CGP) via pRRophetic package²³ and then compared between high- and low-risk cohorts of training set using Wilcoxon test (P < 0.05).

2.7 Single-cell analysis

In single-cell data, Seurat package (v 4.3.0)²⁴ was applied for single-cell data and utilized to filter cells with less than 200 genes and genes covered by less than 3 cells firstly. The data were then quality controlled (QC) for nFeature_RNA number (greater than 300 and less than 6,000) and nCount_RNA number (less than 1,500,000). Next, the data was normalised using NormalizeData function. Then, highly variable genes were extracted by FindVariableFeatures function. After data normalization by means of the ScaleData function, Principal Component Analysis (PCA) downscaling was performed to select Principal Components (PCs). Unsupervised cluster analysis was performed to gain clusters according to Seurat standard procedure and the outcome was demonstrated with UMAP (resolution = 0.5). With the Human Primary Cell Atlas Data employed as the reference dataset, clusters were annotated for cell subtypes using singleR package (v 1.0.6)²⁵, and marker genes and cell types were identified further through literature²⁶ in conjunction with the singleR results. And histograms were used to show the proportions of different cell types in GBM samples. In addition, the expression of prognostic genes on different cell types was also shown using DotPlot function. Pseudotemporal on key cell types were simulated using monocle3 package (v 2.24.1) ²⁷. Furthermore, gene set variation analysis (GSVA) analysis of annotated cell types was utilized to obtain score and pathways based on both hallmark gene set and KEGG gene set using GSVA package (v 1.44.5)²⁸.

2.8 Statistical analysis

R language was applied to execute all analyses. And P < 0.05 was considered as statistically significant.

3.1 Dual analysis of melatonin-related genes expression identifies two distinct subgroups of GBM

To find out whether genetic variants in melatonin related genes(MT-RGs)are associated with GBM, 34 MT-RGs were included in this study according to the previous study¹⁵, which are summarized in Table S1., To explore the survival significance of MT-RGs, 12 genes were selected as key factors from the above 34 MT-RGs by univariate Cox regression analysis (Fig. 1A) Subsequently, Spearman correlation analysis revealed complex relationships among the 12 prognosis MT-RGs, for example, ECE1 had a strongly negative correlation with MAP2 (Fig. 1B). To stratify GBM based on their relative expression levels of MT-RGs, we adopted the classic unsupervised clustering algorithm to perform remodeling analysis of TCGA-GBM. We chose k = 2 as the optimal cluster number through the R package and confirmed it using the PAC algorithm. GBM samples were specifically finalized into 2 subtypes, Cluster 1 and Cluster 2 (Fig. 1C-E), the GBM samples K-M survival curves revealed that the survival rate of Cluster 1 was lower than that of cluster 2 in the training set (Fig. 1F). Besides, the differential analysis showed that there were 1,373 DEGs inter-Clusters, among which 876 were up-regulated and 497 were down-regulated (Fig. S1).

3.2 Difference of biological features between two MT-RTs clusters

We compared the gene expression and clinical information of the two clusters, found that the expression of MT-RGs differed significantly, with CYP1B1, ECE1, IRAK1 and NFKB1 were up-regulated and notably expressed, as well as ACHE, CSNK1E, MAP2 and SIRT1 were down in Cluster 1(Fig. 2A). Meanwhile, the PPI network demonstrated 16 interactions with 8 proteins corresponding to differential prognostic MT-RGs, such as MAP2-ACHE, ECE1-NFKB1, and IRAK1-SIRT1 (Fig. 2B). Furthermore, the heat map showed the expression profiles of 8 differential prognostic MT-RGs in two clusters and various clinical characteristics (Fig. 2C). Besides, differential prognostic MT-RGs were enriched in 391 items in Gene ontology (GO) analysis (‘negative regulation of NF-kappa B transcription factor activity’, ‘cellular response to tumor necrosis factor’, ‘intrinsic apoptotic signaling pathway in response to oxidative stress’, ‘rhythmic process’, etc.) and 92 pathways in KEGG (‘NF-kappa B signaling pathway’, ‘Neurotrophin signaling pathway’, ‘microRNAs in cancer’ and so on) (Fig. 2D-E).

3.3 Construction and validation of the prognostic MT-RGs score model

To predict the outcome of each patient with GBM, we created a scoring model based on prognosis-related DEGs between two MT-RGs clusters, called MT-RGs score (risk score). Totally 5 prognostic genes (ECE1, CYP1B1, IRAK1, SIRT1 and ACHE) were selected out by LASSO (λ = 0.053) according to 8 differential prognostic MT-RGs (Fig. 3A-B). Meanwhile, risk score was counted for constructing the risk model (Table S2). Subsequently, risk curve and survival distribution in high-risk cohort (n = 76) and low-risk cohort (n = 91) were charted in Fig. 3C. K-M curves demonstrated that GBM patients in high-risk cohort had lower survival probability (P < 0.0001) (Fig. 3D). At the same time, the PCA findings indicated a separation between two risk cohorts (Fig. 3F). By ROC curves, AUC values of 3-, 5- and 7-year were separately 0.689, 0.721 and 0.900, revealing a good predictive ability of the risk model for specific survival times (Fig. 3G). To validate the risk model, patients in validation set were also divided into a high-risk cohort (n = 62) and a low-risk cohort (n = 23) (Fig. S2A). As with the training set, patients in high-risk cohort had lower survival probability (P = 0.013) (Fig. S2B), and there was a clear distinction between two cohorts by PCA (Fig. S2C). In addition, AUC values of 3-, 5- and 7-year were all exceeded 0.600, indicating the risk model could predict survival favorably (Fig. S2D).

3.4 Association Between the Risk Score and Clinical Features

Next, we explored whether the MT-RGs score were associated with GBM patient prognosis. Univariate Cox regression analysis demonstrated age, MGMT Promoter Methylation, and risk score were associated with patient survival and prognosis (Fig. 4A). After PH test, age, MGMT Promoter Methylation, and risk score were used to construct a multivariate Cox model (Fig. 4B). Ultimately, risk score (P = 0.004) and age (P = 0.003) were considered to be independent prognostic factors for the plotting of nomogram (Fig. 5A). Calibration curves demonstrated the ability of the nomogram to predict survival at 3-, 5-, and 7-year, suggesting that the nomogram might be a useful tool for predicting survival in GBM patients (Fig. 5B-D). Further, analysis of the relationship between age, gender, subtype, MGMT Promoter Methylation and risk score showed that there were notable differences in risk score among different subgroups in subtype (Fig. S3A-D). In addition, GSEA indicated that risk cohorts enriched for ‘cell cycle’, ‘chemokine signaling pathway’, ‘cytokine-cytokine receptor interaction’ and so on in KEGG (Fig. S4).

3.5 Impact of ACHE on GBM Immunity.

The abundance of the 22 immune cells was proportionally distributed across high- and low-risk cohorts (Fig. 6A). After filtering out, variance analysis of the remaining 17 immune cells revealed that neutrophils were markedly more abundant in high-risk cohort than in low-risk cohort, while abundance of activated NK cells was significantly higher in low-risk cohort (Fig. 6B). In addition, Spearman correlation analysis was performed to calculate the correlation between 17 immune cells (Fig. 6C). The correlation between resting NK cells and regulatory T cells (Tregs) was positive (cor = 0.46) and notable (P < 0.05), which represented a certain synergy between these two types of immune cells in their functioning. Resting mast cells were negative (cor=-0.79) and remarkable (P < 0.05), representing some antagonism when functioning. The results of correlation between 5 prognostic genes and 17 immune cells showed that ACHE was the strongest positively and significantly correlated with M0 macrophages (cor = 0.278) and the strongest negatively and significantly correlated with M2 macrophages(cor=-0.293) (Fig. 6D).

GBM is highly resistant to chemotherapeutic or immunotherapy drugs. To further explore the drug sensitivity between two risk cohorts. By drug sensitivity analysis, a total of 198 drugs were identified between two risk cohorts. Temozolomide (Fig. 6E) along with 173 drugs (e.g. Bleomycin, Cispiatin, Etoposide) exhibited higher sensitivity in high-risk cohort compared to low-risk cohort (Fig. S5A), whereas only 13 drugs were sensitive lowly in high-risk cohort (Fig. S5B), which the effect of chemotherapeutic drugs was better in the low-risk groups.

3.6 Seven cell types were annotated in single-cell dataset

In GSE84465, 3,571 cells and 20,809 genes satisfied the criteria (Fig. 7A), resulting in 3,409 cells and 20,809 genes after QC (Fig. 7B). Then 2000 highly variable genes and top 20 PCs were selected for subsequent analysis (Fig. 7C-E). Cells were clustered into 15 clusters using UMAP (Fig. 7F) Combined with singleR and literature [PMID: 29091775], 7 cell types and corresponding marker genes were further identified (Fig. 7F-H), which were Myeloid (PTPRC), Neoplastic (EGFR), Oligodendrocytes (MOG), Astrocytes (AGXT2L1), Vcascular (DCN), OPCs (GPR17), and Neurons (STMN2), and Myeloid was the most abundant population (Fig. 7I).

Expression of 5 prognostic genes in the cell types showed that ECE1 was higher expressed in Neoplastic, IRAK1 in Myeloid, CYP1B1 and SIRT1 in Vcascular, as well as ACHE in Neurons (Fig. 8A-B). Pseudotemporal analysis was conducted on key cells (Myeloid, Vascular and Neoplastic) revealing that ECE1, CYP1B1 and ACHE expression steadily decreased, whereas IRAK1 expression exhibited a gradual increase followed (Fig. 8C-D). Additionally, SIRT1 expression displayed a gradual increase in pseudotemporal. Moreover, GSVA revealed that the score of Myeloid was highly enriched in ‘allograft rejection’, ‘IL6 JAK STAT3signaling’, ‘inflammatory response’ and so forth in Hallmark (Fig. S6) and ‘ribosome’, ‘cytokine-cytokine receptors interactions’, ‘glycosphingolipid biosynthesis globo series’ etc. pathways in KEGG (Fig. S7).

GBM is the most aggressive type of brain tumour with a poor prognosis²⁹. Hence, the development of molecular mechanisms of GBM is essential. Melatonin axis signaling plays a wide role in the development and progression of cancer³⁰. As reported previously, melatonin reduces the morbidity and associated mortality of fatal viral infections, including COVID-19, due to its antioxidant, anti-inflammatory, and immunomodulatory activities³¹. The relationship between melatonin and GBM has rarely been explored, moreover, melatonin-based molecular subtypes in tumor repertoires was less frequent. In the present study, we made a comprehensive summarization of melatonin regulator genes via in silico analysis. We found that these signatures were significantly dysregulated in GBM and were associated with genomic and epigenetic modification. According to the expression pattern of these signatures, GBM patients could be subgrouped into two significantly different melatonin-associated subtypes (Cluster 1 and Cluster 2). Compared with Cluster 1, the Cluster 2 exhibited higher tumor mutation burden, and poorer survival. Furthermore, we established a reliable risk model based on subtype biomarkers and validated the model using external datasets.

Melatonergic system axis signaling have a synergistic anti-tumor effect with targeted therapy or immunotherapy. We compared the gene expression of the two clusters, we created a scoring model, called MT-RGs score (risk score), indicating that the higher the MT-RGs score, the worse the prognosis. Our results show remarkable differences in genomic alterations between the low-risk and high-risk groups, with the expression of ECE1, CYP1B1, IRAK1, being higher in the high-risk group, while the expression of ACHE and SIRT1 was higher in the low-risk group. As, previous studies, CircNDC80 promotes GBM multiforme tumorigenesis via the miR-139-5p/ECE1 pathway³², and ECE1 along with PODXL, ICAM1, ALCAM1, CD97, and CD44 contributing to breast cancer chemotherapy evasion and metastatic seeding³³. Qimei Lin et al. reported that increased CYP1B1 expression in prostate cancer cells could promote tumor progression³⁴. A study by Chao Xi et al. showed that miR-146a-5p enhances tumor epithelial-mesenchymal transition through disinhibition of the TRAF6-IRAK1 complex and IKK-dependent NF-κB signaling pathway³⁵. Besides its role in synaptic transmission, ACHE also regulates multiple oncogenic signaling pathways involved in the classic function of tumors, as Kai Wang et al. proposed expression of the ACHE gene can significantly inhibit the proliferation of liver cancer cells in xenograft model³⁶. A study by Tan et al. showed that SRT1720 inhibits bladder cancer growth by inhibiting the SIRT1-HIF pathway³⁷. These results suggest that these MT-RGs have certain carcinogenic or anti-tumor functions to some extent. The effectiveness of MT-RGs score was also confirmed using CGGA datasets. Nonetheless, more datasets are needed in the future to verify the effectiveness of the model.

Both univariate and multivariate Cox regression analysis showed that MT-RGs score was an independent predictor of survival outcome in GBM patients. The ROC validated its predictive robustness for 3-, 5- and 7-year OS. Thus, MT-RGs score may have reliable predictive power for GBM patient prognosis. We also analyzed the correlation between MT-RGs score and immune cell infiltration. The results suggested that MT-RGs score was strongly correlated with immune cells. Neutrophils were markedly more abundant in high-risk cohort, while abundance of activated NK cells was significantly higher in low-risk cohort. Those cells played complex roles in tumor immunity. For example, a study suggested that enhanced neutrophils activity was associated with increased levels of IL-12p70 and correlated with worse patient outcome³⁸. Ji Liang et al. also found that neutrophils promote glioma growth via induction of S100A4 expression in glioma cells³⁹. NK cells have crucial roles in the innate immunosurveillance of cancer. The ex vivo activation, expansion and genetic modification of NK cells can greatly increase their anti-tumor activity and equip them to overcome resistance⁴⁰. Laughney AM et al. found metastatasizing lung adenocarcinoma cells that express high levels of the transcription factor SOX9 can better evade NK cell-mediated killing, due to increased MHC-I expression, whereas SOX2high cells were more sensitive to NK cell attack, through downregulation of MHC-I and upregulation of activating receptor ligands⁴¹. Poznanski SM et al. found expansion of NK cells with IL-21 or treatment with an activator of the Nrf2 antioxidant pathway promote Warburg-like metabolism with increased glycolysis and reduced oxidative phosphorylation, resulting in substantially improved resistance to the oxidatively stressed and nutrient-deprived TME and enhanced antitumour responses⁴². The above results revealed that immune monitoring function of patients in high-risk group was weakened, which was conducive to immune escape.

The resistance to chemotherapy in GBM is the most important cause of recurrence^43,44. Melatonin possesses the capability to sensitize various cancers, such as breast, lung, colon, hepatic, and hematologic cancers to the effects of antineoplastic agents like anthracyclines, alkylating agents, tyrosine kinase inhibitors, endogenous factors such as TNF-α, TRAIL, and FAS, as well as ionizing radiation⁴⁵. The development of resistance to chemotherapy in GBM is often a problem for physicians and patients. GBM patients with lower MT-RGs score were more sensitive to Temozolomide, Cisplatin, Cytarabine, Doxorubicin, Eesclomol, Eoposide, Gemcitabine, Methotrexate, Mitomycin, Nilotinib, Temsirolimus, Viblastine, Vorinostat and multiple chemotherapeutic drugs. As in previous studies, co-treatment of melatonin with vorinostat promotes the therapeutic sensitivity of GBM by inhibiting the expression of transcription factor EB⁴⁶. Our findings may provide more evidence for the follow-up study of MT-RGs and tumor resistance, which may help to reduce drug resistance and improve clinical outcomes.

In this work, we systematically explored the function of melatonin regulator genes and identified two molecular subtypes of GBM. There are also some limitations of our study. First, the main findings were based on bioinformatics analyses, which need to be further verified in subsequent experiments. Second, the risk model may be disturbed by some confounding factors. Thus, further independent and prospective datasets are warranted to confirm this risk model validity.

In summary, activating melatonergic system axis signaling may be an appropriate approach for the treatment of GBM. The melatonergic system axis may alter the immunosuppressive microenvironment and promote anti-tumor immune responses. Hopefully, the results of this study will facilitate individualized assessment of patient prognosis and provide guidance for drug therapy in clinical.

Author Contribution

MZ and ZJH conceived the study. MZ and SC performed data analyses. CHX and WXD prepared figures S5.MZ wrote the first draft of this manuscript. PYW and YGQ critically revised the manuscript. All the authors have approved the final version of the manuscript. All authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Berger, T. R., Wen, P. Y., Lang-Orsini, M. & Chukwueke, U. N. World Health Organization 2021 Classification of Central Nervous System Tumors and Implications for Therapy for Adult-Type Gliomas: A Review. JAMA Oncol 8, 1493-1501 (2022). https://doi.org:10.1001/jamaoncol.2022.2844
Schaff, L. R. & Mellinghoff, I. K. Glioblastoma and Other Primary Brain Malignancies in Adults: A Review. Jama 329, 574-587 (2023). https://doi.org:10.1001/jama.2023.0023
Pfister, S. M. et al. A Summary of the Inaugural WHO Classification of Pediatric Tumors: Transitioning from the Optical into the Molecular Era. Cancer Discov 12, 331-355 (2022). https://doi.org:10.1158/2159-8290.Cd-21-1094
Shergalis, A., Bankhead, A., 3rd, Luesakul, U., Muangsin, N. & Neamati, N. Current Challenges and Opportunities in Treating Glioblastoma. Pharmacol Rev 70, 412-445 (2018). https://doi.org:10.1124/pr.117.014944
Acuña-Castroviejo, D. et al. Extrapineal melatonin: sources, regulation, and potential functions. Cell Mol Life Sci 71, 2997-3025 (2014). https://doi.org:10.1007/s00018-014-1579-2
Venegas, C. et al. Extrapineal melatonin: analysis of its subcellular distribution and daily fluctuations. J Pineal Res 52, 217-227 (2012). https://doi.org:10.1111/j.1600-079X.2011.00931.x
Kim, K. J. et al. Melatonin suppresses tumor progression by reducing angiogenesis stimulated by HIF-1 in a mouse tumor model. J Pineal Res 54, 264-270 (2013). https://doi.org:10.1111/j.1600-079X.2012.01030.x
Vriend, J. & Reiter, R. J. Melatonin as a proteasome inhibitor. Is there any clinical evidence? Life Sci 115, 8-14 (2014). https://doi.org:10.1016/j.lfs.2014.08.024
Papagiannakopoulos, T. et al. Circadian Rhythm Disruption Promotes Lung Tumorigenesis. Cell Metab 24, 324-331 (2016). https://doi.org:10.1016/j.cmet.2016.07.001
Papantoniou, K. et al. Rotating night shift work and colorectal cancer risk in the nurses' health studies. Int J Cancer 143, 2709-2717 (2018). https://doi.org:10.1002/ijc.31655
Martín, V. et al. Involvement of autophagy in melatonin-induced cytotoxicity in glioma-initiating cells. J Pineal Res 57, 308-316 (2014). https://doi.org:10.1111/jpi.12170
Colombo, J., Maciel, J. M., Ferreira, L. C., RF, D. A. S. & Zuccari, D. A. Effects of melatonin on HIF-1α and VEGF expression and on the invasive properties of hepatocarcinoma cells. Oncol Lett 12, 231-237 (2016). https://doi.org:10.3892/ol.2016.4605
Vimalraj, S., Saravanan, S., Raghunandhakumar, S. & Anuradha, D. Melatonin regulates tumor angiogenesis via miR-424-5p/VEGFA signaling pathway in osteosarcoma. Life Sci 256, 118011 (2020). https://doi.org:10.1016/j.lfs.2020.118011
Carbajo-Pescador, S. et al. Inhibition of VEGF expression through blockade of Hif1α and STAT3 signalling mediates the anti-angiogenic effect of melatonin in HepG2 liver cancer cells. Br J Cancer 109, 83-91 (2013). https://doi.org:10.1038/bjc.2013.285
Jiang, A. et al. Deciphering the Role of Melatonin-Related Signatures in Tumor Immunity and the Prognosis of Clear Cell Renal Cell Carcinoma. Oxid Med Cell Longev 2023, 3077091 (2023). https://doi.org:10.1155/2023/3077091
Ramsay, I. S. et al. Model selection and prediction of outcomes in recent onset schizophrenia patients who undergo cognitive training. Schizophr Res Cogn 11, 1-5 (2018). https://doi.org:10.1016/j.scog.2017.10.001
Wilkerson, M. D. & Hayes, D. N. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 26, 1572-1573 (2010). https://doi.org:10.1093/bioinformatics/btq170
Gu, Z. & Hübschmann, D. Make Interactive Complex Heatmaps in R. Bioinformatics 38, 1460-1462 (2022). https://doi.org:10.1093/bioinformatics/btab806
Yu, G., Wang, L. G., Han, Y. & He, Q. Y. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics 16, 284-287 (2012). https://doi.org:10.1089/omi.2011.0118
Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014). https://doi.org:10.1186/s13059-014-0550-8
Li, Y., Lu, F. & Yin, Y. Applying logistic LASSO regression for the diagnosis of atypical Crohn's disease. Sci Rep 12, 11340 (2022). https://doi.org:10.1038/s41598-022-15609-5
Heagerty, P. J., Lumley, T. & Pepe, M. S. Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics 56, 337-344 (2000). https://doi.org:10.1111/j.0006-341x.2000.00337.x
Mihajlović, M. et al. Utilization of phosphate rock from Lisina for direct application: release of plant nutrients in the exchange-fertilizer mixtures. J Agric Food Chem 62, 9965-9973 (2014). https://doi.org:10.1021/jf502704j
Ritchie, M. E. et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43, e47 (2015). https://doi.org:10.1093/nar/gkv007
Aran, D. et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol 20, 163-172 (2019). https://doi.org:10.1038/s41590-018-0276-y
Darmanis, S. et al. Single-Cell RNA-Seq Analysis of Infiltrating Neoplastic Cells at the Migrating Front of Human Glioblastoma. Cell Rep 21, 1399-1410 (2017). https://doi.org:10.1016/j.celrep.2017.10.030
Trapnell, C. et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol 32, 381-386 (2014). https://doi.org:10.1038/nbt.2859
Hänzelmann, S., Castelo, R. & Guinney, J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 14, 7 (2013). https://doi.org:10.1186/1471-2105-14-7
Miller, K. D. et al. Brain and other central nervous system tumor statistics, 2021. CA Cancer J Clin 71, 381-406 (2021). https://doi.org:10.3322/caac.21693
Zhang, J. et al. Pan-cancer analyses reveal genomics and clinical characteristics of the melatonergic regulators in cancer. J Pineal Res 71, e12758 (2021). https://doi.org:10.1111/jpi.12758
Tan, D. X. & Hardeland, R. Targeting Host Defense System and Rescuing Compromised Mitochondria to Increase Tolerance against Pathogens by Melatonin May Impact Outcome of Deadly Virus Infection Pertinent to COVID-19. Molecules 25 (2020). https://doi.org:10.3390/molecules25194410
Wang, Y. et al. CircNDC80 promotes glioblastoma multiforme tumorigenesis via the miR-139-5p/ECE1 pathway. J Transl Med 21, 22 (2023). https://doi.org:10.1186/s12967-022-03852-3
Dashzeveg, N. K. et al. Dynamic Glycoprotein Hyposialylation Promotes Chemotherapy Evasion and Metastatic Seeding of Quiescent Circulating Tumor Cell Clusters in Breast Cancer. Cancer Discov 13, 2050-2071 (2023). https://doi.org:10.1158/2159-8290.Cd-22-0644
Lin, Q. et al. CYP1B1-catalyzed 4-OHE2 promotes the castration resistance of prostate cancer stem cells by estrogen receptor α-mediated IL6 activation. Cell Commun Signal 20, 31 (2022). https://doi.org:10.1186/s12964-021-00807-x
Xu, C. et al. MiR-146a-5p deficiency in extracellular vesicles of glioma-associated macrophages promotes epithelial-mesenchymal transition through the NF-κB signaling pathway. Cell Death Discov 9, 206 (2023). https://doi.org:10.1038/s41420-023-01492-0
Wang, K. et al. Protein liposomes-mediated targeted acetylcholinesterase gene delivery for effective liver cancer therapy. J Nanobiotechnology 19, 31 (2021). https://doi.org:10.1186/s12951-021-00777-9
Tan, P. et al. SRT1720 inhibits the growth of bladder cancer in organoids and murine models through the SIRT1-HIF axis. Oncogene 40, 6081-6092 (2021). https://doi.org:10.1038/s41388-021-01999-9
Rahbar, A. et al. Enhanced neutrophil activity is associated with shorter time to tumor progression in glioblastoma patients. Oncoimmunology 5, e1075693 (2016). https://doi.org:10.1080/2162402x.2015.1075693
Liang, J. et al. Neutrophils promote the malignant glioma phenotype through S100A4. Clin Cancer Res 20, 187-198 (2014). https://doi.org:10.1158/1078-0432.Ccr-13-1279
McArthur, D. R. Canines as removable partial denture abutments. Part II: Rest and undercut location for retainers. J Prosthet Dent 56, 445-450 (1986). https://doi.org:10.1016/0022-3913(86)90386-0
Laughney, A. M. et al. Regenerative lineages and immune-mediated pruning in lung cancer metastasis. Nat Med 26, 259-269 (2020). https://doi.org:10.1038/s41591-019-0750-6
Poznanski, S. M. et al. Metabolic flexibility determines human NK cell functional fate in the tumor microenvironment. Cell Metab 33, 1205-1220.e1205 (2021). https://doi.org:10.1016/j.cmet.2021.03.023
Lathia, J. D., Mack, S. C., Mulkearns-Hubert, E. E., Valentim, C. L. & Rich, J. N. Cancer stem cells in glioblastoma. Genes Dev 29, 1203-1217 (2015). https://doi.org:10.1101/gad.261982.115
Bao, S. et al. Glioma stem cells promote radioresistance by preferential activation of the DNA damage response. Nature 444, 756-760 (2006). https://doi.org:10.1038/nature05236
Asghari, M. H., Ghobadi, E., Moloudizargari, M., Fallah, M. & Abdollahi, M. Does the use of melatonin overcome drug resistance in cancer chemotherapy? Life Sci 196, 143-155 (2018). https://doi.org:10.1016/j.lfs.2018.01.024
Sung, G. J. et al. Inhibition of TFEB oligomerization by co-treatment of melatonin with vorinostat promotes the therapeutic sensitivity in glioblastoma and glioma stem cells. J Pineal Res 66, e12556 (2019). https://doi.org:10.1111/jpi.12556

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Role of Melatonin-Related Signatures in the Prognosis of Glioblastoma

Status:

Version 1

Abstract

Background

Methods

Results

Conclusion

Figures

1 Introduction

2 Materials and methods

2.1 Data source

2.2 Acquisition of differential prognostic MT-RGs

2.3 Construction and validation of risk model

2.4 Independent prognostic analysis

2.5 Enrichment analysis of risk cohorts

2.6 Analysis of immune cell infiltrates and drug sensitivity

2.7 Single-cell analysis

2.8 Statistical analysis

3 Results

3.1 Dual analysis of melatonin-related genes expression identifies two distinct subgroups of GBM

3.2 Difference of biological features between two MT-RTs clusters

3.3 Construction and validation of the prognostic MT-RGs score model

3.4 Association Between the Risk Score and Clinical Features

3.5 Impact of ACHE on GBM Immunity.

3.6 Seven cell types were annotated in single-cell dataset

4 Discussion

5 Conclusion

Declarations

Author Contribution

References

Additional Declarations

Supplementary Files

Status:

Version 1