The results from the randomized phase II BELOB trial provided evidence for a potential benefit of bevacizumab (beva), a humanized monoclonal antibody against circulating VEGF-A, when added to CCNU chemotherapy in patients with recurrent glioblastoma (GBM). In this study, we performed gene expression profiling (DASL and RNA-seq) of formalin-fixed, paraffin-embedded tumor material from participants of the BELOB trial to identify patients with recurrent GBM who benefitted most from beva+CCNU treatment. We demonstrate that tumors assigned to the IGS-18 or “classical” subtype and treated with beva+CCNU showed a significant benefit in progression-free survival and a trend toward benefit in overall survival, whereas other subtypes did not exhibit such benefit. In particular, expression of FMO4 and OSBPL3 was associated with treatment response. Importantly, the improved outcome in the beva+CCNU treatment arm was not explained by an uneven distribution of prognostically favorable subtypes as all molecular glioma subtypes were evenly distributed along the different study arms. The RNA-seq analysis also highlighted genetic alterations, including mutations, gene fusions, and copy number changes, within this well-defined cohort of tumors that may serve as useful predictive or prognostic biomarkers of patient outcome. Further validation of the identified molecular markers may enable the future stratification of recurrent GBM patients into appropriate treatment regimens. Cancer Res; 76(3); 525–34. ©2016 AACR.

Glioblastomas (GBM) are the most common and aggressive type of glial brain tumors (1). The current standard of care for GBM patients includes surgical resection followed by combined chemoirradiation with temozolomide (2). However, GBMs invariably relapse and when this progression occurs, treatment options are limited. Nitrosoureas (in particular lomustine and carmustine), retreatment with (dose intense) temozolomide, re-irradiation, and re-resection are treatments that are often used but have limited activity (3). Progression-free survival (PFS) of recurrent GBM patients is in the range of 2 to 4 months, and post-progression survival of about 6 to 8 months, with conventional chemotherapy (4).

GBMs are histologically characterized by endothelial proliferation and necrosis. At the molecular level, they are characterized by hypoxia-induced and HIF1α-dependent upregulation of VEGF production by tumor cells, which subsequently induces new blood vessel formation (5). Because of the extensive endothelial proliferation that characterizes GBM, soon after the discovery of VEGF and its significance for the angiogenesis of tumor growth it has been hypothesized that GBM would provide a good target for anti-angiogenic treatments.

In spite of some initially promising results in uncontrolled trials however, recent data suggest that the effects of angiogenesis inhibition on overall survival (OS) is limited in both primary and recurrent GBM patients (6–9). For example, bevacizumab (beva), a humanized monoclonal antibody against circulating VEGF, failed to demonstrate significant improvement of survival in newly diagnosed GBM patients in two large randomized phase III clinical trials (10, 11). Similarly, cediranib, a pan VEGFR inhibitor, did not improve outcome in recurrent GBM patients (12). Although the overall effects of angiogenesis inhibition may be disappointing, it is possible that subsets of patients do obtain some clinical benefit (for examples in glioma, see refs. 13–16). For bevacizumab, gene expression analysis provided preliminary evidence of benefit in distinct molecular subtypes of GBMs (17, 18)

Recently, results were reported for the BELOB trial, a trial from the Dutch Neuro-Oncology Group (LWNO), in which patients were randomly assigned to treatment with lomustine (CCNU), bevacizumab or a combination of CCNU with bevacizumab (beva/CCNU; ref. 19). The trial results were intriguing as they indicate a potential survival benefit of beva/CCNU in recurrent GBM patients. In this study, we have performed gene expression profiling and RNA-sequencing (RNA-seq) on tumor samples derived from patients treated within the BELOB trial to identify recurrent GBM patients who benefit from combined beva/CCNU treatment. Our results show that of the different molecular subtypes of GBM (20, 21) “classical” GBMs or those assigned to IGS-18 showed a trend toward benefit from treatment; other subtypes did not show such benefit. When validated in an independent dataset, our data will allow selection of recurrent GBM patients that benefit from beva/CCNU combination treatment.

Patients were eligible for the BELOB trial if they were ≥18 years and had a first recurrence of GBM after temozolomide and radiotherapy treatment. Details of the study have been described previously (19). We also selected 37 paired FF–formalin-fixed, paraffin-embedded (FFPE) samples from the Erasmus University Medical Center glioma tumor archive, with FF and the FFPE samples taken as parallel biopsies from the same tumor. All of the FF–FFPE sample pairs were described previously to assess the performance of HuEx_1.0_St arrays (Affymetrix; refs. 20, 22). Use of patient material was approved by the Institutional Review Board of the respective hospitals and patients provided written informed consent according to national and local regulations for the clinical study and correlative tissue studies.

Total RNA extraction, purification, and quantification from FF and FFPE material were reported previously (22). Purified RNA (250 ng) was used for labeling and hybridization on DASL beadchips (run by Service XS); 500 ng was used for RNA-seq on an Illumina TruSeq and approximately 35 to 40 million 40-base paired end-reads were generated per sample. RNA-seq (n = 96) was run by Expression Analysis (Quintiles). DNA extraction was performed using the Allprep FFPE DNA/RNA FFPE Kit (Qiagen) according to the manufacturers' instructions. RNA expression profiles were then assigned to one of six intrinsic molecular subtypes of glioma (omitting IGS-0), or to one of four glioma subtypes as defined by The Cancer Genome Atlas (TCGA), using the ClusterRepro R package (http://crantastic.org/packages/clusterRepro; ref. 23). Expression data are available at the NCBI Geo datasets, accession number GSE72951.

Gene expression levels (Ref-seq genes) were extracted from the RNA-seq data using featureCounts (24), after alignment on hg19 with Tophat2 (25) of clipped/trimmed reads as provided by the manufacturer. VarScan 2 was used to identify SNVs and indels in the RNA-seq data (26–28), Annovar was used to annotate the SNVs (29). Candidate SNVs were filtered for SNVs that are absent in the 1000 genomes database and for changes that would result in a change in the primary protein sequence. We further focused on those that are either absent from dbSNP138 or are present in the COSMIC database (30, 31). ANKRD36 was removed from this analysis: 497 mutations in 94 samples were identified in this gene, which is likely due to misalignment of homologous and/or non-refseq genes (e.g., ANKRD36B or FLJ54441). Candidate fusion genes were detected using ChimeraScan V0.4.5 on hg19 (32). RT-PCR was used to confirm fusion candidates, and Sanger sequencing was performed to confirm genetic changes.

