The EGFRvIII transcriptome in glioblastoma: A meta-omics analysis.

BACKGROUND: EGFR is among the genes most frequently altered in glioblastoma, with exons 2-7 deletions (EGFRvIII) being among its most common genomic mutations. There are conflicting reports about its prognostic role and it remains unclear whether and how it differs in signaling compared with wildtype EGFR.
METHODS: To better understand the oncogenic role of EGFRvIII, we leveraged 4 large datasets into 1 large glioblastoma transcriptome dataset (n = 741) alongside 81 whole-genome samples from 2 datasets.
RESULTS: The EGFRvIII/EGFR expression ratios differ strongly between tumors and range from 1% to 95%. Interestingly, the slope of relative EGFRvIII expression is near-linear, which argues against a more positive selection pressure than EGFR wildtype. An absence of selection pressure is also suggested by the similar survival between EGFRvIII-positive and -negative glioblastoma patients. EGFRvIII levels are inversely correlated with pan-EGFR (all wildtype and mutant variants) expression, which indicates that EGFRvIII has a higher potency in downstream pathway activation. EGFRvIII-positive glioblastomas have a lower CDK4 or MDM2 amplification incidence than EGFRvIII-negative (P = .007), which may point toward crosstalk between these pathways. EGFRvIII-expressing tumors have an upregulation of "classical" subtype genes compared to those with EGFR-amplification only (P = 3.873e-6). Genomic breakpoints of the EGFRvIII deletions have a preference toward the 3'-end of the large intron-1. These preferred breakpoints preserve a cryptic exon resulting in a novel EGFRvIII variant and preserve an intronic enhancer.
CONCLUSIONS: These data provide deeper insights into the complex EGFRvIII biology and provide new insights for targeting EGFRvIII mutated tumors.

CDK4 and MDM2 amplifications appear less frequently in EGFRvIII+ compared with EGFRvIII− but EGFR-amplified GBM.

Transcriptomes of EGFR-amplified GBM differ marginally between EGFRvIII+ and EGFRvIII−.

EGFRvIII breakpoints preferentially retain an intronic enhancer.

Glioblastoma is the most prevalent and aggressive form of malignant primary brain tumors, often characterized by EGFR mutations for which no effective treatments is available. We aimed to understand the role of its most common mutation, EGFRvIII (in-frame deletion of exons 2-7). By exploiting 6 combined datasets, we show the interplay between pan-EGFR and EGFRvIII levels, find no positive selection towards EGFRvIII expression and demonstrate that EGFRwt and EGFRvIII largely activate similar pathways. However, significant and unique EGFRvIII mutation-specific associations were found with Cell Cycle (eg, CDK4) and RTK/RAS/PI3K genes (eg, MDM2) which provide new insights for tumor targeting. A preference in breakpoint location in intron-1 not only results in a distinct variant of EGFRvIII but also preserves an enhancer region, and so provides new insights into EGFR(vIII) gene regulation.

Glioblastoma is the most prevalent and aggressive form of malignant primary brain tumor in adults with a short median survival time of 14.6 months.[1] Extensive research on the genetic makeup of glioblastoma has revealed recurrent genetic changes typically involving the RTK/RAS/PI3K, p53, and RB signaling pathways.[2-4] Although the diverse genetic features of glioblastoma have become increasingly better understood, no effective treatment options are currently available that specifically target the most common mutations. One of the most frequently altered genes in glioblastoma encodes the epidermal growth factor receptor (EGFR). EGFR is amplified in ~50% of all glioblastomas,[4-7] typically within small circular extrachromosomal DNA (ecDNA) copies.[8] The most common mutation on top of this amplification is an in-frame deletion of exons 2-7 (EGFRvIII), found in ~50% of the EGFR-amplified glioblastoma patients.[9] EGFRvIII is a constitutively, but low-level, active form of EGFR that is independent of ligand for its activation,[9] likely due to the partially deleted extracellular ligand-binding receptor domain. The EGFRvIII variant results from a genomic deletion, not from alternative or aberrant splicing. Unfortunately, treatments aimed at targeting EGFRvIII have thus far not provided clinical benefit to patients.[10,11]

