Treatment-Altered VEGF-Associated Network in GBM BVZ-Responsive Subtypes: Bioinformatics Case Studies

Objective: As an anti-angiogenetic treatment, bevacizumab (BVZ) therapy for GBM patients targets vascular endothelial growth factor A (VEGF), which can significantly reduce tumor sizes in the early stage based on imaging studies. However, the related mechanism of the VEGF-associated network has not yet been elucidated. Recently, subtypes of GBM BVZ response were classified and detected, so the treatment effects and related mechanisms could be reassessed individually for these different subtypes. In this study, we want to further explore possible mechanisms, effects


Introduction
Glioblastoma multiforme (GBM) is the deadliest brain tumor with a median lifespan of approximately 15 months after diagnosis [1]. As GBM is a heterogeneous and multivariate disease, despite decades of investigation into the pathogenesis of GBM, the various genetic factors involved in the development and recurrence of the disease and their regulatory pathways remain poorly understood [2][3][4][5]. At present, there are many treatment methods for GBM patients, including resection, radiation, chemotherapy, immunotherapy, and antiangiogenic therapy, such as bevacizumab (BVZ) treatment (trade name Avastin or BEV). BVZ treatments target vascular endothelial growth factor A (VEGF), which can significantly reduce the early stage tumor diameters of some GBM patients based on imaging studies [6,7].However, BVZ alone or in combination with treatment had a significant effect on the transcriptional changes and tumor size decreases in its responders, with approximately 30-33% of all GBM patients; but changes in nonresponders were small [8][9][10].
Angiogenesis, as a VEGF-related process, plays an important role in the development and recurrence of GBM, and despite numerous reports,it remains unclear in these processes [11]. It requires advanced and precise investigation. Recently, the classification and detection of GBM BVZ-responsive subtypes using machine learning approaches and miRNA biomarkers, including miR-21, miR-10b, and miR-197, allowed further and preciseinvestigation of functional pathways related to VEGF and angiogenesis. Compared with patients with the GBM BVZnonresponsive subtype, patients with the GBM BVZ-responsive subtype had shorter survival times, suggesting that the effect of BVZ may be better than the previously calculated; and lower levels of VEGF methylation, indicating its high expression levels and the high rate of microvascular cell proliferation. Furthermore, differentially expressed (DE) mRNAs and different functional pathways were found between these two subtypes [10,12]. Therefore, it is worth re-evaluating and precisely exploring existing datasets in light of new classifications of regulatory angiogenesis mechanisms.
To explore the mechanisms of various genes and the VEGFassociated network in GBM, on the basis of the classification of GBM BVZ subtypes and the comparison of before and after BVZ treatment for GBM patients, multiple bioinformatics and statistical methods and mRNA expression datasets were used to analyze the functional pathways and network between these subtypes, as well as DE mRNAs.

Experimental datasets
The datasets (GSE79671) of GBM patients with BVZ treatment were downloaded from from gene expression omnibus GEO. A sequencing dataset from GBM tissue before and after BVZ treatment was trimmed with Trimmomatic and mapped to the human genome (hg19) [13]. After processing, 15630 genes were obtained and expressed as fragments per kilobase per million reads (FPKM) [9]. The data of the same patient before and after BVZ treatment were paired for statistical study.

Patients' characteristics
This study analyzed BVZ-related network and DE gene expression using 426 GBM patient profiles from which references with complete gene profiles from microarray or high through the sequence. In Table 1, the numbers and percentages of subtypes of patients who are classified in BVZ-responsive or nonresponsive subtypes are shown, along with the classification methods used. a: These patients have gene expression profiles from both before and after BVZ treatment. b: These patients have complete miRNA and mRNA expression profiles for machine learning (ML) analysis.