Differences between the Kaplan–Meier survival curves were calculated using the log-rank (Mantel–Cox) test using GraphPad Prism version 5.00 for Windows (GraphPad Software). For all analysis, OS and PFS were calculated from the point of first recurrence. The significance of prognostic factors was determined with a multivariate analysis using Cox regression. In this analysis, treatment was coded as (i) beva/CCNU, (ii) CCNU, and (iii) beva. Comparisons between frequencies were assessed by the Fisher exact test or the χ2 test (where indicated). SAM analysis was performed using SAMR, an R package. The SAM approach to identify genes associated with treatment response is similar to previously reported using methylation arrays (33). Pathway analysis was performed using Ingenuity IPA (Qiagen) using genes differentially expressed between IGS-18 and all other subtypes at P < 0.01 and >2-fold change in expression level.

DASL and RNA-seq performance on FFPE samples and sample assignment

As RNA isolated from samples that are fixed in FFPE is degraded, we first ran a series of tests to determine the suitability of using DASL arrays and RNA-seq on FFPE-isolated RNA. In this analysis, we compared performance of the two platforms between technical and biologic replicates, and compared the performance between snap-frozen and archival samples that were stored >15 years in paraffin. Results are highlighted in Supplementary Figs. S1 to S3. Our experiments demonstrate that both DASL beadchips and RNA-seq can be used to perform gene expression profiling and assignment to specific molecular subtypes using FFPE tissues, even when the tissues were stored up to 20 years in paraffin.

Molecular subtypes of glioma in BELOB trial tumor samples

All available BELOB trial samples (114/152) were assigned to molecular subtypes as defined by “Gravendeel” (IGS-9, 17, 18, 22, and 23) and “Verhaak” (proneural, neural, classical and mesenchymal; Table 2; refs. 20, 21). The patient characteristics of tumor samples included in the current study did not differ from the entire BELOB patient cohort with respect to age, sex, performance score, MGMT promoter methylation and survival (Table 1). Material was unavailable for the remaining 38 tumor samples. As can be expected, most samples were assigned to the prognostically unfavorable subtypes IGS-18 and IGS-23; GBMs are often assigned to these molecular subtypes (Supplementary Table S3). Few samples were assigned to the prognostically favorable subtypes IGS-9 and IGS-17. All subtypes, including the prognostically favorable samples, were evenly distributed along the different study arms, and the improved outcome in the beva/CCNU arm therefore is not explained by a skewed distribution of these samples. When samples were assigned according to the TCGA classification of GBMs, most samples are assigned to the TCGA “classical” and “mesenchymal” subtypes and few to the “proneural” and “neural” subtypes (Table 2). From the point of first recurrence, the time point used for all analysis in current article, the survival (both OS and PFS) between molecular subtypes was not significantly different (Supplementary Fig. S4).

Similar to previously reported, specific genetic changes segregate into distinct intrinsic subtypes (Supplementary Fig. S5). For example, IDH1 mutations were present in eight BELOB tumor samples (expression data are available for seven), and were identified predominantly in prognostically favorable subtypes IGS-9 (2/4, 50%), IGS-17 (3/9, 33%), and at much lower frequency in IGS-18 (1/72) and IGS-23 (1/28; P < 0.001 IGS-9+17 vs. IGS-18+23). Similarly, IDH1 mutations were predominantly found in “proneural” GBMs (3/8). EGFR amplification is a common event in IGS-18 and in “classical” GBMs (20, 34) and is associated with mRNA upregulation and with an accumulation of genetic changes in this locus (34). In the current dataset, samples assigned to either IGS-18 or “classical” GBMs indeed have a significantly higher expression of EGFR (14.2 vs. 12.8 for IGS-18 vs. any other subtype on DASL, P < 0.001 Students t test; data are identical for “classical” vs. any other subtype). Samples assigned to IGS-18 (or “classical” GBMs) also have a significantly higher frequency of genetic changes in the EGFR locus as determined by RNA-seq (31/47 in IGS-18 vs. 5/29 in any of the other subtypes, P < 0.001, Fisher exact test). Similar for “classical” vs. other GBMs: 34/51 vs. 2/24 (P < 0.001). The assignment of tumors to specific molecular subtypes therefore corresponds to the histologic diagnosis (i.e., predominantly prognostically poor subtypes) and to the type of genetic changes found (i.e., specific genetic changes segregate within defined molecular subtypes).

Response to treatment

The molecular subtypes of gliomas identified in BELOB tumor samples were then stratified according to treatment arm. In this retrospective analysis, samples assigned to IGS-18 showed a trend toward benefit from the combination of beva and CCNU (median OS 7.9, 8.3, and 11.9 months in the CCNU, beva, and beva/CCNU arms, respectively; P = 0.09; Fig. 1). This trend was absent when grouping all samples not assigned to IGS-18 (median OS 10.1, 8.9, and 6.3 months in the CCNU, beva, and beva/CCNU arms, respectively; P = 0.85); too few samples were assigned to IGS-9, IGS-17, IGS-22, or IGS-23 to assess treatment response in individual subtypes. A significant benefit in PFS was observed for tumor samples assigned to IGS-18: survival of samples assigned to IGS-18 was 1.4, 2.9, and 4.2 months in the CCNU, beva, and beva/CCNU arms, respectively (P = 0.0004); for non–IGS-18 samples it was 2.8, 4.1, and 3.7 months (P = 0.23; Fig. 1).

Analysis of subtypes defined by the TCGA showed similar results: “classical” GBMs showed a trend toward improved survival from the combination of beva and CCNU (median OS 8.2, 8.3, and 11.9 months in the CCNU, beva, and beva/CCNU arms, respectively; P = 0.097; median PFS 1.7, 3.0, and 4.2 months, P = 0.001; Supplementary Fig. S6). This trend was absent when grouping all “non-classical” GBMs (median OS 8.7, 7.1, and 5.9 months in the CCNU, beva, and beva/CCNU arms, respectively; P = 0.90; median PFS 2.7, 3.0, and 3.7; P = 0.29); too few samples were “mesenchymal,” “proneural,” or “neural” to assess treatment response in individual subtypes. These data suggest that recurrent GBM patients whose tumors are assigned to IGS-18 or are “classical” GBMs may benefit from the combination of bevacizumab and CCNU, whereas tumors assigned to other subtypes do not show such benefit.