It is assumed that EGFRvIII typically is a late event that arises after chromosome 7 amplification and after EGFR high-copy amplifications and is therefore considered subclonal.[12] However, even as subclonal mutation, it is highly prevalent in glioblastoma and contributes to and alters the biology of the tumor. EGFRvIII has been shown to reduce apoptosis and increase proliferation and invasiveness,[9] key features of tumor progression. Protein levels of EGFRvIII vary widely across and spatially within glioblastoma tumors.[13-15] Moreover, recent observations show changes in EGFRvIII levels during tumor evolution after initial resection,[6,16,17] including cases with complete loss of EGFRvIII over time. That such a common presumed driver mutation gets lost, or levels get reduced during tumor evolution is paradoxical and will complicate targeting it for clinical benefit.

In this study, we aim to unravel EGFRvIII-specific mechanisms related to glioblastoma tumorigenesis. We examined EGFRvIII expression, its genomic breakpoints, and co-occurrence with other genetic changes using a large combined dataset.

Methods

Sequencing Data

Sequencing of the Intellance-2[18] (paired-end; 2 × 151 bp total RNA + paired-end; 2 × 76 bp TruSight Tumor 170 panel) and BELOB (single-end; 50 bp)[19] data was described elsewhere. For G-SAM, RNA extraction was performed using the AllPrep DNA/RNA FFPE kit or the RNeasy FFPE kit (Qiagen, Venlo, the Netherlands). G-SAM samples were sequenced (150 bp paired-end reads) on the Illumina NovaSeq at the GIGA-Genomics Core Facility University of Liège. Each of these datasets was non-poly(A)+-enriched and thus also include non-polyadenylated transcripts.[20] Raw sequencing data are available (BELOB: EGAS00001004570, Intellance-2: EGAS00001005437; G-SAM: EGAS00001005436). TCGA-GBM (poly(A)+ RNA and DNA mutations) and CPCT-02 and PCAWG DNA data were obtained from their public repositories.

Human Specimens

Tissue and metadata from the G-SAM and Intellance-2 studies were accrued through the pan-European Organisation for Research and Treatment of Cancer network.[6,18] Informed written consent was obtained from all patients. The study design was approved by the institutional review board of Erasmus MC (Rotterdam, the Netherlands), and conducted according to institutional and national regulations.

Sequencing Data Processing

