Early 2-Factor Transcription Factors Associated with Progression and Recurrence in Bevacizumab-Responsive Subtypes of Glioblastoma

1. Introduction

The early 2-factor (E2F) family consists of eight key members in humans; E2F1-8, these members have highly similar DNA-binding domains that directly interact with consensus sequences [1]. The members are categorized into three subgroups based on sequence similarity and functional activity: activator proteins (E2F1-3), atypical repressors (E2F7-8), and canonical repressors (E2F4-6). Each member has distinct expression and functional patterns that align with their subgroup classification [2,3]. E2Fs play a crucial role in the cell cycle by forming the core transcriptional axis. In addition to cell cycle regulation, E2Fs also influence various biological pathways that contribute to malignant progression, such as apoptosis, angiogenesis, and metabolism in several cancers including glioblastoma [4,5,6].Growing evidence indicates that elevated expression of E2F family members is significantly associated with malignant progression in multiple tumor types [7,8,9]. Notably, E2F1 and E2F2 are associated with increased expression of pro-angiogenic genes in breast cancer, potentially contributing to a more aggressive phenotype [10].Recently, the role of E2F8 in cancer has attracted attention due to its aberrant expression, which is associated with angiogenesis and cancer cell proliferation, leading to adverse clinical outcomes in most human cancers [11]. Therefore, E2F members may serve as promising therapeutic targets and potential biomarkers for specific cancers.

Recent research, including our own, has used miRNA biomarkers and machine learning techniques to classify different subtypes of GBM that respond to Bevacizumab (BVZ) treatment [12]. BVZ therapy targets vascular endothelial growth factor (VEGF) and inhibits angiogenesis within the tumor. Current opinion suggests that BVZ therapy may prolong progression-free survival (PFS) in GBM treatment, but there is no evidence of its effect on overall survival (OS). However, our observation indicated that the BVZ-responsive subtype of GBM tends to exhibit greater malignancy compared to non-responsive subtypes, resulting in a shorter overall survival of two or three months. It is possible that the effect of BVZ treatmenton is masked by this increased malignancy. Furthermore, some patients with genetic alterations, such as EGFR vIII-positive GBM show better outcomes with BVZ treatment [13]. Therefore, it is important to accurately classify, treat, and investigate GBM patients.

Despite responding to treatment, including BVZ therapy, most GBM cases eventually relapse. The underlying mechanisms behind this relapse have been extensively studied. After BVZ treatment, glioblastomas can stimulate the growth of tumor vessels by expressing various angiogenic factors such as bFGF (Fibroblast Growth Factor 2) [14]. In addition, research has shown that post-treatment recurrent GBM is also driven by pre-existing, treatment-resistant stem-like cells present in the tumor microenvironment [15]. Moreover, genetic reprogramming and differential oncogene activation, and regulation of self-renewal, tumorigenicity, and metastasis in the relapsed tumor cells and their microenvironments contribute to GBM recurrence [16].By studying these aspects of relapse mechanisms, we can gain a better understanding of GBM development and progression, leading to the identification of new targets and effective treatments.

Previously, our research identified multiple transcription factors, including members of the E2F family, involved in GBM treated with BVZ [17]. Despite the impact of BVZ on patient prognosis, overall survival (OS) remains unaffected, and tumor recurrence continues to be a major challenge in GBM treatment. Therefore, this study aims to investigate the influence of E2F family members on the malignancy of BVZ-responsive GBM subtypes before and after treatment. The analysis will be based on the clinicopathological characteristics of the patients. The goal of this study is to gain insights into the molecular mechanisms underlying recurrence and progression following BVZ treatment in GBM patients, and to provide potential guidance for future therapeutic and translational research strategies.

2. Materials and Methods

2.1. Data Collection