We next aimed to identify genes and molecular pathways that were associated with treatment response on OS. We first performed SAM analysis and identified two genes (FMO4 and OSBPL3) that are associated with OS, specifically in the beva/CCNU treatment arm (using a false discovery rate cutoff of 0.05). SAM analysis failed to identify genes associated with survival in the CCNU and bevacizumab monotherapy arms of this study. As the two genes are associated with survival only in the beva/CCNU treatment arm, these genes may be considered specifically associated with response to treatment (Fig. 2). Indeed, clustering samples based on FMO4 and OSBPL3 expression divided samples in two major subtypes (Fig. 2). OS between the two subtypes was similar for patients treated with either CCNU or bevacizumab monotherapy (median survival of 8.3 vs. 8.5 and 6.1 vs. 10.1 months for CCNU and bevacizumab monotherapy, respectively; P = 0.57). However, survival between the two subtypes significantly differed for patients treated with beva/CCNU (median survival of 6.1 vs. 12.4 months; P < 0.0001). There is no difference in gene expression levels of these two genes between the different molecular subtypes (Supplementary Fig. S8).

Because tumors assigned to IGS-18 show a trend toward benefit from beva/CCNU, we also performed pathway analysis on genes differentially expressed between IGS-18 and all other molecular subtypes. Pathway analysis was performed using our previously reported dataset containing snap-frozen tumor samples (20). A total of 241 differentially expressed genes and 17 differentially activated pathways were identified. Pathways included those that are involved in, among others, cellular assembly, cancer, cellular growth and proliferation. We then used the entire set of genes of each individual pathway (also including genes not differentially expressed between IGS-18 and other molecular subtypes) to screen whether the pathway is associated with response to beva/CCNU treatment. Four of these pathways were associated with response to beva/CCNU treatment, the other 13 pathways showed little association (Supplementary Table S4; Supplementary Fig. S7a and b). Functions of the four pathways associated with response to beva/CCNU combination treatment were often overlapping and include functional descriptors such as “cellular assembly and organization,” “growth and proliferation,” and “tissue/organ morphology.”

Univariate analysis highlighted that in BELOB samples, MGMT promoter methylation (P = 0.01), IDH1 mutation status (P = 0.04), treatment (P = 0.04), and FMO4/OSBPL3 expression (P < 0.0001) were factors significantly correlated with OS. It should be noted however that the study was not powered to perform comparative analysis. Multivariate analysis including these factors showed MGMT promoter methylation, treatment and FMO4/OSBPL3 expression as independent prognostic factors associated with OS (Table 3; Supplementary Table S5). The interaction between FMO4/OSBPL3 expression and treatment was significant, confirming the notion that high FMO4/OSBPL3 expression is associated with treatment response in the bev/CCNU arm.

Factors associated with PFS in a univariate analysis were MGMT promoter methylation (P = 0.004), IDH1 mutation status (P = 0.017), performance status (P = 0.013), intrinsic subtype (P = 0.018), treatment (P = 0.006), and FMO4/OSBPL3 expression (P = 0.0004). Multivariate analysis including these factors showed MGMT promoter methylation, IDH1 mutation status and treatment as independent prognostic factors associated with PFS (Table 4). The interaction between FMO4/OSBPL3 expression and treatment was also significant for PFS in this analysis, confirming the association of FMO4/OSBPL3 expression with bev/CCNU treatment.

In patients of the TCGA database that received bevacizumab at any time during the treatment, a nonsignificant tendency toward increased survival was observed in tumors with high expression of FMO4 or OSBPL3 (strongest in GBMs assigned to IGS-18: median survival of 1.21 vs. 1.50 years for both FMO4 and OSBPL3; Supplementary Fig. S9). High expression of FMO4 or OSBPL3 was not associated with increased survival in the entire TCGA GBM dataset. We should stress that these data should be interpreted with caution: the number of patients receiving bevacizumab (and assigned to IGS-18) is small (n = 18), patients did not receive uniform treatment, and bevacizumab may have been given at a time point different from the BELOB study.

Mutations and structural rearrangements