Different gene expression patterns between GBM BVZ subtypes before and after BVZ treatments
Gene expression patterns were significantly different in GBM BVZ-responsive and nonresponsive patients before and after BVZ treatment, as shown in Figure 1. The treatment response was assessed according to the RANO criteria and confirmed on the subsequent follow-up MRI [9], including complete or partial response (CR + PR).After dividing GBM patients into GBM BVZ-responsive and nonresponsive subtypes, paired t-tests before and after BVZ treatment were performed for these two subtypes. For GBM BVZresponsive patients following BVZ treatment, 469 genes were significantly expressed based on apaired t-test (p<0.05), in which 391 gene expression levels were increased, 78 gene expression levels were decreased, and their expression histogramsare shown in Figure 1A and B. In the graph, major genes are located in the nearby regions, indicating specific expression patterns. Alternatively, for GBM BVZnonresponsive patients following BVZ treatment, 169 genes were significantly expressed, 60 gene expression levels were increased, and 109 gene expression levels were decreased; their expression histogramsare shown in Figure1B, where there are no specific expression regions like in Figure 1A. For gene ontology (GO) and KEGG pathway analysis of 469 DE mRNAs from GBM BVZ-responsive patients, 84 transcription factors (TF), 10 biological processes (BP), 5 cellular components (CC), and 3 KEGG pathways crossed the threshold (p=0.01) and were involved in BVZ treatment for GBM BVZ-responsive patients, as shown in Figure 1C. The VEGF pathway exists as one of the KEGG pathways involved. The most involved part was TF, which may be an unexpected result. In contrast, using the same analysis procedure for 169 DE mRNAs from GBM BVZ-nonresponsive patients, no functional pathways crossed the threshold (p=0.05) and were involved in BVZ treatment, suggesting that DE mRNAs expression was not specific.
This result indicated that after BVZ treatment, the gene expression patterns and related functional pathways in GBM BVZresponsive patients were specific and beneficial to the patients; in contract, in the treatment for GBM BVZ-nonresponsive patients, the gene expression patterns and related functional pathways were not specific and may not be beneficial to the patients, suggesting the need for prescreened for BVZ therapy.

VEGF-associated and -unassociated angiogenesis pathways
Angiogenesis plays a critical role in the development and recurrence of GBM with VEGF-associated and unassociated pathways. According to Venn diagram analysis, for GBM BVZresponsive patients, angiopoietin1 (Angpt1) and Sh2 domain containing 2A (SH2D2A) existed among angiogenesis (315 genes), VEGF (101 genes), and the differential expression (DE) mRNA profile (469 genes) of BVZ treatment, as shown in Figure  2A, indicating that they are VEGF-associated genes. In the graph, STRING provided 315 angiogenesis genes and 101 ranked with high confidence VEGF-associated genes, and 5 differentially expressed genes were obtained by paired t-test between pre-and post-BVZ treatment for GBM BVZ-responsive subtypes. The angiogenesis pathway existed among these genes, these two genes connected directly with VEGF with thick lines suggested confident not predicted connections, as shown in Figure 2B. Angpt1 has been reported for its complex angiogenesis function [20], but there is almost no report about the angiogenesis function of SH2D2A although it was decreased seven times (p=0.04, n=5) after BVZ treatment.
Lymphoid enhancer-binding factor 1 (LEF1), phosducin-like3 (PDCL3), and a thrombospondin type 1 domain containing 7A (THSD7A) are angiogenesis genes but not in the VEGFassociated gene group, so they belong to the VEGF-unassociated angiogenesis pathway. As shown in Figure 2B, they are not directly connected with VEGF but through other genes, THSD7A was not even found to connect with any other VEGF-related genes in this network. There were no reports about the angiogenesis function of three genes except LEF1, which was decreased approximately 1.5 times (p=0.037, n=5) after BVZ treatment. A recent study in zebrafish showed that the transcription factor LEF1 positively regulated the specification of embryonic endothelial cells during embryonic cell growth through the Wnt pathway [21]. This result suggested that LEF1 may be positively correlated with the growth of vascular endothelial cells, so BVZ treatment suppressed the expression level of LEF1, indicating that it may be beneficial to inhibit the growth of tumor blood vessels. As shown in Figure 2B, LEF1 is not directly associated with VEGF, so its effects may be VEGF indirect. Comparing the expression levels of Angpt1, SH2D2A, and LEF1 between GBM BVZ-responsive and nonresponsive patients before treatment as shown in Figure 2C, there were no significant differences among these three genes. However, after treatment, their expression levels changed significantly, among which Angpt1 expression levels are shown in Figure 2D, suggesting these genes may serve as therapeutic indicators to monitor treatment.

Figure 3:
Aging and VEGF related network in GBM BVZ-nonresponsive patients. A) Venn diagram showing the relationship of agingassociated genes, VEGF-associated genes, and DE mRNAs before/after BVZ treatment in GBM BVZ-nonresponsive patients. B) Network between VEGF, aging-associated and DE mRNAs before/after BVZ treatment in GBM BVZ-nonresponsive patients.

