Introduction
E2F transcription factors are referring
E2F1 to
E2F8 eight genes [
1]. The E2F family genes are critical to the development of cancer by regulations of DNA replication and cell cycle progression [
2‐
4]. Expression levels of E2F family genes were associated with the clinical outcomes of bladder cancer, prostate cancer, lung cancer, colon cancer or breast cancer [
5‐
7]. Previously, using The Cancer Genome Atlas datasets, we had analyzed the expressions and prognosis of E2F transcription factors across different tumor types [
8]. Our results suggested the unfavorable prognostic effects of E2F transcription factors in liver hepatocellular carcinoma (LIHC) and lung adenocarcinoma (LUAD). However, the analysis was focusing on adult tumor types. The expressions and prognosis of E2F transcription factors in pediatric cancer types, particular in pediatric neuroblastoma are unclear.
Neuroblastoma is a malignant pediatric disease with poor prognosis [
9,
10].
MYCN amplification is detected in approximate 25% neuroblastoma patients [
11].
MYCN amplification [
12] and
MYCN high expression [
13] represent adverse prognostic factors in neuroblastoma. E2F transcription factors are involved in the development of neuroblastoma by regulating
MYCN expression [
14]. Our previous results showed that
E2F1 was a target of
MYCN amplification and associated with the poor overall survival of neuroblastoma [
15]. Moreover,
MYCN amplification induced
E2F5 expression to promote neuroblastoma progression through regulation of cell cycle pathway [
16]. The higher expression levels of
E2F3 were associated with the worse clinical outcomes of stage 4S neuroblastoma [
17]. However, the expressions and prognosis of E2F transcription factors have not yet been studied in neuroblastoma in a comprehensive manner.
With the improvements of technologies, multiple neuroblastoma cohorts had been studied in transcriptional analysis. Previously, using integrated neuroblastoma cohorts, we analyzed the transcriptional features of neuroblastoma associated with
MYCN amplification [
15] and age at diagnosis ≥ 18 months [
18]. Here, using four independent neuroblastoma cohorts from Therapeutically Applicable Research to Generate Effective Treatments (TARGET), Gene Expression Omnibus (GEO) and European ArrayExpres datasets, we studied the expressions and prognosis of E2F transcription factors in neuroblastoma. Our results suggested that
E2F1 and
E2F3 were prognostic makers of neuroblastoma independent of
MYCN amplification and age of diagnosis.
Materials and methods
Data collection and processing
TARGET datasets were downloaded from
https://ocg.cancer.gov/ [
19]. GSE16476 [
20‐
22] and GSE85047 [
23] was downloaded from the GEO website (
www.ncbi.nlm.nih.gov/geo). E-MTAB-1781 dataset [
24] were downloaded from
https://www.ebi.ac.uk/arrayexpress/ website. Neuroblastoma patients with event free survival or overall survival were selected for further studies. All the data was processed using R software. The matrix files were annotated by corresponding platforms. Averaged expression values of same gene symbol were used for further studies.
Single sample Gene Set Enrichment Analysis (ssGSEA)
E2F_01 associated gene set (derived from c3.tft.v7.2.symbols gene sets) and cell cycle signaling pathway associated gene set (derived from c2.cp.kegg.v7.2.symbols gene sets) were downloaded from GSEA website. The scores of E2F_01 and cell cycle signaling pathway were determined using ssGSEA assay “GSVA” package in R software [
25]. “GSVA” is a non-parametric, unsupervised method for estimating variation of gene set enrichment through the expression dataset, thereby evaluating the scores of pathways or transcription factors in each sample.
Univariable and multivariable cox regression analysis
Univariable and multivariable cox regression analysis were carried out using “survival” and “survminer” packages “coxph” method in R software. The forest plots were generated using “forestplot” and “ggforest” packages in R software. The hazard ratio (HR) and P values were determined using cox regression survival analysis.
Kaplan–Meier survival analysis
Kaplan–Meier plots were created using “survival” and “survminer” packages in R software. Pediatric neuroblastoma patients were divided into “high” or “low” sub-groups based on the optimal cutoff points using “survminer” package “surv_cutpoint” method. “surv_cutpoint” is an outcome-oriented method by calculating all possible cutoff values between the lower and upper sub-groups. Cutoff points that were most significantly associated with the clinical outcomes were selected as the best cutoff points. P values were determined using log-rank test in “survival” package.
Risk score plot
The risk score plots were generated using “ggrisk” and “rms” packages in R software. The risk score was calculated based on the cox regression in “survival” package by summing up the variables in the cox model weighted by the corresponding regression coefficients. The cutoff was determined by the median of the risk score.
Heatmap presentation
The expression features of the genes associated with E2F1 and E2F3 expressions in neuroblastoma were clustered using “pheatmap” package in R software. The “average” method and “correlation” were used to determine the clustering scale and clustering distance, respectively.
Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway and transcription factor enrichment analysis
E2F1 and
E2F3 associated KEGG signaling pathways and transcription factors were determined using The Database for Annotation, Visualization and Integrated Discovery (DAVID) website (
https://david.ncifcrf.gov) [
27,
28]. The False Discovery Rate was used for multiple hypothesis testing. Signaling pathways or transcription factors with
P values < 0.05 were significantly enriched.
Statistical analysis
Statistical analysis was performed using the two tails paired student’s t test in GraphPad Prism software. P value < 0.05 was chosen to be significantly different.
Discussion
Neuroblastoma is a heterogeneous disease, and studies from distinct neuroblastoma cohorts may result different conclusions. Integrated analysis across different neuroblastoma cohorts may provide more robust and consistent results. In this study, using four neuroblastoma datasets, we showed that, compared with other E2F factors, E2F1 and E2F3 shared similar expressions and prognosis in pediatric neuroblastoma. Higher expression levels of E2F1 or E2F3 were associated with the worse prognosis of neuroblastoma. Also, E2F1 and E2F3 were associated with MYCN amplification and age of neuroblastoma diagnosis. Moreover, E2F1 and E2F3 were independent neuroblastoma prognostic factors. Combinations of E2F1 expression, E2F3 expression with MYCN amplification or age of diagnosis achieved better prognosis in neuroblastoma. On the contrary, other E2F transcription factors had no similar prognostic effect in neuroblastoma.
Furthermore,
E2F1 and
E2F3 also shared similar downstream transcriptional features in pediatric neuroblastoma. Genes associated with
E2F1 high expressions were also associated with
E2F3 high expressions. Those results suggested the functional redundancy of
E2F1 and
E2F3 in the regulation of neuroblastoma development. Similarly, in liver cancer,
E2F1 and
E2F3 were over-expressed, and synergistically induced the spontaneous development of liver cancer in mice [
30,
31]. In the further development of therapeutic strategies of neuroblastoma by targeting on E2F transcription factors,
E2F1 and
E2F3 should be simultaneously inhibited to achieve better clinical outcomes.
Uncontrolled cell cycle and DNA replication are important hallmarks of cancer [
32,
33]. The abnormal expressions of
E2F1 and
E2F3 may confer the uncontrolled cell cycle and DNA replication in neuroblastoma [
34]. We also showed that higher scores cell cycle signaling pathway were associated with the worse overall survival of neuroblastoma.
MYCN amplification mediated the reprogramming of metabolism is another hallmark of neuroblastoma [
35‐
37]. Our results suggested that genes associated with
E2F1 and
E2F3 expressions were enriched in pyrimidine metabolism and purine metabolism signaling pathways. Previously, our results had shown that pyrimidine metabolism signaling pathway was associated with the progression of LUAD [
38]. It was interesting to determine the prognosis of pyrimidine and purine metabolism signaling pathways in neuroblastoma. Moreover, the mechanisms of
E2F1 and
E2F3 in the regulations of metabolism signaling pathways in neuroblastoma should be further studied.
Except
MYCN amplification and age of diagnosis, more reliable prognostic markers are required to predict the clinical outcomes of neuroblastoma [
39]. In this paper, we highlighted that
E2F1 and
E2F3 were prognostic makers of neuroblastoma independent of
MYCN amplification and age of diagnosis. The prognostic effects of
E2F1 and
E2F3 may relate to the regulations of cell cycle signaling pathway and metabolism signaling pathways. However, those conclusions were generated from published datasets and lack of validations using further experiments. Therefore, functions of
E2F1 and
E2F3 should be further studied in neuroblastoma cells by using in vivo or in vitro experiments.
Conclusions
E2F1 and E2F3 were prognostic factors in neuroblastoma, independent of MYCN amplification and age of diagnosis. The expression levels of E2F1 and E2F3 were higher in neuroblastoma patients with MYCN amplification or age at diagnosis ≥ 18 months. Combinations of E2F1 expression, E2F3 expression with MYCN amplification or age of diagnosis achieved better prognosis of neuroblastoma. Genes associated with E2F1 and E2F3 expressions in neuroblastoma were significantly enriched in cell cycle signaling pathway. Higher scores of cell cycle signaling pathway were correlated with the adverse prognosis of neuroblastoma.
Open AccessThis article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit
http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (
http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.