Mutation analysis using RNA-seq data identified a total of 45 mutations in 13 genes (out of the 71 genes frequently mutated in GBMs (34) in 78 BELOB trial samples (Table 5). Genes include EGFR (n = 29, with often multiple mutations in a single sample), PTEN (n = 10, of which 4 result in a premature stop codon), TP53 (n = 4), and NF1 (n = 2). The frequency of most mutations is slightly lower than previously reported (34) and can be explained because many genes are expressed at low levels, especially when one allele is lost as is the case in many tumor suppressor genes. In addition, homozygous deletions will not be identified by RNA-seq. However, a low frequency for several genes also reflects the high number of “classical” GBMs in our dataset: PIK3CA, PIK3R1, and TP53 mutations are common to “proneural” GBMs, and NF1 mutations are common to “mesenchymal' GBMs (21). The mutations per sample, along with the clinical data and molecular subtype, are shown in Fig. 3. The type and frequency of mutations identified corresponds to those reported previously for GBMs, although RNA-seq is likely to underestimate the mutation frequencies in genes expressed at low level, particularly in tumor suppressor genes (34, 35).

We also used RNA-seq to identify (large) chromosomal losses: regions with chromosomal loss will be represented by RNA-seq as regions absent in heterozygous SNPs. Allele-specific expression will also be reflected as a region absent in heterozygous SNPs (such as occurs in females on the X-chromosome), but this occurs a minority of genes only (36). Our algorithm was tested by comparing control samples, of which, SNP data were also present from snap-frozen samples (Supplementary Fig. S10; see ref. 37). The only frequent LOH observed was LOH of chromosome 10 and was observed at significantly higher frequency in samples assigned to IGS-18 (34/43, 79%) compared with IGS-23 (6/16, 38%) or other molecular subtypes (1/9, 11%, P < 0.001 χ2 test).

ChimeraScan identified a total of 8,879 candidate fusion genes (range, 12–209). Within the technical and biologic replicates, the overlap in identified candidates was low (39.6 ± 9.7%) but improved significantly when additional filters were applied (at least 10 reads covering the breakpoint, of which at least two chimeric reads), to 95% ± 9.2%. When applying these filters to the entire dataset, and removing falsely called fusion transcripts due to alternative splicing events and those that occur in repetitive sequences, 28 candidate fusion genes remained (Supplementary Table S6). These fusion candidates include the known fusion genes FGFR3TACC (n = 2) and EGFRSEPT14 (n = 2; refs. 38–40). WIF1 and VSTM2A were identified as a novel fusion partner of EGFR. All remaining fusion genes were unique events, though some fusion partners (NUP107, VOPP1, GRB10, LANCL2, MLLT3, and VSTM2A) were identified in two samples. Fusion genes are incorporated, alongside with clinical and other molecular data, in Fig. 3.

Fusion genes involving NUP107 surround the MDM2 locus, and it is possible that they occur secondary to a high copy amplification of this locus (both samples express high levels of MDM2). There are four samples with NUP107 fusion genes in the TCGA dataset, all of which have high copy MDM2 amplification (34). By analogy, both PSPH and VOPP1 are located close to the EGFR locus and, in the TCGA dataset, fusions involving PSPH or VOPP1 are either associated with EGFR amplification or are direct fusion partner of EGFR.

Analysis of intronic reads can identify the genomic breakpoint of fusion genes

Sequencing of total RNA results in a large proportion of intronic reads (41). This is likely caused by the use of random primers for cDNA synthesis; in random primed cDNA synthesis both the mature and unspliced mRNAs are converted into cDNA. In general, introns are spliced out only after transcription of the entire intron. Therefore, on a population level, the 5′ end of an intron will have a higher read-depth than the 3′ end (41). This is particularly visible in transcripts with long introns (42), as the half-life of an intron is determined by its length and by the rate of RNA polymerase II transcription, approximately 3.5 to 4 kb/h (43, 44). Of note, we have observed exceptions to this rule, arguing for intra-intronic splicing events (see Supplementary Fig. S11 for some examples). Additional calculations on the rate of transcription and the level of expression are shown in Supplementary Fig. S12.

We hypothesized that the presence of the pre-mRNA transcripts can be used to identify intronic genomic breakpoints of fusion genes: where in intact introns there is a relatively stable coverage of pre-mRNA levels, the levels of pre-mRNA may suddenly change at the exact genomic breakpoint of a fusion gene. Indeed, such sudden decreases/increases of pre-mRNA levels are often present in the introns of the fusion genes identified by ChimeraScan. PCR using primers spanning the putative breakpoint confirmed the presence and location of the genomic break (Supplementary Fig. S13).

In this study, we have performed gene expression profiling and RNA-seq on tumor material derived from patients treated within the BELOB trial in order to identify recurrent GBM patients who benefit from combined CCNU and bevacizumab treatment. Our results indicate that patients with a specific molecular subtype of glioma, IGS-18, or “classical GBMs,” may show more benefit from beva/CCNU treatment. In particular, the expression of FMO4 and OSBPL3 are correlated with benefit from this combination treatment. It should be noted however that our data analysis is post hoc, and confirmation in an independent dataset is therefore required. The in-depth analysis of the transcriptome by RNA-seq also highlighted genetic changes (mutations, indels, gene fusions—including identification of the exact genomic breakpoint—and copy number changes, albeit with limited resolution) within this well-defined cohort of tumors.

Recently, data were reported on two large randomized phase III clinical trials that investigated the role of the addition of bevacizumab to temozolomide chemoradiotherapy in newly diagnosed GBM patients (10, 11). Unfortunately, the results show that the combination treatment did not result in an increased OS of patients. Subsequent translational research however, did identify subgroups of patients that showed benefit from beva + temozolomide treatment (17, 18). For example, Sulman and colleagues provided evidence, in the RTOG 0825 study, that “mesenchymal” GBMs performed particularly poor on the combination treatment (18). Although these data were generated on newly diagnosed GBMs, our data on recurrent GBMs largely corroborate these findings: no trend toward treatment benefit was identified in “non-classical” GBMs (most of which were mesenchymal tumors, though numbers are too small to draw firm conclusions). Recently, Phillips and colleagues, in the AvaGlio study, demonstrated that G-CIMP- “proneural” GBMs benefitted from beva + temozolomide (17). Our dataset, however, contains very few proneural GBMs, which may be related to the observation that these tumors have the worst prognosis of all GBM subtypes and may simply not qualify for second-line treatment (34). Based on data provided in this article, it would be interesting to see whether in those studies tumors assigned to IGS-18 also show benefit from the addition of bevacizumab to chemotherapy. Of note, in a recent retrospective analysis, patients with tumors assigned to IGS-18 or classical GBMs had worse outcome than those assigned to IGS-22/23 when treated with a combination of bevacizumab and irrinotecan (45).

SAM analysis identified two genes that were associated with survival specifically in the beva/CCNU combination arm. Flavin-containing monooxygenase 4 (FMO4) is a protein that catalyzes the NADPH-dependent oxygenation of drugs, pesticides, and xenobiotics (46). FMO4 is part of a protein family (FMO1-5) that, after cytochrome P450, is the second largest protein family involved in drug metabolism. FMO oxidation increases the polarity of (nitrogen-containing) substrates, which aids excretion and detoxification, but may also catalyze drugs into more active forms. Because FMO4 expression is involved in drug metabolism and is associated with survival in the beva+CCNU arm, it is possible that FMO4 expression may render tumors more sensitive to CCNU treatment (when in combination with bevacizumab).

OSBPL3 is one of the 12 members of the oxysterol binding protein (OSBP) related protein (ORP) family that play a role in lipid metabolism, vesicle trafficking, and cell signaling (47). OSBPL3 binds to the phosphoinositides PIP2 and PIP3 and can interact with the small GTPase R-RAS (48, 49). OSBPL3 was shown to play a role in the regulation of the actin cytoskeleton and cellular adhesion (49). Because a set of tumors respond to bevacizumab by increased migration/invasion and vessel co-option (50, 51) and VEGF directly inhibits glioma invasion (52), it is possible that increased OSBPL3 expression increases cellular adhesion and so reduces the migratory capacity of tumor cells. Perhaps a combination of altered migration and drug metabolism renders tumors in the beva+CCNU arm more sensitive to this treatment.

One of the peculiarities of our dataset is that it contained an overrepresentation of tumors assigned to IGS-18 or “classical” GBMs where, in general, GBMs are roughly equally distributed across the various molecular subtypes (20, 21, 53, 54). We are confident however of the molecular classification as described in this article: There was a high correlation in tumor assignment between snap-frozen and FFPE samples and between different methods (DASL and RNA-seq) and genetic changes in the EGFR locus (and 10q LOH) are a hallmark of IGS-18 or “classical” GBMs and most tumors assigned to IGS-18 or “classical” GBMs indeed harbor genetic changes in this gene. This segregation of genetic changes was also observed for IDH1. The higher frequency of tumors assigned to IGS-18 or “classical” GBMs may reflect a sample bias of current study. However, in a study analyzing EGFR amplification in matched primary and recurrent tumors, we also find a high frequency of EGFR amplification at initial diagnosis (∼70%; van den Bent and colleagues; submitted for publication). It is therefore possible that (EGFR-amplified) tumors assigned to IGS-18 are more prone to receiving re-treatment at tumor recurrence. For example, G-CIMP proneural GBMs are selected against as they have poorest prognosis of all GBM subtypes (34).

A limitation of our study is the use of primary tumor tissue for classification, whereas the recurrent tumor was treated. In 2006, a study identifying three molecular subtypes of glioma included analysis of 26 (unselected) matched primary and recurrent tumor pairs (54). This study showed that at progression, gliomas shift from a proneural toward a mesenchymal subtype. However, it also demonstrated that 18/26 (70%) tumors remain in the same subtype (a similar retention of molecular subtype is observed when these tumors to molecular subtypes as defined by Gravendeel and colleagues; Supplementary Table S7). More specifically, within the eight tumors assigned to IGS-18 at initial diagnosis in this study, six remained in IGS-18 at tumor recurrence. A study performed by our group included seven repeat samples and both tumors assigned to IGS-18 at initial diagnosis remained in IGS-18 at tumor recurrence (20). Moreover, we have recently performed analysis of the EGFR locus in 55 uniformly treated primary recurrent tumor pairs (van den Bent and colleagues; submitted for publication). EGFR amplification can be used as surrogate marker for IGS-18 tumors: retention of EGFR amplification status at tumor recurrence therefore is suggestive of retention of the intrinsic glioma subtype. The data from this study indicate that the EGFR status of the tumor remains similar at the point of recurrence in 47 of 55 matched tumor pairs. These studies suggest that in many tumors, the molecular subtype largely remains identical at tumor progression (even though they may shift toward a more mesenchymal phenotype).

A second limitation of this study is that material was available for only a subset (114/152) of BELOB trial samples. Although we did not identify differences in patient characteristics of tumor samples included versus not included (Table 1), the use of only a subset of patients may have introduced a sample bias.

Recurrent fusion genes identified in our dataset (FGFR3TACC3 and EGFRSEPT-14) were identified by others with similar frequencies (34, 38, 39, 55). A recent analysis identified the genomic landscape of fusion genes in 13 tumor types, including GBMs (40). Two identical fusion genes (LANCL2–SEPT14 and ZZEF1–ANKFY1) can be detected when overlaying the identified fusions in that study with those reported here. Based on the amplification status of MDM2 and EGFR and the fact that they are present as double minutes when amplified, we hypothesize fusions involving NUP017, PSPH, or VOPP1 occur secondary to high copy gene amplification. If so, this implies that GBMs harbor few functional fusion genes.

Apart from the recurrent fusions, many of the 5′ and 3′ fusion partners were also identified by others (e.g., GRB10, LANCL2, MLLT3, ASCC3, VOPP1, TERT, VSTM2A, PSPH, NUP107, and LEMD3), some of which are found at significant frequency (GRB10, NUP107, VOPP1, TERT, PSPH; refs. 38–40, 55). One of the potentially functional fusions partner is TERT: TERT fusions can represent a method to increase telomerase activity, as alternative to the frequent TERT promoter mutations in GBMs (56). A total of 7 samples (of which 2 GBMs) were recently identified to contain TERT fusion genes. In all cases identified, TERT acted as the 3′ fusion partner, and all contained the telomerase transcriptase domain (transcripts from exon 2 or 3 onward), which argues against random occurrence of the fusion.

Many of the reads identified by RNA-seq of FFPE isolated material mapped to intronic regions and represent pre-mRNA species. Similar to previously reported, our data show that the 5′ end of an intron has a higher read-depth than the 3′ end, at least in large introns (41). The 5′ to 3′ decline in intronic reads is linear and is explained by a mechanism where introns are first entirely transcribed before being spliced out (41). However, some genes have introns that are apparently spliced out before transcription of the entire intron (i.e., intra-intronic splicing). Apart from a basic insight into the mechanism of RNA maturation, we show that intronic reads can be used to map genomic breakpoints of fusion genes.

In summary, our results show that tumors assigned to IGS-18 or classical GBMs showed a significant benefit in PFS and a trend toward benefit in OS from beva+CCNU treatment; other subtypes did not show such benefit. Expression of FMO4 and OSBPL3, genes involved in drug metabolism and cell signaling, regulation of the actin cytoskeleton and cellular adhesion were specifically associated with treatment response. When validated in an independent dataset, our data will allow selection of recurrent GBM patients that benefit from beva+CCNU treatment.

M.J. van den Bent reports receiving commercial research grants from Abbvie and Roche and is consultant/advisory board member for Roche, Cavion, Novocure, Celldex, Abbvie, Actelion, and BMS. H.M. Oosterkamp reports receiving commercial research grant from Roche and is a consultant/advisory board member for the same. P.J. French is consultant/advisory board member for Roche. No potential conflicts of interest were disclosed by the other authors.

Conception and design: M.J. van den Bent, W. Taal, L.V. Beerepoot, A.H. Honkoop, P.J. French

Development of methodology: M.J. van den Bent, Y. Hoogstrate, P. van der Spek, W. Taal, P.J. French

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M.J. van den Bent, M. de Wit, W. Taal, H.M. Oosterkamp, A. Walenkamp, L.V. Beerepoot, M.C.J. Hanse, J. Buter, A.H. Honkoop, R.M. Vernhout, J.M. Kros

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): L. Erdem-Eraslan, M.J. van den Bent, Y. Hoogstrate, H. Naz-Khan, A. Stubbs, P. van der Spek, R. Böttcher, A. Walenkamp, J.M. Kros, P.J. French