For each RNA-seq sample, FASTQ files were cleaned using fastp (https://github.com/OpenGene/fastp), aligned to hg19 using STAR[21] and then de-duplicated with sambamba. For the Dr. Disco[20] pipeline, samples were first FASTQ-level de-duplicated level using HTStream deduper (https://github.com/ibest/HTStream). EGFRvIII and EGFRwt expression was estimated directly from BAM files using junction-reads (https://github.com/yhoogstrate/egfr-v3-determinerv0.7.4:--spliced-reads-only). Reads considered EGFRvIII spanned the splice junction of exons 1-8, and reads considered EGFRwt exons 1-2. Samples with <10 such reads were excluded, except for TCGA-GBM, where EGFRwt read counts for EGFRvIII-negative samples were missing. Junction read counts of replicated samples were merged by summing the spliced read counts. The EGFRvIII percentage was defined as the average percentage when matching data from both sequencing assays were present (Intellance-2). Gene level read counts were obtained using featureCounts and Gencode v31. EGFRvIII counts from TCGA-GBM were taken from elsewhere.[4] Junction-counts involving non-canonical exons A, B, and C were determined using egfr-v3-determiner with modified exon annotations. Genomic events were taken from processed WES data or public resources (Supplementary Methods).

Expression Analysis

Samples with an EGFRwt + EGFRvIII read count ≥10 were eligible for EGFRvIII status and percentage determination and for differential gene expression (DE) analysis. For DE analysis, only genes with on average ≥3 reads per sample were included. Only genes marked as “protein_coding” were included. DE analysis was performed using DESeq2 (Wald test),[22] in which EGFRvIII was excluded in estimateSizeFactors to avoid redundant counts. The FDR-adjusted P-value reflects the q-value. For the tests with 4 datasets combined, the intersected protein-coding genes with on average ≥3 reads per sample, per dataset, were included. Normalized expression levels were estimated using DESeq2 followed by the VST transformation (blind=TRUE) to ensure homoscedasticity.[22] A batch correction was performed for DE and for correlation analysis to correct per-dataset differences (DESeq2 for count data; limma::removeBatchEffect[23] for VST-transformed data). Volcano plots were generated with the EnhancedVolcano package (https://github.com/kevinblighe/EnhancedVolcano). Kaplan-Meier analysis was performed using R’s survival package. Survival analysis on EGFRvIII expression was performed with a Cox Proportional Hazard survival using R’s survival package on the normalized VST-transformed expression values. Because the Depatux-M antibody binds EGFRvIII with high affinity[24] and the Intellance-2 trial reported a benefit from Depatux-M in EGFRvIII-positive samples,[18] Depatux-M arms were excluded from survival analysis.

Breakpoint Analyses

Non-poly(A)+-enriched RNA-seq samples include relatively large proportions of intronic reads derived from actively transcribed pre-mRNA. This allows detection of genomic breakpoints when corresponding introns are sufficiently covered.[19,20] Settings for Chimeric alignment are given in Supplementary Methods.

The 100-vertebrates-phastCons track was obtained from UCSC and smoothened by a running mean of 200 bp fixed windows. H3K27ac Chip-Seq data were obtained from GSM3382305,[25] GSM3670052, GSM3670055, and GSM3670058.[26] Actual genomic enhancer locations were not provided in the original manuscript.[25] Their raw CRISPRi-assay data (GSM4141363 + GSM4141364) were used to reproduce their findings according to their described methodology (Supplementary Table 3).

Exon-B Variant Experiments

To confirm the EGFRvIII exon-B variant, 10 samples positive for the variant (RNA-seq) with remaining isolated RNA leftover from sequencing were chosen (Supplementary Table 2). cDNA was synthesized in a buffer of 1 µl random primers, 1 µl dNTP mix, 1000 ng RNA, and 13 µl dH2O. The mixture was heated to 65°C for 5 minutes and incubated on ice for 1 minute. After brief centrifugation, the contents were collected, and the following was added: 4 µl 5× First-Strand buffer, 1 µl 0.1 M DTT, 0.5 µl RNaseOUT, 1 µl Superscript III. The mixture was incubated at 25°C for 5 minutes, at 50°C for 45 minutes, and inactivated at 70°C for 15 minutes. Partial sequences spanning the exon-B splice junction were PCR-amplified using 6 primer combinations (2× exon-B, 1× exon-8, 2× exon-9). For each reaction, the buffer consisted of: 7.9 µl nuclease-free water, 3 µl 5× GoTaq buffer, 0.8 µl 10 mM dNTPs, 1 µl 10 µM forward primer, 1 µl 10 µM reverse primer, 1 µl cDNA, and 0.3 µl GoTaq polymerase. Denaturation of cDNA was performed at 98°C for 30 seconds, followed by 40 cycles of 30 seconds at 98°C, 30 seconds at 60°C, and 30 seconds at 72°C. The final extension was performed at 72°C for 5 minutes and brought back to 12°C. Of the 10 samples, 6 showed bands of the expected size on agarose gel. Of these 6, 4 were sent out for Sanger Sequencing to Macrogen Europe B.V., Amsterdam (Supplementary Table 2). Three of the four samples showed good per-base quality. Forward and reverse reads were assembled into consensus contigs using UGENE.

Constructs were generated to evaluate the function of EGFR variants initiating from exon-B. Because exon-B lacks a translation initiation site, we generated these constructs using the first in-frame ATG in exon-2 or exon-8 (in the case of EGFRvIII). A total of 16 different constructs were made: those that initiated translation in exons 2 or 8 with (i) either an in-frame eGFP (located C-terminal to the transmembrane region[27]) or eGFP co-expressed via an IRES sequence; (ii) with and without the L858R activating mutation (to compare the activation state of the novel isoforms with a constitutively active isoform); and (iii) without/with a canonical Kozak sequence (to ensure optimal translation of the latter). Constructs were generated by in-fusion cloning into a piggyback vector. Constructs were stably transfected in HeLa cells, imaged using an Opera Phenix (PerkinElmer, Hamburg, Germany) high content imager and analyzed using Harmony software (PerkinElmer, Hamburg, Germany).

Results

We have collected glioblastoma data from the following cohorts: BELOB,[19] Intellance-2,[18] G-SAM,[6] TCGA-GBM,[4] CPCT-02,[28] and PCAWG.[29] The compiled results are available as a study dataset: https://zenodo.org/record/4792445.

Molecular Differences of EGFRvIII-Expressing Tumors

To determine EGFRwt (spliced across exons 1-2) and EGFRvIII (spliced across exons 1-8) expression, we first developed egfr-v3-determiner (publicly available, see Methods). Out of the 839 available RNA-seq samples, we included samples with a combined EGFR junction read count (spliced across exons 1-2 and 1-8) of ≥10 into our combined study RNA dataset: n = 741 from 622 patients; BELOB (n = 69/92), Intellance-2 (n = 224/239), G-SAM (n = 285/345) complemented with all primary TCGA-GBM samples (n = 163). In this combined dataset, 464/741 (62.6%) samples had EGFR gene amplification or upregulation if copy-number data were absent. Using the transcript-specific junction-counts, we calculated the ratio EGFRvIII (. Of the EGFR-amplified samples, 225/464 (48.5%) were considered EGFRvIII-expressing (, consistent with observations in the literature.[9,30] These ratios revealed a high dynamic from 1% to 95%, consistent in all datasets (Figure 1). Lower EGFRvIII expression ratios were slightly over-represented (1%-10%; P = 3.2e−9; Wilcoxon test on the first derivative of the ordered percentages). The total EGFR expression levels are on average lower for samples with higher EGFRvIII percentages, implying that EGFRvIII is more potent in EGFR signaling (Figure 2, Supplementary Figure S1D).

Fig. 1

Range of EGFRvIII percentages relative to total EGFR. Results are split per dataset (top) and combined (bottom). Gray vertical lines (Intellance-2) indicate levels determined by both full and panel-based RNA-seq where the actual percentages reflect their mean. Mutation statuses are indicated underneath. N/A-values are indicated in black.

Fig. 2

EGFRwt/EGFRvIII correlations. (A) EGFRwt and EGFRvIII correlation and (B) total EGFR and percentage of correlation, per dataset. (A) The correlations between EGFRwt and EGFRvIII are negative; (B) y-axis represents a surrogate for the total EGFR level (VST-transformed sum of EGFRwt + EGFRvIII junction-reads because the full gene EGFR read count is negatively affected by exons missing in EGFRvIII). Correlations are negative, indicating that tumors with higher proportions of EGFRvIII have lower levels of both variants combined.

Several reports have indicated that EGFRvIII and EGFRwt activate different signal transduction pathways.[9,31] To assess if such differences are reflected in their transcriptomes, we performed DE analysis comparing the transcriptomes of EGFRvIII-positive (using 2 cutoffs: ≥1.0% or ≥10.0%) with EGFRvIII-negative (<1.0%) but EGFR-amplified tumors. Tests were performed for all 4 datasets separately, to correlate the logarithmic fold changes (LFC) of the genes between the datasets. Markedly higher LFC correlations were found across the datasets using ≥10% EGFRvIII as cutoff (Supplementary Figure S2), which suggests lower percentages (1-10) harbor a limited EGFRvIII response signal.

We therefore proceeded with the combined dataset using only ≥10% EGFRvIII as cutoff (n = 368) and found 213 genes significantly (q-value < 0.01, |LFC| > 0.5) differentially expressed (Figure 3A). They showed enrichment in genes related to microtubule-, cilium-, and axoneme-related pathways (Supplementary Figure S3).

Fig. 3

(A) DE analysis between EGFR-amplified samples with (≥10%) and without EGFRvIII (<1%), with batch correction for the 4 datasets (Intellance-2, G-SAM, BELOB, and TCGA-GBM). 213/15.617 protein-coding genes were differentially expressed, including DLX1, DLX2, TSPAN31, TMPRSS7, PPBP, and DPT. Classical subtype genes are marked black. Overall LFCs were more often negative while the majority of the classical subtype genes had a positive LFC. (B) First two components of a supervised principal component analysis (213 DE genes). (C) Z-scores of Pearson correlation tests between genes and the relative EGFRvIII (x-axis) and EGFRwt (y-axis) levels, in samples with ≥10% EGFRvIII. Values near 0 represent no correlation, negative values represent a negative correlation, and positive values represent a positive correlation. Classical subtype genes are marked black. Genes with a significant difference (t test; q-value < 0.01) are marked purple. Genes showing a trend (q-value < 0.1) are marked blue.

The 187 significantly downregulated genes in ≥10% EGFRvIII included CDK4 and MDM2, genes that are frequently hyper- and co-amplified in glioblastoma. Their observed differences were not a result of consistent down-regulation of CDK4 or MDM2 across all ≥10% EGFRvIII-positive patients but were caused by a lower proportion of tumors with extremely high CDK4 or MDM2 expression levels (Supplementary Figure S4A and B). Integration with copy-number data confirmed the negative association: CDK4 or MDM2 DNA amplifications appeared in significantly fewer tumors expressing EGFRvIII (P = .007, Fisher exact test, Supplementary Figure S4C). TP53 mutations were indeed[32] less frequently present in EGFR-amplified tumors, both EGFRvIII-positive and -negative (Figure 1, Supplementary Figure S1). Similarly, TACC3-FGFR3 fusions were indeed[33] exclusively present in EGFR non-amplified tumors. The overall transcriptome differences did not show a strong separation between EGFRvIII-positive and -negative tumors, indicating the overall differences are modest (Figure 3B).

Glioblastomas are classified into 3 transcriptional subtypes: mesenchymal, proneural, and classical. Classification is based on genes that are exclusively upregulated within their subtype.[34,35] The classical subtype is characterized by EGFR amplifications.[35] We observed that almost all classical subtype genes tend to be upregulated in EGFRvIII-positive tumors (P = 3.873e−6; Fisher exact test on positive/negative LFC, Figure 3A) compared with EGFRvIII-negative tumors, all harboring EGFR amplifications. The classical subtype, therefore, is at least partly defined by EGFRvIII-specific signaling. While certain neuronal precursor and stem cell marker, sonic hedgehog pathway, and notch pathway member genes are highly expressed in the classical subtype,[36] these individual pathways did not differ across EGFRvIII-positive/-negative tumors (Supplementary Figure S5A–C). According to a pathway-based glioblastoma classification,[37] two subtypes, proliferative/progenitor (PPR) and mitochondrial (MTC), are associated with RTK pathway amplifications such as EGFR and PDGFRA. Of these, PPR is associates positively with EGFRvIII (Supplementary Figure S5D and E).

In addition to the DE analysis using a defined EGFRvIII expression cutoff, we interrogated the linear correlation between the expression of all genes to the EGFRvIII expression. This analysis was performed within the same ≥10% EGFRvIII-positive samples. CDK4 and MDM2 expression levels did not linearly correlate with EGFRvIII expression. That there is a significant difference in CDK4 and MDM2 expression levels across EGFRvIII-positive and -negative tumors while their expression levels do not correlate with EGFRvIII, is in concordance with the difference in hyper-amplification incidence. To identify genes that correlate differently between EGFRvIII and EGFR, we performed the same test against EGFRwt (Figure 3C). We then calculated per gene to what extent the correlation with EGFRwt and EGFRvIII differs, and tested which differences were beyond what may be expected by chance (Supplementary Methods). This revealed 6 additional genes that significantly differ in their correlation to EGFRwt in contrast to EGFRvIII (NSG1, GALNT15, RFWD3, NCAPD3, ARHGEF26, and PHF19; q < 0.01). RFWD3 was positively correlated with EGFRvIII (coef = 0.33) while negatively correlated with EGFRwt (coef = −0.20). Similar to using a defined EGFRvIII expression cutoff, we found that the classical subtype genes correlate positively stronger with EGFRvIII compared with EGFRwt (P = 1.1e−9; 2-sided t test on Z-score difference).

The difference in correlation between EGFRvIII and EGFRwt and the difference in gene expression by EGFRvIII presence, showed correlation (Spearman coef = 0.4, Supplementary Figure S6). For classical subtype genes, this correlation was stronger (Spearman coef = 0.7), indicating consistency in the outcome of the tests. In particular, genes that showed strong concordant results were PHF19, NSG1, and Sprouty/Spred family members SPRED2, SPRY4, and SPRY2 (Supplementary Figure S6B). Furthermore, PTPRZ1, occasionally found in glioma as donor partner in fusions such as PTPRZ1-ETV1 and PTPRZ1-MET,[38] positively associates with EGFRvIII.

EGFRvIII Prognostic Value

There have been conflicting data on the association of EGFRvIII with prognosis.[9] We interrogated the patient survival between EGFRvIII-positive and -negative patients in the BELOB, G-SAM, and TCGA-GBM and Intellance-2 (control arm) datasets. Within patients with EGFR-amplified tumors, there was no significant difference in overall survival between patients with EGFRvIII-positive and -negative tumors (n = 327) in each dataset or combined (Figure 4, Supplementary Figure S7). There was no significant association between relative EGFRvIII expression levels and patient survival (HR: BELOB = 1.1, G-SAM = 0.96, Intellance-2 = 1.2). In summary, we found no evidence for an association of EGFRvIII with survival in patients with EGFR-amplified tumors.

Fig. 4

Kaplan-Meier survival plots of patients with EGFR amplification with/without EGFRvIII. Patients included were from the BELOB trial, Intellance-2 TMZ/control arm, primary G-SAM tumors, and primary TCGA-GBM tumors. Difference in patient survival between EGFR-amplified glioblastoma patients with/without EGFRvIII (≥1% and ≥10%) was not significant and neither in each dataset separately (Supplementary Figure S7).

EGFRvIII Breakpoints Preferentially Retain Intronic Enhancer

With a transcription rate of 1-6 kb/min,[39] transcription of the ~120 kb EGFR intron-1 can take up to 2 hours. The closer the breakpoint of the causal EGFRvIII deletion is to exon-1, the shorter its intron. Given the large size of intron-1, breakpoints at the beginning of the intron (early breakpoints) may provide an energetic and temporal benefit over breakpoints at the end of the intron (late breakpoints). We screened ≥1% EGFRvIII-positive samples for their genomic EGFRvIII breakpoints based on the presence of pre-mRNA.[19,20] We found 44 breakpoints within our datasets (Supplementary Table 1; Figure 5). One sample harbored 2 unique EGFRvIII breakpoints. We complemented these breakpoints with those identified from CPCT-02 (8/41 patients)[28] and PCAWG (11/40 patients)[29] whole-genome sequencing datasets. In several samples, we observed multiple, unique EGFRvIII breakpoints (Supplementary Figure S8A) that could not have evolved from a tumor-specific ancestor EGFRvIII variant. In these cases, EGFRvIII thus has independently reoccurred within the same tumor.

Fig. 5

Overview of genomic EGFR locus (exons 1-11) and EGFRvIII breakpoints. From bottom to top: chr7, transcript annotations, late and early breakpoint regions, conservation (purple), H3K27ac intensity in GSC23 cells[25] (green), and the actual breakpoints (blue and mustard). Genomic EGFRvIII breakpoints are indicated with mustard (RNA detected) and blue (DNA detected) bars on top.

Interestingly, the genomic breakpoints found in intron-1 show a difference in breakpoint density, where the region close to exon-1 contains 3.63 times fewer breakpoints per base than the region close to exon-2 (P = 4.9e−13; Fisher exact test; decision-boundary: chr7:55.182.397). Genomic breakpoints between exons 7-8 were more uniformly distributed (Supplementary Figure S8B). The breakpoint preference in intron-1 may suggest preserving functional regions that confer a selective advantage to the tumor. Upon closer inspection, EGFR intron-1 contains 3 non-canonical exons[40] of which their expression is only rarely observed. We refer to these as exons A, B, and C. The EGFRvIII breakpoint preference region is located 3′ of exon-B (Figure 5) and thus preserves this exon at the genomic level. All datasets examined revealed junction-reads that initiated in exon-B and were spliced to exon-2 (EGFRwt) or exon-8 (EGFRvIII) (Figure 6). However, the fraction of transcripts containing exon-B was low compared to those initiating in exon-1 (≤1.05%; Supplementary Figure S8C), indicating exon-B expression is driven by a weak promoter. Transcripts spliced from exons A or C to exon-2 were extremely rare.

Fig. 6

Exon-B expression. Spliced read counts for exon-B (exon-B→exon-2: bars up and exon-B→exon-8: bars down) in tumors with RNA detected genomic EGFRvIII breakpoint. Tumors with a “late” intron-1 breakpoint (≥chr7:55.182.397) are marked with a square and “early” with a cross. Regular (A) and high (B) depth datasets were split. EGFRvIII exon-B variant reads (exon-B→exon-8) are only present in tumors with a late EGFRvIII breakpoint, which retains exon-B.

In samples with breakpoints retaining exon-B, a novel exon-B-exon-8 EGFR(vIII/B) variant is created (Figure 6). This variant was confirmed with RT-PCR in 6 out of 10 tested tumor samples (Supplementary Table 2). We verified the presence of the exon-B→exon-8→exon-9 sequence in 3 samples (GenBank: MZ484953, MZ484954, and MZ484955). EGFR transcripts that initiate in exon-B lack part of the extracellular domain on protein level as the translation initiation sites are located in exon-2 or exon-8. To test the potential functional role of exon-B variants, we created constructs of EGFR starting in exon-B and spliced to either exon-2 or exon-8.

Even after optimizing the Kozak sequence surrounding the translation initiation site, we failed to see the expression of “exon-B” variants in any of the 16 constructs generated. This absent expression was not due to a potential lethality of exon-B constructs as (1) RT-PCR did show expression of the EGFR transgene and (2) biscistronic constructs (in which eGFP was independently translated from EGFR constructs as they were separated by an IRES sequence) did express eGFP. These data argue for an inferior protein translation of exon-B transcripts. Given the inadequate translation into protein combined with the low level of transcripts incorporating exon-B, we deemed it unlikely these constructs significantly impact the tumor biology.

We explored the possibility that late breakpoints retain regulatory sequences further. H3K27ac ChIP-seq data from recent studies on EGFR enhancers[25,26] were plotted onto the EGFR locus (Figure 5, Supplementary Figure S9). An enhancer, previously referred to as “E3,” [25] located just 5′ to the EGFRvIII breakpoint preference region and is thus more often preserved. This region is conserved across 100 vertebrates. Previous experiments using CRSIPRi demonstrated its functional relevance in cell fitness. Unfortunately, too few samples with detected breakpoints and combined RNA-seq and DNA-seq data were available to determine whether late breakpoints have a higher fractional EGFRvIII expression (Supplementary Figure S10).

Discussion

EGFR is commonly amplified, mutated, and activated in glioblastoma, resulting in increased cell invasion and proliferation.[41]EGFRvIII is a specific tumor marker often present in glioblastoma, that has been intensively investigated.[9,18] Here, we report on this genomic mutation using a large glioblastoma EGFRvIII omics dataset. To maximize statistical power, analysis was performed across a combined cohort of 4 RNA datasets and 2 independent whole-genome sequencing datasets. Previous data on the prognostic value of EGFRvIII were conflicting, with some suggesting EGFRvIII is a negative[42,43] or a positive[44] prognostic marker, where other studies also suggested it did not affect survival.[45,46] Here, we demonstrate that within patients with EGFR-amplified glioblastoma, we observed no difference in survival between EGFRvIII-positive and -negative tumors. Because EGFRvIII is known to be spatially heterogeneously distributed,[13,14]EGFRvIII-positive tumors can, therefore, through sampling, be marked EGFRvIII-negative by omics analysis. Tumor sampling is therefore a limitation potentially influencing this survival analysis.

The expression levels of EGFRvIII and EGFRwt were anti-correlated and the total EGFR levels were generally lower when higher levels of EGFRvIII were present. This is in agreement with the hypothesis that EGFRvIII lowers the tumors’ dependency on high EGFR amplification levels.[18]

Transcriptomes of EGFRvIII-positive and -negative tumors showed only minor differences (Figure 3B). A possibly related factor of this limited difference may be the ability of EGFRvIII to alter expression in EGFRvIII-negative tumor cells.[15] Within EGFR-amplified tumors, those with ≥10% EGFRvIII were found to have significantly lower expression of CDK4 and MDM2 due to a lower incidence of respective amplifications. This inverse correlation may point toward crosstalk or redundancy between these pathways. Of the genes correlated positively to EGFRvIII expression, RFWD3 can form a complex with MDM2, known for regulating p53.[47] Furthermore, Sprouty/Spred family genes were consistently associated with EGFRvIII presence and subsequent expression and are known for their inhibiting role in Ras/Raf/ERK[48] and involvement in EGF/EGFR signaling. The presented results are not supporting the standpoint that EGFRvIII is causing large distinct changes in downstream gene expression compared with EGFRwt amplifications.

Overall, our molecular analysis demonstrated that glioblastomas expressing EGFRvIII show a distinct but limited difference in their transcriptome compared with EGFRwt. The clearest observed signal is an increased correlation with classical glioblastoma subtype genes, which may indicate that the constitutively active EGFRvIII is, in the context of EGFR-amplified glioblastoma, stronger in downstream EGFR signaling than (amplified) EGFRwt. This is in line with the lower total EGFR levels for tumors having higher EGFRvIII levels.

We found a broad range of EGFRvIII/total EGFR expression levels (1%-95%). Such range is puzzling because, if EGFRvIII is only a variant that is stronger in activating downstream EGFR signaling, it is possible that EGFRvIII would simply outcompete the EGFRwt ecDNA copies. This would likely take place relatively quickly since ecDNA amplifications are notorious for increasing tumor heterogeneity.[8] However, the presence of extrachromosomal EGFRwt copies lasts in virtually all analyzed EGFRvIII-positive tumors. An explanation could be that EGFRvIII depends on the presence of EGFRwt,[49] for instance, to form dimers to complete EGFRvIII phosphorylation[50] or in an inter-cellular context, for instance by EGFRvIII-dependent secretion of cytokines.[15] Such dependencies would likely come with a preferred EGFRvIII/EGFRwt ratio. Alternatively, the linear slope is indicative for an absence of selection pressure to retain EGFRvIII over EGFRwt. This absence can explain the highly heterogeneous spatial and temporal expression pattern of the mutant. It may also explain the near-identical survival between EGFRvIII-positive and -negative glioblastoma patients. However, if there is no selection pressure to retain EGFRvIII, it remains puzzling why this particular mutant is found at such a high frequency. EGFR signaling in glioblastomas is highly complex as the tumor can adopt various methods to enhance its pathway activation. Multiple mutations can co-exist in the same tumor, sometimes subclonal and with reported longitudinal differences, with a unique, different ligand dependency.

An earlier study proposed defining samples with a read count of at least 1% or 10% EGFRvIII compared with total EGFR as EGFRvIII-positive.[4] We recommend similarly rather than using the presence of any EGFRvIII read, as mapping artifacts and index hopping/switching derived reads are common in multiplexed RNA-seq and because higher EGFRvIII percentages showed a stronger response signal.

Determination of the subclonal breakpoints in pre-mRNA data was more complicated than in datasets where breakpoints were clonal.[20] Breakpoints were found predominantly in samples with high fractions of EGFRvIII. The median EGFRvIII percentage in samples with detected breakpoints was 55%, whereas 29% in samples without.

Intriguingly, we find a minority of EGFR transcripts starting with a cryptic exon preferentially preserved in EGFRvIII-expressing tumors. The first translation initiation site is located in exon-8, but the total exon-B read count is very low and, combined with a weak Kozak sequence, we did not consider this variant to be the main reason for a breakpoint preference.

Recently, the promotor and functional enhancers specifically retained in extrachromosomal EGFR fragments in glioblastoma and neuroblastoma cells have been interrogated.[25] These enhancers, including “E3,” were discovered using 4C-seq, H3K27ac ChIP-seq, and a CRISPRi knock-down proliferation dropout assay. The E3 enhancer also showed H3K27ac in an independent dataset.[26] The preferential retention of intragenic enhancer E3 in EGFRvIII is in line with these observations. As the E3 enhancer is also conserved across vertebrates, it likely results in higher EGFR transcription rates. Unfortunately, both absolute and relative EGFRvIII levels differ essentially between samples, which combined with a high level of EGFRwt heterogeneity makes it difficult to confirm this hypothesis.

In summary, using the largest combined EGFRvIII omics dataset to date, we find that the expression profiles of EGFRvIII-positive tumors differ only marginally from EGFRvIII-negative tumors. The results suggest that EGFRvIII mainly performs a similar role as EGFRwt but with a stronger affinity to activate EGFR downstream pathways, possibly linked to persistent activity independent of ligand(s). Furthermore, genomic breakpoints in intron-1 retain an enhancer that likely increases the expression of EGFRvIII transcripts. In this retrospective setting, no prognostic difference was found between EGFRvIII-positive patients compared with those harboring EGFRwt amplifications. However, associations between EGFRvIII and genes such as CDK4, MDM2, and PTPRZ1 suggest that the relation between EGFR and EGFRvIII is not fully understood and further research is needed, ideally to find therapies targeting both isoforms.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.