In this study, mRNA and microRNA (miRNA) expression profiles, along with clinico-pathological data of glioblastoma (GBM), were sourced from the Cancer Genome Atlas (TCGA) pilot study. A total of 426 GBM samples were acquired from the TCGA data portal (http://cancergenome.nih.gov/dataportal/) [18]. Additionally, data on GBM patients treated with Bevacizumab (BVZ) were downloaded from the Gene Expression Omnibus (GEO), dataset GSE79671. Sequencing data from GBM tissues, collected before and after BVZ treatment, were refined using Trimmomatic and subsequently mapped to the human genome (hg19) [19,20]. These gene expressions were quantified as fragments per kilobase million (FPKM). For analytical purposes, data from the same patients, pre- and post-BVZ treatment, were paired to facilitate a comparative statistical analysis.

2.2. Bioinformatics Analysis

To explore the potential functions of the eight E2Fs in GBM, several bioinformatics tools were applied to determine the associations between E2Fs expression patterns and the clinical characteristics of samples from the TCGA and GSE79671 datasets. A protein-protein interaction (PPI) network of the eight E2Fs and other differentially expressed (DE) genes was constructed using the STRING database, Cytoscape software, and Bioinformatics [17,21]. For E2Fs expression in GBM patients and normal individuals, DoSurvive and its specimens were used [22]. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed by DE genes between pre- and post-BVZ treatment for BVZ responsive subtypes of GBM.

2.3. Software and Statistical Methods

In this study, the software and online tools used included R-project (version 4.0.4, www.r-project.org), MatLab2023a (https://www.mathworks.com), MeV (version 4.9.0, MEV, LLC., Walnut Creek, CA, http://www.tm4.org/mev.html), g:Profiler (ELIXIR, Tartu, Estonia, https://biit.cs.ut.ee/gprofler/), Bioinformatics (https://bioinformatics.com.cn), and Prism 9.0 (GraphPad Software, Inc). The statistical methods used included Student’s t-test, the standard Bonferroni adjusted t-test, paired t-test, and the Wilcoxon-Mann–Whitney test if the data distribution was not normal.

2.4. Variable Selection and TreeBagger Model Construction

We employed TreeBagger (TB), one of the random forest algorithms, for further prediction and assessment of patient survival after BVZ treatment for GBM patients[23].The TCGA dataset has been classified by miRNA biomarkers and SVM algorithm, and the GSE79671 dataset has been classified by CT scan before and after BVZ treatment for GBM patients.E2F1-8 expression data were extracted from TCGA and GSE79671 and classified into BVZ responsive subtype and the non-responsive subtype of GBM, and the pre- and post-BVZ treatment for the BVZ responsive subtypes of GBM.

Based on the clinical survival data, the survival criteria were classified as:

Criteria = [1, OS≤450 days;2, OS>450 days].

All parameters in TB were optimized, including tree number, leaf splits, and number to even. One process of selecting tree numbers has been shown in the results section. After the TB model has been built, the TCGA dataset was used as the training set and GSE79671 was used as the testing and predicted set. In these processes, Pearson’s correlation was also used to select variables.

3. Results

3.1. Differential Transcription Factor Expression in Bevacizumab-Responsive Glioblastoma Post-BVZ Treatment

To explore the impact of bevacizumab (BVZ) treatment on BVZ-responsive glioblastoma (GBM) subtypes, we performed a comprehensive analysis of functional assays including Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and transcription factors (TF).Our analysis involved comparing differentially expressed (DE) genes obtained from 426 TCGA datasets between BVZ-responsive subtypes, identified using miRNA biomarkers and machine learning approaches [12], whereas DE genes from 17 GBM patients using CT scan before and after BVZ treatment to identify BVZ-responsive subtypes[20]. While typical changes were observed in GO and KEGG pathways, striking alterations were detected in TFs expression, particularly evident in Figure 1A.

When analyzing the 557 DE genes from GSE79671(p=0.05 for paired t-test), 106 TFs surpassed the significance threshold (p=0.01), as illustrated in Figure 1A. Notably, only 36 TFs met the threshold using the 1000DE genes from the TCGA dataset (Figure 1B). Conversely, no significant TFs were identified between pre- and post-BVZ treatment in BVZ non-responsive patients using the 326 DE genes (p=0.05 for paired t-test), not even other terms except one (Figure C), indicating non-specific treatment. All three analyses were set at the same parameters for those functional analyses and rechecked them on the same day (5-8-2024). These substantial differences underscore the importance of further investigation into TF expression dynamics, especially in the BVZ responsive subtype of GBM.

Furthermore, the protein-protein interaction (PPI) network of these 106 TFs, including 7 of E2F1-8, namely E2F1, E2F2, E2F3, E2F4, E2F6, E2F7, and E2F8, revealed these E2Fs as hub genes, indicating that they play a key role in malignant cells. Supplementary Figure 1 shows the interaction of these TFs in BVZ-responsive GBM patients treated with BVZ. Details of significant TFs, including names, sequences, and p-values, are provided in Supplementary Table S1.

3.2. E2F Expression and Bevacizumab-Responsive Subtypes of Glioblastoma

To explore the potential involvement of E2Fs in the malignancy of BVZ-responsive GBM subtypes, we conducted a systematic analysis of the expression patterns of the eight E2F transcription factors and their associations with various clinicopathological features before and after BVZ treatment. Our analysis encompassed 426 GBM patients from the TCGA dataset, classified using miRNA biomarkers and a machine learning algorithm, and 17GBM patients from the GSE79671 dataset, classified into BVZ-responsive subtypes based on CT scans conducted pre- and post-BVZ treatment.

All eight E2Fs exhibited significantly altered expression in GBM compared to normal controls, indicating their crucial association with GBM malignancy. Specifically, we observed significant upregulation of all eight E2Fs in GBM (Table 1), based on GBM and normal brain tissue specimens from the DoSurvive dataset. The log2(FC(GBM/Normal)) values of E2F2, E2F5, E2F7, and E2F8 are 5.9, 2.4, 5.5, and 6.9, respectively.

Using TAGA dataset, patients’ risk scores were calculated based on E2F1-8 expression, OS, and BVZ response subtypes of GBM. The risk score is shown in Figure 2A, with high and low risks represented in pink and blue. The patients’ OS is shown in Figure 2B, and the expression of E2F1-8 is shown in Figure 2C. There were no significant differences between BVZ-responsive and non-responsive GBM subtypes in the expression of E2F1-8. The hazards of E2F1-8 were calculated using cox regression method, and there is no significant difference, as shown in Figure 2D.

Similar to the TCGA dataset, there was no significant difference between BVZ-responsive and non-responsive GBM subtypes in GSE79671 before treatment, as shown in Figure 2E. However, post-BVZ treatment analysis in the GSE79671 dataset of BVZ-responsive patients showed increased expression of E2F2 and E2F8, with particular note being given to E2F8 (p-value = 0.006) and E2F2 (p-value = 0.052) (Figure 2F). In contrast, despite the presence of other DE genes, there were no significant changes before and after BVZ treatment in BVZ-nonresponsive patients. These findings highlight the complex involvement of E2F in BVZ treatment response and the underlying biological mechanisms of GBM recurrence.

3.3. Predicting Survival in Bevacizumab-Responsive Subtypes of Glioblastoma Using TreeBagger Analysis

Machine learning algorithms have become widely utilized for prediction and prognostic purposes in cancer research. Previously, we employed Support Vector Machine (SVM) and miRNA biomarkers to classify BVZ-responsive subtypes of GBM. However, SVM is typically suited for scenarios involving two or three variables. Our research aimed to predict the survival terms of patients with glioblastoma (GBM) who are responsive or nonresponsive to bevacizumab (BVZ) treatment using an advanced machine learning algorithm. In this study, we chose TreeBagger, a Random Forest ensemble learning method, for its ability to manage a high-dimensional dataset and to model complex, nonlinear interactions.

The study incorporated all eight members of the E2F transcription factor family as predictive variables within the TreeBagger model. Despite decades of efforts to optimize and combine glioblastoma treatments [24], outcomes for patients with GBM remain poor, with a median life expectancy of approximately 15 months after diagnosis [25]. In this study, we set the threshold as 15 months (450 days) to predict patient survival terms based on E2Fs expression and the machine learning analysis. Using the TCGA dataset, 80% and 20% of the dataset were set as the training set and test set respectively. The accuracy for predicting patient survival term reached 75%. The confusion matrix illustrating the performance is presented in Figure 3A. The area under the curve (AUC) was determined to be 60%, as depicted in Figure 3B. For the calculation, the optimization of tree number is shown in Figure 3C. Employing 5-fold cross-validation yielded an average accuracy of 67%, indicating no overfitting.

To validate our model, we utilized the GSE79671 dataset and applied our trained TreeBagger model to forecast survival terms before and after BVZ treatment. Prediction of patient survival before and after BVZ treatment using the GSE79671 dataset and trained TreeBagger by the TCGA dataset. When comparing the prediction of BVZ-responsive versus non-responsive GBM subtypes, this analysis did not show a significant improvement in survival after BVZ treatment, consistent with clinical observations. However, due to the small sample size, these findings necessitate further investigation.

Additionally, we explored the predictive potential of TreeBagger using a subset (E2F1, E2F2, E2F3, E2F7, and E2F8) of E2F members obtained from Pearson's correlation analysis as shown in Figure 3D. Specifically, we investigated the performance of TreeBagger with 150 trees and 25 leaf splits to predict patient survival using five E2F members. The results closely paralleled those obtained with all eight E2F members, with the highest accuracy reaching 73.8% and an AUC of 50%.

Considering the involvement of 106 transcription factors, these findings underscore the significant role of the E2F family in tumor progression and patient survival. Notably, based on the confusion matrix, E2Fs and TreeBagger algorithms showed better predictive ability in identifying survival probabilities of BVZ-responsive versus BVZ-non-responsive GBM subtypes because, as mentioned earlier[12], compared with BVZ-non-responsive patients, these patients have lower OS.

3.4. Post-Treatment Functional Analysis of E2Fs on GBM BVZ Response Subtypes

GBM frequently relapses after BVZ treatment, even in patients who respond to BVZ, and experience prolongation of progression-free survival (PFS), but do not show a significant prolongation of overall survival OS [13,26]. However, there are some conflicting reports. To explore possible mechanisms of GBM recurrence in BVZ-responsive subtypes, the functional regulation of E2F2 and E2F8 after BVZ treatment was analyzed.

We identified genes that were significantly differentially expressed before and after BVZ treatment in the BVZ-responsive subtype of GBM. We further analyzed the biological processes and cellular components associated with malignancy, as well as the genes associated with enriched miRNA in this process. In Figure 4A, the names of functional processes, gene counts, and p-values are shown. We observed enrichment of genes involved in sprouting angiogenesis (GO:0002040), cell cycle comprising mitosis without cytokinesis (GO:0033301), RNA polymerase II transcription regulator complex (GO:0090575), and miRNAlet-7a. These GO terms and pathway analyses further revealed that important genes involved in angiogenesis, cell cycle mitosis, and cell proliferation were significantly enriched in the BVZ-responsive subtype of GBM after BVZ treatment.

Furthermore, the regulatory gene network of E2F8 was explored by integrating DE gene expression in the BVZ responsive subtype of GBM before and after treatment into the current BioGrid network, as shown in Figure 4B. The protein-protein interaction (PPI) network data were organized and visualized using Cytoscape and 25 nodes and 78 edges were found in the PPI network (All nodes are shown in Supplementary Table S2). The size and color depth of the nodes in the network correspond to their degree labeled by logFC values. The expression of STUB1 (STIP1 Homology And U-Box Containing Protein 1, E3 Ubiquitin) and YWHAQ (Tyrosine 3-Monooxygenase/Tryptophan 5-Monooxygenase) increases (appears red) after BVZ treatment, indicating their involvement in the E2F8 network and potential regulation by E2F8. Other nodes do not appear in the DE gene list, are shown in light blue. The colors of STUB1 and YWHAQ are shown in red, in this network, indicating that their expression levels increase and overlap with DE genes of the BVZ-responsive subtype of GBM after BVZ treatment.

Because of the involvement of angiogenesis GO:0002040, the expression of VEGF in BVZ-responsive GBM subtypes after BVZ treatment was examined. There were 17 paired datasets from a total of 21 GBM samples before and after BVZ treatment, including 11 BVZ non-responsive data (Non) and 6 BVZ responsive data (BVZ). The log(FC) of (VEGFpost / VEGFpre) was calculated and presented as a bidirectional bar graph. To our surprise, 80% of patients who did not respond to BVZ treatment had increased or no change in VEGF expression after treatment, as shown in Figure 4C, while 67% of GBM patients who responded to BVZ treatment had increased VEGF expression, as shown in Figure 4D. This fact suggests that BVZ treatment in unresponsive GBM patients may worsen disease progression, and accurate BVZ treatment in GBM patients is strongly recommended.

4. Discussion and Conclusions

This study focused on the E2Fs family to explore potential relapse mechanisms in patients with GBM who experience remission after BVZ treatment. BVZ treatment has been shown to have significant impact on transcription factors and their hub genes, including E2Fs, which promoted our investigation. By analyzing clinicopathological characteristics and utilizing machine learning bioinformatics approaches, we gained insights and identified targets related to angiogenesis and proliferation through E2F8 regulation after BVZ treatment. Most importantly, BVZ therapy may worsen disease progression in patients with unresponsive GBM because of increased VEGF expression. Additionally, we have identified the involvement of STUB1 and YWHAQ, suggesting that their potential roles in GBM recurrence and progression in patients who are responsive to BVZ treatment. These findings may guide future therapeutic strategies. However, further experimentation is needed to confirm these results derived from bioinformatics and artificial intelligence analyses.

Previous studies have indicated that E2F7 and E2F8, both considered as atypical E2Fs, function as transcriptional regulators of VEGF expression [27]. Specifically, they directly bind to and stimulate the VEGF promoter, independent of canonical E2F binding elements, thereby controlling angiogenesis. BVZ, an antibody that neutralizes and inhibits VEGF, thus suppressing angiogenesis, has been found to potentially increase E2F8 expression. This increase in E2F8 expression may reduce the effectiveness of BVZ treatment, rapidly restore VEGF expression and angiogenesis, and increase its hazard of risk, which is lower than that of other E2Fs before treatment. However, the underlying mechanisms responsible for BVZ enhancing E2F8 expression remain unknown. Furthermore, studies on mouse and zebrafish with non-functional E2F7/8 mutant have demonstrated varying degree of defects in angiogenesis and apoptosis, highlighting the crucial role of E2F7/8 in angiogenesis. In mouse model, the inhibition of E2F8 has been shown to suppress the development of various tumors by inhibiting angiogenesis [28]. Therefore, incorporating E2F8 inhibition into BVZ therapy may potentially enhance the inhibition of angiogenesis and lead to improved treatment efficiency.

E2F8 plays a significant role in cancer by affecting the proliferation and differentiation of cancer stem cells (CSCs). CSCs are a small group of cancer cells that can self-renew and give rise to a diverse lineage of tumor-forming cells. They are largely responsible for tumor growth, progression, and recurrence [29]. Studies have shown that E2F8 is upregulated in various cancers, including glioblastoma, ovarian cancer, and hepatocellular carcinoma, where it promotes cell proliferation and other effects [11,30]. Moreover, the regulation of E2F8 throughout the cell cycle demonstrates extensive coordination between phosphorylation, ubiquitination, and transcription in the mammalian cell cycle. E2F8 is involved in the regulation of cyclin D1 and promotes the accumulation of S-phase cells in hepatocellular carcinoma, further highlighting its role in cell proliferation [31].Additionally, E2F8 is one of the E2F family of transcription factors and is crucial for various cellular processes, including those associated with cancer stem cells [3]. Therefore, targeting E2F8may present a promising therapeutic approach to disrupt cancer stem cells and hinder tumor growth and progression.

MicroRNAs, as post-transcriptional regulators, play a vital role in the initiation, recurrence, and progression of GBM [32,33,34], particularly let-7a, which regulates cancer stemness [35]. In this study, we found that let-7a is involved in the E2F8 regulated network after BVZ treatment in its responded patients. Generally, highly expressed miRNAs act as negative regulators; however, it has been reported thatlet-7 expression levels were significantly reduced in tumor initiating cells (TICs) but increased in differentiated tumor cells in breast cancer [36]. In tumor cells, Inhibiting let-7 by antisense oligonucleotides enhanced the self-renewal capacity of the cells, indicating the role of let-7 in suppressing tumorigenicity by downregulating CSC properties. On the other hand, CSCs have a strong ability to adapt their metabolism to survive and proliferate under adverse environmental conditions like hypoxia and acidosis because cancer cells heavily rely on glycolytic fermentation for energy metabolism even when oxygen is available. Related studies have shown that the let-7/Lin28 axis regulates cancer progression and CSC activity through aerobic glycolysis in hepatocellular carcinoma cells (HCC) with high Lin28 and low let-7 expression [37].This suggests that the miRNA let-7 may also be a target to block or inhibit GBM recurrence by inhibiting CSC growth.

Furthermore, STUB1, an E3 ligase, plays a vital role in cancer by regulating multiple pathways and participating in the E2F8 regulatory network. The study reveals that STUB1 inhibitors suppress interferon gamma (IFNγ) responses by degrading IFNγ receptor 1 (IFNGR1). As a result, STUB1 inhibitors enhance tumor response to checkpoint immunotherapy[38]. Additional studies demonstrate that in lung cancer, knocking down STUB1 increases the likelihood of metastasis to the lungs of mice when injected intravenously or subcutaneously. In breast and lung cancer, STUB1 interacts with various proteins through its ubiquitination-dependent proteasome activity. It acts as a negative regulator associated with different proteins, affecting cell cycle progression and metastasis [39,40].

YWHAQ belongs to the 14-3-3 (YWHA) proteins, which are phospho-serine/-threonine binding proteins that are highly conserved and participate in many important cellular processes, including metabolism, protein transport, signal transduction, apoptosis and cell cycle regulation. YWHAQ is upregulated in patients with amyotrophic lateral sclerosis (ALS) [41]. While its expression in cancer has not been extensively reported, other members of the family have been studied in cancer. For example, YWHAZ, also known as 14-3-3ζ, is frequently upregulated in various cancers and acts as an oncogene, promoting tumor progression through multiple cellular activities. In breast cancer tissues, the expression levels of DAAM1(Dishevelled-associated activator of morphogenesis 1) and YWHAZ were significantly upregulated and colocalized, with YWHAZ binding to DAAM at its phosphorylated site regulating cancer cell migration and strongly associated with poor prognosis [42]. The World Health Organization (WHO) now recognizes four categories of endometrial stromal tumor, with YWHAE, another 14-3-3 protein, translocation identifying high-grade endometrial stromal sarcoma (HG-ESSs). This highlights the impact of the YWHA family in molecular genetics in cancer and supports for a new classification system [43].

Accurate BVZ treatment is highly recommended. Our previous study showed that BVZ treatment in its unresponsive patients of GBM resulted in multiple side effects, including aging [17]. In this study, we found that BVZ treatment is most likely to increase VEGF expression, particularly in patients with BVZ unresponsive GBM. This increase may lead to GBM recurrence and affect the OS of GBM patients treated with BVZ.

In conclusion, a closer examination of E2F members provides new insights into the mechanisms of glioblastoma recurrence after bevacizumab (BVZ) treatment, specifically associated with E2F8. However, there is still much more to explore. It is important to note that although the E2F family positively influences angiogenesis, proliferation, and cancer progression, they also regulate some negative regulators or tumor suppressors. This results in complex regulation of tumor growth and diverse responses to cancer treatments. Therefore, targeting hub genes like E2F8 or multiple targets simultaneously may yield favorable outcomes in cancer treatment.