Writing, review, and/or revision of the manuscript: L. Erdem-Eraslan, M.J. van den Bent, Y. Hoogstrate, A. Stubbs, Y. Gao, M. de Wit, A. Walenkamp, L.V. Beerepoot, J. Buter, A.H. Honkoop, B. van der Holt, P.A.E. Sillevis Smitt, J.M. Kros, P.J. French

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M.J. van den Bent, B. van der Holt, J.M. Kros

Study supervision: M.J. van den Bent, P. van der Spek, W. Taal, A. Walenkamp, R.M. Vernhout, P.A.E. Sillevis Smitt, P.J. French

The study was supported by the Stichting Stophersentumoren.nl 2013 and by grant nr DDHK 2010-4678 from the “KWF Kankerbestrijding” (Dutch Cancer Society). RNA sequencing was performed by EA/Quintiles as part of an RNA-seq research grant awarded by EA/Quintiles and Illumina, Inc. The BELOB clinical trial was financially supported by Roche Netherlands.

1.
Louis
DN
,
Ohgaki
H
,
Wiestler
OD
,
Cavenee
WK
. 
WHO classification of tumours of the central nervous system, 4th edition
.
Lyon
:
World Health Organization
; 
2007
.
2.
Stupp
R
,
Mason
WP
,
van den Bent
MJ
,
Weller
M
,
Fisher
B
,
Taphoorn
MJ
, et al
Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma
.
N Engl J Med
2005
;
352
:
987
96
.
3.
Weller
M
,
Cloughesy
T
,
Perry
JR
,
Wick
W
. 
Standards of care for treatment of recurrent glioblastoma–are we there yet
?
Neuro Oncol
2013
;
15
:
4
27
.
4.
Gorlia
T
,
Stupp
R
,
Brandes
AA
,
Rampling
RR
,
Fumoleau
P
,
Dittrich
C
, et al
New prognostic factors and calculators for outcome prediction in patients with recurrent glioblastoma: a pooled analysis of EORTC brain tumour group phase I and II clinical trials
.
Eur J Cancer
2012
;
48
:
1176
84
.
5.
Fischer
I
,
Gagner
JP
,
Law
M
,
Newcomb
EW
,
Zagzag
D
. 
Angiogenesis in gliomas: biology and molecular pathophysiology
.
Brain Pathol
2005
;
15
:
297
310
.
6.
Friedman
HS
,
Prados
MD
,
Wen
PY
,
Mikkelsen
T
,
Schiff
D
,
Abrey
LE
, et al
Bevacizumab alone and in combination with irinotecan in recurrent glioblastoma
.
J Clin Oncol
2009
;
27
:
4733
40
.
7.
Vredenburgh
JJ
,
Desjardins
A
,
Herndon
JE
 2nd