Possible side effects for GBM BVZ-nonresponsive patients following BVZ treatments
The side effects may be caused when BVZ treats GBM BVZ-nonresponsive patients. Several genes were significantly altered in the pathways of apoptosis, inflammation, and cellular proliferation, but they were altered by no more than three genes in the aging pathway. After multiple assessments, significant decreases in Ephrin type A receptor1 (EPHA1), endothelial cellspecific molecular 1 (ESM1), and gremin 1 (GREM1) were found, as shown in Figure 3A. Compared with pretreatment, the expression levels of EPHA1, ESM1, and GREM1 were reduced by 5.25, 5.22, and 2.17 fold (p=0.004, 0.04, and 0.04, n=12) after BVZ treatment, respectively. They were predicted to link to VEGF, but there is no confidence for this connection, and are represented by thin lines rather than thick lines between genes, as shown in Figure 3B. However, these results require further experiments.

Discussion
According to statistical and bioinformatics analyses, this study explored and compared the effects of BVZ treatments for GBM BVZ-responsive and nonresponsive patients, revealing that the therapeutic effects of BVZ on GBM BVZ-responsive patients may go through VEGF-associated and non associated angiogenesis pathways, and BVZ treatment for GBM BVZ-nonresponsive patents may lead to aging-related side effects. Therefore, prescreening for BVZ therapy in GBM BVZ-responsive patients is recommended to avoid side effects.
The anti-angiogenesis effects of BVZ treatment on GBM BVZ-responsive patients may occur through the inhibition of Angpt1. Angpt1, one of four family members that bind to the Tie2 tyrosine kinase receptor with varying outcomes, is the main ligand for Tie2 or Tie-4 and an agonistic ligand [22,23]. In synergy with VEGF, Angpt1 had been shown to enhance angiogenesis in an animal aorta model and increased vessel density in a corneal implant experiment [24,25]. In this study, Angpt1 expression levels were significantly reduced after BVZ treatment in GBM BVZ-responsive patients but not in GBM BVZ-nonresponsive patients, suggesting the inhibition of tumor angiogenesis. These angiogenic effects of Angpt1 were also demonstrated by animals with its null mutations, as Angpt1 null embryos were unable to form complex vascular networks and exhibited reduced mural cell support for vessels [20]. However, Angpt1 has complex roles in tumor angiogenesis, and a growing number of reports indicate that Angpt1 stimulation of mural cells with endothelial cells may lead to the stabilization of newly formed blood vessels, which in turn may limit other continuous angiogenesis in tumors, thereby inhibiting tumor growth. Therefore, these Angpt1 double-effects may result in different effects at different times.
Aging-related side effects in GBM BVZ-nonresponsive patients may be caused by BVZ treatment, possibly through significant expression changes inEPHA1 and ESM1 based on current studies, suggesting the need of predetectionPre detection of GBM BVZ-responsive patients for BVZ treatment. An aged animal study showed that after 12 weeks of exercise, marked overexpression of ESM1 was observed in exercised animals compared with control animals, and significantly improved diastolic function was also observed in these animals [26]. After BVZ treatment for GBM BVZ-nonresponsive patients, the expression levels were dramatically decreased, possibly indicating the loss of diastolic function or one of the aging phenotypes. In addition, EphA1 belongs to the family of ephrin receptors that are involved in developmental events, especially in the development of the nervous system. Recent Drosophila studies showed that EphA1 is anAlzheimer's disease-associated gene, asRNAimediated knockdown of ephrin significantly impaired fly memory [27].
Given that patients who received preoperative and postoperative BVZ are known to be at risk of increased morbidity and mortality caused by side effects, patients who do not respond to BVZ treatment should not be taken at risk. From a practical standpoint, patients treated with BVZ are at risk of wound complications including dehiscence, CSF leak, infections, etc. Based on a study of 209 GBM patients, significant healing complications occurred in 44% of patients treated with preoperative BVZ compared to 9% of untreated patients, which significantly increased the rate of morbidity and mortality in this patient population [28]. Careful selection of patients for initiation of BVZ treatment is very important, taking into account laboratory and clinical side effects, and it should be done as soon as possible.

Conclusion
This study demonstrated the existence of VEGF-associated and nonassociated angiogenesis pathways during BVZ treatment for GBM BVZ-responsive patients and nonspecific gene expression patterns. Most importantly, aging-related side effects may exist in GBM BVZ-nonresponsive patients after BVZ treatment, in addition to the side effects of healing complications caused by BVZ treatment, thus prescreening GBM BVZ-responsive patients is strongly recommended for BVZ therapy to avoid side effects. Many factors exposed here are still unclear, so more experiments are needed to further explore the mechanism of angiogenesis in GBM.