,
Dowell
JM
,
Reardon
DA
,
Quinn
JA
, et al
Phase II trial of bevacizumab and irinotecan in recurrent malignant glioma
.
Clin Cancer Res
2007
;
13
:
1253
9
.
8.
Vredenburgh
JJ
,
Desjardins
A
,
Herndon
JE
 2nd
,
Marcello
J
,
Reardon
DA
,
Quinn
JA
, et al
Bevacizumab plus irinotecan in recurrent glioblastoma multiforme
.
J Clin Oncol
2007
;
25
:
4722
9
.
9.
Khasraw
M
,
Ameratunga
MS
,
Grant
R
,
Wheeler
H
,
Pavlakis
N
. 
Antiangiogenic therapy for high-grade glioma
.
Cochrane Database Syst Rev
2014
;
9
:
CD008218
.
10.
Chinot
OL
,
Wick
W
,
Mason
W
,
Henriksson
R
,
Saran
F
,
Nishikawa
R
, et al
Bevacizumab plus radiotherapy–temozolomide for newly diagnosed glioblastoma
.
N Engl J Med
2014
;
370
:
709
22
.
11.
Gilbert
MR
,
Dignam
JJ
,
Armstrong
TS
,
Wefel
JS
,
Blumenthal
DT
,
Vogelbaum
MA
, et al
A randomized trial of bevacizumab for newly diagnosed glioblastoma
.
N Engl J Med
2014
;
370
:
699
708
.
12.
Batchelor
TT
,
Mulholland
P
,
Neyns
B
,
Nabors
LB
,
Campone
M
,
Wick
A
, et al
Phase III randomized trial comparing the efficacy of cediranib as monotherapy, and in combination with lomustine, versus lomustine alone in patients with recurrent glioblastoma
.
J Clin Oncol
2013
;
31
:
3212
8
.
13.
Hegi
ME
,
Diserens
AC
,
Gorlia
T
,
Hamou
MF
,
de Tribolet
N
,
Weller
M
, et al
MGMT gene silencing and benefit from temozolomide in glioblastoma
.
N Engl J Med
2005
;
352
:
997
1003
.
14.
van den Bent
MJ
,
Brandes
AA
,
Taphoorn
MJ
,
Kros
JM
,
Kouwenhoven
MC
,
Delattre
JY
, et al
Adjuvant procarbazine, lomustine, and vincristine chemotherapy in newly diagnosed anaplastic oligodendroglioma: long-term follow-up of EORTC brain tumor group study 26951
.
J Clin Oncol
2013
;
31
:
344
50
.
15.
Erdem-Eraslan
L
,
Gravendeel
LA
, de
Rooi
J
,
Eilers
PH
,
Idbaih
A
,
Spliet
WG
, et al
Intrinsic molecular subtypes of glioma are prognostic and predict benefit from adjuvant procarbazine, lomustine, and vincristine chemotherapy in combination with other prognostic factors in anaplastic oligodendroglial brain tumors: a report from EORTC study 26951
.
J Clin Oncol
2013
;
31
:
328
36
.
16.
Cairncross
JG
,
Wang
M
,
Jenkins
RB
,
Shaw
EG
,
Giannini
C
,
Brachman
DG
, et al
Benefit from procarbazine, lomustine, and vincristine in oligodendroglial tumors is associated with mutation of IDH
.
J Clin Oncol
2014
;
32
:
783
90
.
17.
Phillips
H
,
Sandmann
T
,
Li
C
,
Cloughesy
TF
,
Chinot
OL
,
Wick
W
, et al
Correlation of molecular subtypes with survival in AVAglio (bevacizumab [Bv] and radiotherapy [RT] and temozolomide [T] for newly diagnosed glioblastoma [GB])
.
J Clin Oncol
2014
;
32
:
5S:suppl abstract 2001
.
18.
Sulman
EP
,
Won
M
,
Blumenthal
DT
,
Vogelbaum
MA
,
Colman
H
,
Jenkins
RB
, et al
Molecular predictors of outcome and response to bevacizumab (BEV) based on analysis of RTOG 0825, a phase III trial comparing chemoradiation (CRT) with and without BEV in patients with newly diagnosed glioblastoma (GBM)
.
J Clin Oncol
2013
;
31
:
suppl; abstr LBA2010
.
19.
Taal
W
,
Oosterkamp
HM
,
Walenkamp
AM
,
Dubbink
HJ
,
Beerepoot
LV
,
Hanse
MC
, et al
Single-agent bevacizumab or lomustine versus a combination of bevacizumab plus lomustine in patients with recurrent glioblastoma (BELOB trial): a randomised controlled phase 2 trial
.
Lancet Oncol
2014
;
15
:
943
53
.
20.
Gravendeel
LA
,
Kouwenhoven
MC
,
Gevaert
O
, de
Rooi
JJ
,
Stubbs
AP
,
Duijm
JE
, et al
Intrinsic gene expression profiles of gliomas are a better predictor of survival than histology
.
Cancer Res
2009
;
69
:
9065
72
.
21.
Verhaak
RG
,
Hoadley
KA
,
Purdom
E
,
Wang
V
,
Qi
Y
,
Wilkerson
MD
, et al
Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1
.
Cancer Cell
2010
;
17
:
98
110
.
22.
Gravendeel
LA
,
de Rooi
JJ
,
Eilers
PH
,
van den Bent
MJ
,
Sillevis Smitt
PA
,
French
PJ
. 
Gene expression profiles of gliomas in formalin-fixed paraffin-embedded material
.
Br J Cancer
2012
;
106
:
538
45
.
23.
Kapp
AV
,
Tibshirani
R
. 
Are clusters found in one dataset present in another dataset
?
Biostatistics
2007
;
8
:
9
31
.
24.
Liao
Y
,
Smyth
GK
,
Shi
W
. 
featureCounts: an efficient general purpose program for assigning sequence reads to genomic features
.
Bioinformatics
2014
;
30
:
923
30
.
25.
Trapnell
C
,
Pachter
L
,
Salzberg
SL
. 
TopHat: discovering splice junctions with RNA-Seq
.
Bioinformatics
2009
;
25
:
1105
11
.
26.
Koboldt
DC
,
Zhang
Q
,
Larson
DE
,
Shen
D
,
McLellan
MD
,
Lin
L
, et al
VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing
.
Genome Res
2012
;
22
:
568
76
.
27.
Koboldt
DC
,
Chen
K
,
Wylie
T
,
Larson
DE
,
McLellan
MD
,
Mardis
ER
, et al
VarScan: variant detection in massively parallel sequencing of individual and pooled samples
.
Bioinformatics
2009
;
25
:
2283
5
.
28.
Carrara
M
,
Beccuti
M
,
Lazzarato
F
,
Cavallo
F
,
Cordero
F
,
Donatelli
S
, et al
State-of-the-art fusion-finder algorithms sensitivity and specificity
.
Biomed Res Int
2013
;
2013
:
340620
.
29.
Wang
K
,
Li
M
,
Hakonarson
H
. 
ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data
.
Nucleic Acids Res
2010
;
38
:
e164
.
30.
Genomes Project
C
,
Abecasis
GR
,
Altshuler
D
,
Auton
A
,
Brooks
LD
,
Durbin
RM
, et al
A map of human genome variation from population-scale sequencing
.
Nature
2010
;
467
:
1061
73
.
31.
Forbes
SA
,
Bindal
N
,
Bamford
S
,
Cole
C
,
Kok
CY
,
Beare
D
, et al
COSMIC: mining complete cancer genomes in the Catalogue of Somatic Mutations in Cancer
.
Nucleic Acids Res
2011
;
39
:
D945
50
.
32.
Iyer
MK
,
Chinnaiyan
AM
,
Maher
CA
. 
ChimeraScan: a tool for identifying chimeric transcription in sequencing data
.
Bioinformatics
2011
;
27
:
2903
4
.
33.
van den Bent
MJ
,
Gravendeel
LA
,
Gorlia
T
,
Kros
JM
,
Lapre
L
,
Wesseling
P
, et al
A hypermethylated phenotype is a better predictor of survival than MGMT methylation in anaplastic oligodendroglial brain tumors: a report from EORTC study 26951
.
Clin Cancer Res
2011
;
17
:
7148
55
.
34.
Brennan
CW
,
Verhaak
RG
,
McKenna
A
,
Campos
B
,
Noushmehr
H
,
Salama
SR
, et al
The somatic genomic landscape of glioblastoma
.
Cell
2013
;
155
:
462
77
.
35.
Parsons
DW
,
Jones
S
,
Zhang
X
,
Lin
JC
,
Leary
RJ
,
Angenendt
P
, et al
An integrated genomic analysis of human glioblastoma multiforme
.
Science
2008
;
321
:
1807
12
.
36.
Zhang
K
,
Li
JB
,
Gao
Y
,
Egli
D
,
Xie
B
,
Deng
J
, et al
Digital RNA allelotyping reveals tissue-specific and allele-specific gene expression in human
.
Nat Methods
2009
;
6
:
613
8
.
37.
Bralten
LB
,
Kloosterhof
NK
,
Gravendeel
LA
,
Sacchetti
A
,
Duijm
EJ
,
Kros
JM
, et al
Integrated genomic profiling identifies candidate genes implicated in glioma-genesis and a novel LEO1-SLC12A1 fusion gene
.
Genes Chromosomes Cancer
2010
;
49
:
509
17
.
38.
Singh
D
,
Chan
JM
,
Zoppoli
P
,
Niola
F
,
Sullivan
R
,
Castano
A
, et al
Transforming fusions of FGFR and TACC genes in human glioblastoma
.
Science
2012
;
337
:
1231
5
.
39.
Frattini
V
,
Trifonov
V
,
Chan
JM
,
Castano
A
,
Lia
M
,
Abate
F
, et al
The integrated landscape of driver genomic alterations in glioblastoma
.
Nat Genet
2013
;
45
:
1141
9
.
40.
Yoshihara
K
,
Wang
Q
,
Torres-Garcia
W
,
Zheng
S
,
Vegesna
R
,
Kim
H
, et al
The landscape and therapeutic relevance of cancer-associated transcript fusions
.
Oncogene
2015
;
34
:
4845
54
.
41.
Ameur
A
,
Zaghlool
A
,
Halvardson
J
,
Wetterbom
A
,
Gyllensten
U
,
Cavelier
L
, et al
Total RNA sequencing reveals nascent transcription and widespread co-transcriptional splicing in the human brain
.
Nat Struct Mol Biol
2011
;
18
:
1435
40
.
42.
O'Connor
V
,
Genin
A
,
Davis
S
,
Karishma
KK
,
Doyere
V
,
De Zeeuw
CI
, et al
Differential amplification of intron-containing transcripts reveals long term potentiation-associated up-regulation of specific Pde10A phosphodiesterase splice variants
.
J Biol Chem
2004
;
279
:
15841
9
.
43.
Singh
J
,
Padgett
RA
. 
Rates of in situ transcription and splicing in large human genes
.
Nat Struct Mol Biol
2009
;
16
:
1128
33
.
44.
Brody
Y
,
Neufeld
N
,
Bieberstein
N
,
Causse
SZ
,
Bohnlein
EM
,
Neugebauer
KM
, et al
The in vivo kinetics of RNA polymerase II elongation during co-transcriptional splicing
.
PLoS Biol
2011
;
9
:
e1000573
.
45.
Laffaire
J
,
Di Stefano
AL
,
Chinot
O
,
Idbaih
A
,
Gallego Perez-Larraya
J
,
Marie
Y
, et al
An ANOCEF genomic and transcriptomic microarray study of the response to irinotecan and bevacizumab in recurrent glioblastomas
.
Biomed Res Int
2014
;
2014
:
282815
.
46.
Krueger
SK
,
Williams
DE
. 
Mammalian flavin-containing monooxygenases: structure/function, genetic polymorphisms and role in drug metabolism
.
Pharmacol Ther
2005
;
106
:
357
87
.
47.
Lehto
M
,
Olkkonen
VM
. 
The OSBP-related proteins: a novel protein family involved in vesicle transport, cellular lipid metabolism, and cell signalling
.
Biochim Biophys Acta
2003
;
1631
:
1
11
.
48.
Goldfinger
LE
,
Ptak
C
,
Jeffery
ED
,
Shabanowitz
J
,
Han
J
,
Haling
JR
, et al
An experimentally derived database of candidate Ras-interacting proteins
.
J Proteome Res
2007
;
6
:
1806
11
.
49.
Lehto
M
,
Mayranpaa
MI
,
Pellinen
T
,
Ihalmo
P
,
Lehtonen
S
,
Kovanen
PT
, et al
The R-Ras interaction partner ORP3 regulates cell adhesion
.
J Cell Sci
2008
;
121
:
695
705
.
50.
de Groot
JF
,
Fuller
G
,
Kumar
AJ
,
Piao
Y
,
Eterovic
K
,
Ji
Y
, et al
Tumor invasion after treatment of glioblastoma with bevacizumab: radiographic and pathologic correlation in humans and mice
.
Neuro Oncol
2010
;
12
:
233
42
.
51.
Rubenstein
JL
,
Kim
J
,
Ozawa
T
,
Zhang
M
,
Westphal
M
,
Deen
DF
, et al
Anti-VEGF antibody treatment of glioblastoma prolongs survival but results in increased vascular cooption
.
Neoplasia
2000
;
2
:
306
14
.
52.
Lu
KV
,
Chang
JP
,
Parachoniak
CA
,
Pandika
MM
,
Aghi
MK
,
Meyronet
D
, et al
VEGF inhibits tumor cell invasion and mesenchymal transition through a MET/VEGFR2 complex
.
Cancer Cell
2012
;
22
:
21
35
.
53.
Li
A
,
Walling
J
,
Ahn
S
,
Kotliarov
Y
,
Su
Q
,
Quezado
M
, et al
Unsupervised analysis of transcriptomic profiles reveals six glioma subtypes
.
Cancer Res
2009
;
69
:
2091
9
.
54.
Phillips
HS
,
Kharbanda
S
,
Chen
R
,
Forrest
WF
,
Soriano
RH
,
Wu
TD
, et al
Molecular subclasses of high-grade glioma predict prognosis, delineate a pattern of disease progression, and resemble stages in neurogenesis
.
Cancer Cell
2006
;
9
:
157
73
.
55.
Shah
N
,
Lankerovich
M
,
Lee
H
,
Yoon
JG
,
Schroeder
B
,
Foltz
G
. 
Exploration of the gene fusion landscape of glioblastoma using transcriptome sequencing and copy number data
.
BMC Genomics
2013
;
14
:
818
.
56.
Killela
PJ
,
Reitman
ZJ
,
Jiao
Y
,
Bettegowda
C
,
Agrawal
N
,
Diaz
LA
 Jr
, et al
TERT promoter mutations occur frequently in gliomas and a subset of tumors derived from cells with low rates of self-renewal
.
Proc Natl Acad Sci U S A
2013
;
110
:
6021
6
.

Supplementary data