EGFR Variant Heterogeneity in Glioblastoma Resolved through Single-Nucleus Sequencing

Glioblastoma (GBM), also known as diffuse astrocytoma WHO grade IV, is the most common primary brain tumor in adults and a leading example of a tumor containing a high degree of cellular heterogeneity (1). Recurrent genetic alterations in GBM have been identified (2–4), with frequent amplification of several receptor tyrosine kinases (RTK) including EGFR, PDGFRA, and MET. Genomic heterogeneity in tumors of the central nervous system, especially GBM, is an emerging problem with substantial ramifications for treatment strategies (5, 6). One of the most prominent examples is the mosaic pattern of different RTK amplifications present in different subpopulations of tumor cells (7, 8).

As the most frequent RTK amplification in glioblastoma, EGFR is focally amplified in about half of all primary glioblastomas and in 95% of the classic subtype (9). Recent reports by The Cancer Genome Atlas (TCGA) revealed unprecedented complexity at the focally amplified EGFR locus within individual tumors, including clusters of chromosomal breakpoints (10, 11), and the expression of multiple EGFR variants and extracellular domain missense mutations (4, 12). However, outside of the EGFRvIII variant and carboxyl-terminal truncations, it is currently unknown whether many such variants are oncogenic and how such variants coexist at a single-cell level (13–15). Such patterns are important with regard to whether aberrations may functionally cooperate within individual cells or between cells with differing mutations. The functional importance of such distinctions is underscored by recent work showing the differential effects and interactions of variant receptors on EGFR activation within individual cells (16). Furthermore, extrachromosomal amplicons containing EGFRvIII (a constitutively active form of EGFR resulting from the deletion of exons 2–7) were recently shown to be dynamically lost and/or silenced in response to EGFR inhibitor treatment and then to re-emerge following removal of this selective pressure (5). The expression of EGFRvIII is also known to be regionally heterogeneous by IHC and RT-PCR, and it has been proposed that there may be functional synergy between tumor cells expressing the oncogenic variant (EGFRvIII) and those with the slower-growing counterpart (EGFRwt; refs. 5, 14, 17). However, whether this heterogeneity is due to differences in the genomes of single cells has not been fully assessed.

Traditional statistical analysis of bulk average sequencing data cannot resolve the complexity at the EGFR locus because both the high copy number and the nonmitotic segregation of EGFR variants present in extrachromosomal amplicons prevent accurate clonality inference (18, 19). Although FISH can assess the amplifications of different RTK genes (EGFR, PDGFRA, and MET), it is difficult to distinguish between multiple variants of these oncogenes by FISH, as it is generally limited to the characterization of known alterations spanning large genomic regions, and cannot probe single-nucleotide variants (SNV), short deletions, or discover de novo alterations.

Single-cell sequencing has great potential to dissect coexisting genomic alterations. Here, we describe a population-based framework for single-cell genomic analysis and demonstrate its potential by dissecting the subclonal populations in two primary glioblastomas with focal EGFR amplifications. Our methodology uses a novel approach for confident inference of LOH and enables unambiguous characterization of the subclonal tumor cell populations with implications for EGFR-targeted therapies in brain cancer and other cancers.

To characterize the frequency of compound alterations of EGFR in glioblastoma, we examined the RNA-sequencing data of 76 cases of glioblastoma from TCGA with focally amplified EGFR (4, 12), and observed that 71% (54 of 76) had expression of wild-type (WT) EGFR along with at least one aberrant EGFR variant, with 30% (23 of 76) expressing two or more variants. In 34% (28 of 76) of tumors, missense mutations involving the EGFR extracellular domain coexisted with structural variants. We then performed an exhaustive search of structural rearrangements of EGFR from whole-genome sequencing data in the 25 publicly available cases. Joint analysis of the DNA and RNA-sequencing data not only confirmed the presence of multiple variant EGFR transcripts, but also further revealed that the transcript encoding identical oncogenic EGFR variants could have been generated by distinct alterations at the DNA level.

As an example, our analysis suggested that TCGA-19-2624 had expression of wild-type EGFR together with both EGFRvIII (deletion of exons 2–7) and EGFRvII (deletion of exons 14–15). At the DNA level, two distinct genomic rearrangements were found to give rise to EGFRvIII, and two other events were found to generate EGFRvII (Fig. 1A). Coexpression of EGFRvII and EGFRvIII was previously described (20), but the frequency of variant reads suggested that EGFRvII is more prevalent in this sample. In addition, an extracellular domain missense mutation, G63K, was also observed at the DNA and RNA level. On the basis of the five EGFR aberrations alone, TCGA-19-2624 could exist as 32 possible combinations of cellular clones (2⁵). Similar scenarios were discovered in additional samples (Supplementary Fig. S1A–S1D and Supplementary Table S1), with an overall prevalence of 44% (7 out of 16 with focal EGFR amplifications) harboring multiple EGFR mutations, including both rearrangements and single-nucleotide substitutions. Bulk analyses cannot resolve whether such heterogeneity commonly occurs in a single tumor clone, mutually exclusive clones, or other complex combinations of multiple, yet functionally cooperative variants (Fig. 1B). To answer this question, we sought to assess the clonal structure of EGFR aberrations by single-nucleus sequencing (SNS).

We selected for study two primary glioblastomas (BT325 and BT340) with focal EGFR amplification present as extrachromosomal amplicons based on prior clinical analysis. Bulk tumor and SNS were performed to examine tumor heterogeneity within a single sampled region taken from each tumor (∼1 cm³ in size).

The experimental and analytical framework is summarized in Fig. 2. Single nuclei (SN) from tumor resections were isolated by flow cytometry and a limited multiple displacement amplification (MDA)–based whole-genome amplification was performed using a protocol that had been optimized for even genome coverage (see Supplementary Methods). Approximately 100 multiplexed SNS libraries were constructed and sequenced to low depth (∼0.01×) for quality evaluation (Supplementary Methods and Supplementary Fig. S2A and S2B). The quality assessment of single-cell DNA libraries requires an entirely new system of metrics than those used for bulk genomic DNA sequencing due to inherent challenges associated with starting with only two alleles of template DNA (21–23). Here, we described two metrics, both of which can be estimated from low-pass sequencing (median depth, ∼0.01×): one directly measures the amplitude of overamplification (whole-genome amplification pileup metric) and the other measures the unevenness of amplification (modified from Zong and colleagues; ref. 21). Both metrics showed consistent results (Supplementary Fig. S2C).

The 50 to 60 best SN libraries were selected and sequenced to a combined depth of 30×. In parallel, whole-genome sequencing was performed on dissociated nuclear DNA from the tumor to serve as the bulk tumor reference and on matched blood DNA as germline reference. We performed joint analysis of SN and bulk DNA sequencing data to infer and guide the clonality estimates of somatic mutations and the subclonal tumor populations. In comparison with prior studies, which examined only the single-cell genomes, performing joint analysis allowed us to more rigorously establish and confirm the presence of subclonal alterations using methods that have been previously established for bulk genome sequencing and validate findings from SNS.

Sequencing of the bulk BT325 tumor revealed the presence of a 935-kb amplicon containing four genes, including EGFR, with an additional deletion corresponding to EGFRvIII (Fig. 3A and Supplementary Fig. S3A). Bulk tumor analysis suggested an approximate 50% reduction in read density within this region, consistent with a mixture of amplified wild-type EGFR and EGFRvIII deletion. A total of 60 SN libraries from BT325 were sequenced and all the tumor-derived nuclei harbored the EGFR amplification as well as the vIII deletion (Fig. 3B). The common EGFR amplicon boundaries and EGFRvIII deletion breakpoints shared among tumor nuclei suggested that the deletion occurred after the wild-type EGFR amplification. The total EGFR copy number and the wild-type EGFR copy number estimated from the read depth demonstrated that each individual cell harbored amplifications of both wild-type EGFR as well as EGFRvIII, albeit at varying levels. To further validate our findings, we designed an oligonucleotide-based multiprobe FISH assay that specifically discriminates between the exons deleted in the EGFRvIII variant and those retained. Application of this novel FISH approach to BT325 validated the large differences in total EGFR copy number, including the presence of variation in levels of EGFRvIII among individual cells (Fig. 3C). Our observation of highly variable levels of EGFRvIII at the single-cell DNA level with respect to wild-type EGFR further supports a mechanism for local selective pressure of EGFR variants that may be associated with nutrient availability and diversification.

We next performed sequencing of the BT340 tumor. This tumor harbored a 535-kb amplicon containing the EGFR locus (Supplementary Fig. S3B). Bulk sequencing showed a marked decrease in read density in the region of exons 14 and 15 indicative of EGFRvII (Fig. 3D, top). An additional region of decreased read density indicative of a 7.6-kb deletion that included alternative exon 16 was also found. Discordant read pairs in the bulk tumor sequencing data corroborated the presence of these deletions and defined their precise boundaries as being distinct (Supplementary Fig. S3C). We denoted this novel deletion variant as EGFRvII-extended.

Examination of SNS libraries showed three distinct cell populations with respect to exons 14 to 16 of the EGFR locus (Fig. 3D). EGFRvII and EGFRvII-extended segregated into two populations based on read depth and discordant read pairs corresponding to the rearrangements (Fig. 3D, cell A and cell B). Cells that did not harbor any discordant read pairs supporting either EGFRvII variant were inferred to retain these exons (Fig. 3D, cell C). The initial amplicon was derived from a single copy of the wild-type allele (allele A) as assessed by germline SNVs; however, the exact number of copies per amplicon cannot be determined.

We further characterized the copy number of different EGFR variants within each individual cell. The total EGFR copy number was estimated from the read depths in exons 1 to 13 and in exons 17 to 28, and the wild-type EGFR copy number was estimated from the read depth in exons 14 and 15. In 35 of the 48 SN libraries, we were able to confidently assess the absolute number of EGFR copies. The total EGFR copy number showed an extremely broad range, from five copies to more than 200 copies in different nuclei, and was consistent with our FISH analysis (Fig. 3E, gray bars and Supplementary Fig. S3D). From 35 nuclei, we identified 26 containing EGFRvII and three EGFRvII-extended with no overlap. Further examination of the other six nuclei within the non-vII cluster revealed amplified EGFR, with three nuclei showing discordant read pairs spanning C-terminal breakpoints that occur in the presence of wild-type EGFR.

The common genomic boundaries of the overall EGFR amplicon in all three populations in BT340 suggested that the clonal diversification occurred after the common ancestral EGFR amplification. Although the vII-extended and the vII deletion existed in distinct cell populations, their clonal association could not be unambiguously determined (e.g., the larger vII-extended deletion may theoretically be derived from the smaller vII). We therefore sought to characterize the clonal hierarchy based on coexisting alterations unrelated to EGFR to infer the clonal origin of the two EGFR variants.

SNS allows one to reconstruct the lineage evolution history in tumors using copy-number profiling and single-nucleotide variant clustering (24–27). However, the generic clustering approach previously used in these studies cannot identify the exact clonal or subclonal alterations that differentiate subclonal populations as in bulk analysis (19, 28), and can be confounded by artifacts that are introduced during whole-genome amplification (22). To circumvent this limitation, we developed and applied an integrative approach to identify distinct tumor populations based on joint detection of clonal and subclonal events from the bulk tumor and SN genomes. Of particular importance, by requiring true variants to be present in at least two SN libraries, or with evidence in the bulk tumor DNA, we were able to exclude the majority of random MDA artifacts generated within individual libraries, and, at the same time, increase detection sensitivity for subclonal events that are present in the bulk library at low allelic frequencies.

We first focused on clonal events inferred from bulk tumor analysis that could be used as markers for SNS analysis. In both BT340 and BT325, we identified characteristic glioblastoma somatic copy-number alterations (SCNA), including homozygous loss of CDKN2A/CDKN2B, hemizygous loss of chromosome 10, as well as gains of chromosome 7 (Fig. 4A and B and Supplementary Fig. S4A and S4B). BT325 was found to have a common activating mutation within the promoter region of TERT (chr 5:1295228, G→A) that is observed in >80% of glioblastomas (4, 29–31). In BT340, we observed a nonreciprocal translocation of a putative super enhancer within chromosome 10 (32) to the TERT promoter that could be an alternative mechanism for activating TERT in glioblastoma (Fig. 4C).

We then characterized the SNS libraries of BT340 according to the clonal deletion events. In regions of hemizygous deletion or LOH, SN should harbor only one consensus haplotype. We used this feature to cluster individual nuclei based on the differences in their observed genotypes in LOH regions and infer the consensus haplotype within these regions (Supplementary Fig. S5A–S5H and Supplementary Methods). Using the germline heterozygous sites to infer the haplotype information from chromosome 10, we clustered the 48 nuclei derived from BT340 into two groups (Fig. 4D). The 44 tumor-derived nuclei (black circles) showed minimal presence of the deleted haplotype (whole-genome amplification and sequencing errors), and the four normal diploid nuclei (open circles) exhibited a roughly 1:1 ratio between the deleted and the retained haplotypes. The clustering of nuclei was further confirmed by haplotype analysis in hemizygously deleted regions in 9p and from the average copy ratio of chromosome 7, which is polysomic in tumor but disomic in normal nuclei (Supplementary Fig. S6A–S6C).

Haplotyping of 60 SNS libraries derived from BT325 based on deleted regions in chromosome 10 and chromosome 9p identified 56 libraries as tumor derived (Supplementary Fig. S7A–S7C, middle). Of the other four libraries, two showed an approximate 1:1 ratio between the deleted and the retained haplotypes, consistent with having been derived from normal diploid cells (open circles); the other two exhibited a 1:2 ratio between the deleted and the retained haplotype, implying a 1:1 mixture of tumor and normal nuclei (Supplementary Fig. S7D and S7E).

The haplotype-based method for clustering SNS libraries allowed for unprecedented sensitivity and specificity by using a common event in cancer (loss of heterogeneity) rather than sequencing read depths that can be grossly distorted by whole-genome amplification. Moreover, haplotype-based clustering uses the digital nature of single-cell genotypes (in contrast with allele frequencies in bulk sequencing data) and the clonal structure of tumor populations to achieve high accuracy at low sequencing depths (∼1×) even for libraries with high amplification bias.

Chromothripsis has been observed to be a frequent event in glioblastoma and often results in the disruption of regions that undergo recurrent somatic copy-number alterations (10, 33). In both BT325 and BT340, we observed that homozygous inactivation of CDKN2A/CDKN2B resulted from two events, one a simple deletion and the other a complex genomic rearrangement event resembling chromothripsis. In BT325, the chromothripsis event involved chromosomes 1, 3, 8, and 9 and resulted in three derivative chromosomes, one of which was subsequently amplified (Supplementary Fig. S7F). In a 1.25-Mb region encompassing the CDKN2A locus, we observed alternating copy-number events between three states (0, 1, and 2) from the bulk read depth characteristic of chromothripsis. The SNS coverage captured the retained regions after chromothripsis, and the signal was not influenced by the background contribution of contaminating normal cells. In BT340, a complex double-minute amplification involving multiple segments of chromosome 6 and chromosome 9 resulted in homozygous inactivation of CDKN2A/2B (Supplementary Fig. S8A–S8E).

Clonality analysis of SCNA events from bulk tumor sequencing of BT325 suggested three subclonal deletions within chromosomes 1, 4, and 6 at a copy ratio of approximately 0.61 (Supplementary Fig. S4A and S4B). Zooming into these regions, we observed alternating copy-number alterations involving chromosomes 1, 4, and 6 that are characteristic of chromothripsis (Supplementary Fig. S9A–S9E). Of the 56 tumor-derived nuclei, 46 nuclei (red circles) harbored hemizygous deletions in all the chromothriptic regions, whereas six nuclei (blue circles) were heterozygous at these loci (Supplementary Fig. S9B). Four additional SNS libraries had a haplotype ratio suggesting a 1:1 mixture of the tumor clones (magenta circles). This “all-or-none” scenario of all deletions clearly established that these events had been acquired within a single event as implied in chromothripsis.

On the basis of the presence of shared clonal features such as EGFRvIII, TERT promoter mutations, and loss of chromosome 10, but the absence of the chromosome 1, 4, and 6 chromothriptic event in six nuclei, we infer that the minor subclone split from the major subclone at this point (Supplementary Fig. S9F).

Analysis of SCNAs from the BT340 bulk tumor also showed copy ratios of approximately 0.20, 0.64, and 0.86 corresponding to noninteger copy numbers, indicating that such deletions were likely subclonal (Fig. 4A and B).

Regions of subclonal LOH were observed in chromosome 6p11.2 from 53.6 to 57 Mb, in 6q12 from 63.3 to 66.1 Mb, and in 9p21.1 from 27.1 to 30.9 Mb in the bulk tumor sequencing (Fig. 5A and Supplementary Fig. S6D). Clustering analysis of the 44 tumor-derived nuclei at these regions of chromosomes 6 and 9 resulted in two unique populations (Fig. 5B and Supplementary Fig. S6D). Cluster 1 (pink circles), consisting of 35 nuclei, was characterized by hemizygous deletions in 6p11.2 and 6q12 and homozygous deletions in 9p21.1 (27.1–30.9 Mb, 31.8–32 Mb). The second cluster of nine nuclei (green circles) was characterized by the retention of segments from 6p11.2 and 6q12 and hemizygous deletion of segments from 9p21.1. The haplotype ratio between the deleted (D) and the retained (R) in hemizygously deleted regions in 6p11.2 confirmed the hemizygous loss of these regions in cluster 1.

We then compared the subclonal copy numbers statistically inferred from the bulk analysis to the clonal fractions from single-cell analysis. The bulk average allelic copy number of the subclonal deletion in 9p21.1 is approximately 0.2 after correcting for tumor purity (90%) and ploidy (2.1). Clustering of single tumor nuclei found that 20% (9 of 44) retained these regions and 80% (35 of 44) harbored the deletion (Fig. 5C). The same agreement was observed for regions in chromosome 6 in which the bulk average copy number is approximately 0.70. Finally, the copy ratio of 0.86 from the bulk analysis, corresponding to a subclonal deletion of chromosome 16 at approximately 19% clonal fraction, was associated to 7 out of 35 nuclei within cluster 1 with hemizygous loss (Supplementary Fig. S6E). The excellent concordance between the bulk clonality analysis and the single-cell analysis highlights the power of our population-based approach to robustly and efficiently resolve heterogeneous subclonal populations even at low clonal frequencies. The consistency between the clonality analysis of bulk sequencing and population analysis of SN not only validates the statistical rationale behind bulk tumor analysis (34), but further sets a standard paradigm for single-cell–based tumor analysis. Furthermore, single-cell analysis clearly resolved that the subclonal chromosome 16 loss occurred in a different population than the subclone that retained segments from chromosomes 6 and 9, although their clonal fractions are similar (∼19% and ∼20%). Thus single-cell sequencing provides an unambiguous solution to associate independent subclonal alterations with distinct subclonal populations.

Another application of single-cell sequencing in dissecting heterogeneous tumors is demonstrated by correlating alterations in amplified loci with other subclonal events, which is impossible to be inferred from bulk tumor sequencing data. Nucleus-by-nucleus comparison of EGFR mutation status to the clonal hierarchy constructed from subclonal deletions above revealed that EGFRvII and EGFRvII-extended were exclusively associated with only one of the two subclonal populations: EGFRvII was present only in cluster 1 with 6p11.2 and 6q12 deletions, whereas EGFRvII-extended was present only in cluster 2 that retained these chromosomal segments (Fig. 5D). In addition, cluster 2 also contained cells with a known oncogenic C-terminal truncation (ex. 25–28; ref. 15) and a second C-terminal deletion (ex. 26–28), all of which occur in the context of wild-type EGFR amplification. Three cells from cluster 1 and two cells from cluster 2 were found not to contain focally amplified EGFR and, thus, could not be assigned as being a variant. In summary, 89% (39 of 44) of tumor cells were found to contain EGFR variants, with EGFRvII being the most common, followed by EGFRvII-ext., EGFR (del25–28), and EGFR (del25–26). Remarkably, all EGFR variants were mutually exclusive, suggesting that the EGFR mutations were most likely independently derived and expanded. In particular, the convergent evolution of EGFRvII variants suggested its oncogenic potential.

EGFRvII was initially described in GBM two decades ago, but was reported to be nononcogenic (35, 36). Subsequent studies identified the presence of EGFRvII coexisting with EGFRvIII and, thus, EGFRvIII was inferred to be the driving oncogenic variant (20). However, in TCGA-19-2624, the EGFRvII variant was found to be more prevalent than EGFRvIII at both the DNA and RNA levels, arguing in favor of its functional relevance. Furthermore, the convergence of independent clonal diversifications in BT340 resulting in EGFRvII and EGFRvII-extended further underscored its oncogenic potential. We therefore assessed the transformation potential of EGFRvII and the previously characterized EGFRvIII using two model systems (Fig. 6A).

Ba/F3 cells are dependent upon supplemental IL3 for survival but can be transformed to ligand independence by expression of oncogenic EGFR mutants (37, 38). Ectopic expression of EGFRvII and EGFRvIII permitted Ba/F3 cells to grow in the absence of IL3, whereas cells containing the empty parental vector or wild-type EGFR did not survive (Fig. 6B). Constitutive phosphorylation of EGFR^Y1068 and AKT^S473 was observed in serum-starved EGFRvII- and EGFRvIII-expressing cells, demonstrating ligand-independent activation (Fig. 6C).

To evaluate the oncogenic potential of the EGFRvII variant in a genetically relevant model of GBM, we expressed EGFR variants in E14.5 neural stem cells derived from Ink4a/Arf^−/− mice (39). All variants demonstrated constitutive phosphorylation of EGFR as compared with wild-type EGFR-expressing cells (Fig. 6D). Consistent with published reports, the expression of EGFRvIII induced tumor growth when subcutaneously grafted in nude mice (Fig. 6E; ref. 13). Expression of EGFRvII was also tumorigenic, albeit tumors developed at a slower rate (11 weeks for EGFRvII compared with 6 weeks for EGFRvIII; Fig. 6E and F). Histologically, EGFRvIII and EGFRvII both generated high-grade gliomas, but we observed that the EGFRvII tumors exhibited oligodendroglial features and a corresponding lower rate of proliferation. These findings suggest that these similar mutations produce different functional effects.

To examine the drug sensitivities of EGFRvII and EGFRvIII, transformed Ba/F3 cells were tested for sensitivity to a panel of EGFR inhibitors. Ba/F3 cells dependent on EGFRvII were killed by erlotinib, lapatinib, afatinib, and neratinib (Supplementary Fig. S10A–S10D), with the irreversible inhibitors afatinib and neratinib exhibiting the lowest IC₅₀ values for both EGFRvII- and EGFRvIII-expressing cells (Supplementary Fig. S10C and S10D). These results are similar to previous reports for EGFRvIII-expressing Ba/F3 cells (38), and suggested that EGFRvII may be a therapeutic target for kinase inhibition.

Our identification of subclone-specific rearrangements within the EGFR gene expands upon the recent discoveries of heterogeneous genomic alterations in oncogenic kinase drivers in glioblastoma (7, 8, 14, 40). These results further suggest that a substantial fraction of the multiple EGFR variants reported in the TCGA glioblastoma collection may represent distinct subclonal populations within the tumor.

High-resolution genomic approaches have begun to reveal genomic complexity and heterogeneity to a degree not previously understood. These methods include clonality analysis of bulk tumor sequencing (18, 28, 41), regional sampling of tumors (6, 42, 43), and SNS (24–26). Among these techniques, SNS has the unique power to resolve cellular heterogeneity, but can be technically challenging and difficult to interpret. Our methods for LOH detection by allele-specific delineation of haplotypes address limitations in genome coverage due to MDA and expand upon the applications of single-cell DNA analysis. As single-cell sequencing studies expand in size and scope, our methods provide the basis for maximizing sequencing information, given the numbers of cells sequenced. Clonality inference by combining bulk tumor analysis with single-cell population analysis further provides a unique strategy for resolving clonal evolution and coassociation of independent genomic alterations at the single-cell level. The combination of SNS and joint calling with bulk sequencing allowed us to detect minor subclonal populations that were present in as low as 10% of the tumor at relatively low sequencing depth.

Genomic diversification in cancer is known to contribute to targeted and cytotoxic therapy failure and emergence of resistance (5, 44, 45). Large-scale genomic studies of cancer have outlined sets of alterations that affect growth-promoting and tumor-suppressing pathways that are necessary to drive tumorigenesis. In glioblastoma, the RB, TP53, and PTEN tumor suppressors are frequently inactivated directly or indirectly through aberrations in pathway factors such CDKN2A/CDKN2B or MDM2 (4). These genomic alterations serve as the foundation upon which RTKs like EGFR, MET, and PDGFRA can become amplified and evolve as the tumor cells begin to compete with each other for growth advantage. As FISH studies have demonstrated, RTKs within a tumor can be amplified alone or in combination (Fig. 7, left; refs. 7, 8). Presumably, clones with diverse growth signaling pathways will gain an advantage and expand as competition ensues.

Here, our study reveals an additional paradigm for the diversification of growth-signaling pathways through the acquisition of multiple variants of a single RTK gene (Fig. 7, right). The convergence of subclone-specific EGFRvII deletions reflects the dynamic evolution of amplified EGFR copies and pinpoints its oncogenic potential. We further demonstrated that overexpression of EGFRvII does cause transformation in both Ba/F3 and Ink4a/Arf^−/− neural stem cell models, a variant previously thought not to be oncogenic because it co-occurred with EGFRvIII (20, 35). TCGA RNA sequencing even revealed the expression of EGFRvII to be present in 9% (7 of 76) of cases with focally amplified EGFR (4). Our functional characterizations found that constitutive expression of EGFRvII results in downstream activation of AKT signaling, consistent with that of EGFRvIII, but not with enhanced ERK activation. It is possible that EGFRvII activates an alternative pathway such as STAT3/5 and poses a more direct means of inducing transcriptional changes (16). Therefore, the presence of multiple EGFR variants in a single tumor highlights the intratumoral heterogeneity of glioblastoma conferred by the plasticity of EGFR amplicons.

Heterogeneous expression of oncogenic EGFR mutants and cooperativity of EGFR with itself or other RTKs have each been hypothesized to contribute to tumor growth and resistance to therapy (14, 16, 40). It is conceivable that heterogeneity and cooperativity mechanisms play a role in RTK diversification and will affect the efficacy of targeted therapies in cancer, and, in particular, may explain the lack of response in glioblastoma. Further complicating this scenario is recent work demonstrating how EGFRvIII-containing cells can evade EGFR inhibitors through amplicon loss and downregulation of the oncogenic receptor (5). It is striking that small-molecule inhibitors can be so successful in non–small-cell lung cancers, but remain ineffective in glioblastomas, despite EGFR being functionally significant (2, 3, 46–49). This work suggests that combining multiple EGFR inhibitors that act through different mechanisms on the receptor could be required in patients with multiple variants in nonoverlapping subclones and also potentially prevent the activation of a resistance pathway. Such concepts deserve further evaluation in preclinical models.

As cancer genomic research approaches the discovery of the plurality of somatic mutations, it becomes crucial to understand the distribution pattern of somatic mutations among individual tumor cells and clones and how these mutations work together to drive tumor fitness (50). Single-cell analysis has the potential to greatly contribute to our understanding of co-dependency of mutations and tumor heterogeneity. The input requirements of this approach also make it feasible to deliver great insights into biopsies, necrotic tissues, or tumors with great immune infiltration to delineate the full spectrum of cell-specific oncogenic drivers. Once a thorough understanding of the tumor state can be integrated, it is then possible to select combinations of drugs to prevent resistance due to emergence of preexisting subclonal populations of tumor cells not detected in bulk tumor sequencing.

Anonymized tumor and matching blood specimens were collected following informed, written consent and processed under protocols approved by the Institutional Review Boards of the Dana-Farber Cancer Institute and the Broad Institute (protocol 11-433). Two cases of de novo GBM, designated as brain tumor BT325 and BT340, were mechanically dissociated into pools of nuclei as previously described (24). A portion of the dissociated material was selected to represent bulk tumor. Genomic DNA was isolated from the pooled nuclei and matched blood using DNeasy columns (Qiagen).

SN were sorted directly into 96-well PCR plates containing 5 μL of phosphate-buffered saline using the BD FACSAria IIu flow cytometer calibrated to deposit SN into each well. A 4′,6-diamidino-2-phenylindole (DAPI)-based DNA stain was used to select intact nuclei. Empty wells were included on every plate to control for background MDA amplification. See Supplementary Methods for additional details.

We generated SN DNA libraries from approximately 100 SN with the aim of selecting the best libraries from each tumor for deeper sequencing and analysis. Quality assessment was done by multiplexed sequencing of SN libraries to an average depth of 0.01× per library. Two criteria were used to select the best SNS libraries for deeper sequencing.

First, we removed libraries with more than 20% read sequences not aligned to the human reference. We then evaluated the evenness of the remaining MDA libraries by two metrics. The first metric was adapted from the method recently described (21) to measure the evenness of genome-wide coverage by evaluating the distribution of reads in a DNA library across different genomic loci stratified by the sequencing coverage at each locus. The second metric, termed the pileup metric, estimates the percentage of sequenced reads that were accumulated in local regions at densities 10 times higher than the chromosomal average. The spacing between each read and the adjacent reads (left and right) was used to infer if the read was in a pileup region. Reads in pileup regions were then downsampled to the average density estimated from the flanking nonpileup regions. The downsampling was iterated until no pileup region exceeding a certain size (e.g., 100 reads) was remaining. The pileup metric was calculated from the percentage of reads removed during the downsampling procedure.

The two metrics showed a strong correlation between each other (Supplementary Fig. S2C): both can be used to estimate of the evenness of coverage for SNS libraries to guide the optimization of SNS workflow or select better libraries for deeper sequencing. In addition, removing the pileup reads can smooth the sequence coverage and reveal true copy-number profiles for each SNS library, which can be used to identify (tumor) samples of interest.

Whole-genome amplification by MDA is known to generate artifacts, including single-base errors and random chimera, which can be indistinguishable from true somatic mutations. Given these considerations, we adopted a two-step strategy to analyze somatic mutations by SNS of many tumor nuclei. First, we detected all clonal and subclonal somatic variants, including single-nucleotide variants and chromosomal rearrangements, from the pool of SNS complemented with bulk whole-genome sequencing. Second, we searched for the detected variants within each SNS and in the bulk tumor to infer the clonality of these variants and construct the clonal hierarchy of these variants. See expanded methods for additional details.

Whole-genome amplification by MDA from single-nuclear DNA can generate multiple distortions to the original DNA copy-number profile, which pose a considerable challenge to the detection of true somatic copy-number events in SN libraries. We, therefore, did not use the read depth signal from SN libraries to detect copy-number breakpoints but only inferred copy-number breakpoints from read depths in the bulk tumor against the germline reference using SegSeq (51). At 10× sequencing depth, the read depth signal from the bulk tumor should be sufficient to detect somatic copy-number alterations of 10 kb or longer at or above 15% clonal fraction, which is comparable with the smallest subclone we can identify from 50 to 60 single cells.

As the read depth signal often resulted in oversegmentation, we merged the copy-number breakpoints inferred from the bulk tumor read depth with the list of nonreciprocal chromosomal translocation breakpoints detected from the pool of bulk tumor and SN libraries. We then identified true chromosomal breaks as supported both by read depth difference across the breakpoint and by the presence of somatic rearrangement signatures on the higher copy-number side. This final set of chromosomal breakpoints was then used to segment the copy-number profile for the bulk tumor as well as for each SN.

A unique feature of many tumor genomes is the LOH of large chromosomal regions as well as smaller homozygous deletions, leading to the elimination of one or both germline haplotypes (refs. 52, 53; Supplementary Fig. S4A). For cells with LOH in a region of the genome, the retained haplotype, or allelotype (54), can be directly determined from the sequencing data (Supplementary Fig. S5B). If whole-genome amplification from single cells were complete and represented all alleles, then the haploid genotype of cells with LOH could be readily distinguished from the heterozygous genotype of diploid cells (Supplementary Fig. S5B, left). However, single-cell sequencing is inevitably incomplete and not all alleles are represented in the resulting data (Supplementary Fig. S5B, right). Single-base errors introduced during MDA can further confound the inference of true genotype.

We devised a clustering-based approach to dissect the cell populations harboring heterozygous and/or homozygous genotypes in each genomic region. First, the pairwise distance between single-cell genotypes was calculated from the fraction of germline heterozygous sites covered in both libraries with different genotypes (Supplementary Fig. S5C); this distance is large if cells are discordant for heterozygosity/homozygosity, and small if they are both homozygous (concordant). We then performed clustering analysis of cellular populations to identify cells with LOH in a given region, which formed a very tight cluster (Supplementary Fig. S5C, right). Supplementary Figure S5D shows the two clusters formed by 60 SNS libraries from BT325 identified by LOH in chromosome 10, one cluster displaying a unique haploid genotype (LOH), and a second cluster that was heterozygous.

By pooling the genotype data of LOH cells in cluster 1, we can determine both the retained haplotype “R” and the deleted haplotype “D” and assign individual or combined “RD” haplotypes to each nucleus at all covered germline heterozygous sites (Supplementary Fig. S5E). This analysis showed that cluster 1 contains BT325 nuclei with only allele “R,” whereas cluster 2 contains an admixture of alleles “R” and “D” on chromosome 10 (Supplementary Fig. S5F). The ratio of the deleted to retained haplotypes is approximately 0 in LOH nuclei, and close to 1 in diploid nuclei. Interestingly, the D:R ratio is near 0.5 in another population of nuclei (Supplementary Fig. S5G).

We inferred that the population with a D:R ratio of 0.5 represents libraries derived from a 1:1 mixture of tumor and normal nuclei. The normal/tumor mixture was later confirmed by the presence of clonal drivers at subclonal quantities in these samples (Supplementary Fig. S7). Therefore, the haplotype analysis also provides an accurate method to identify and characterize cellular mixtures in presumably single-nuclear samples.

FISH was performed on 5-μm formalin-fixed paraffin-embedded (FFPE) tissues, which were baked at 60°C for at least 2 hours, then deparaffinized and digested using methods described previously (55). See Supplementary Methods for probes used in this study. Tissue sections and probes were codenatured, hybridized at least 16 hours at 37°C in a darkened humid chamber, washed in 2× saline-sodium citrate (SSC) at 70°C for 10 minutes, rinsed in room-temperature 2× SSC, and counterstained with DAPI (Abbott Molecular/Vysis, Inc.). Slides were imaged using an Olympus BX51 fluorescence microscope. Individual images were captured using an Applied Imaging system running CytoVision Genus version 3.9, and all aberrations detected by FISH were reviewed and confirmed by a cytogeneticist (A.H. Ligon).

Bacterial artificial chromosome (BAC) probes were all obtained from the CHORI BACPAC Resources Center (Oakland, CA) and included: clone RP11-815K24 (corresponding to EGFR, at 7p11.2), clone RP11-121K16 (corresponding to 6q13), and clone RP11-482I10 (mapping to 9p21.3). Commercially available probes (all obtained from Abbott Molecular) included CEP6 (6p11.1-q11), CEP7 (7p11.1-q11.1), and CEP9 (9p11-q11). Final homebrew probe concentrations were 100 ng/μL. Final concentrations for commercial reference probes followed the manufacturer's directions.

Specific detection of the EGFRvIII and EGFR wild-type alleles was performed using custom-designed pooled oligonucleotide probe sets (SureFISH), which tiled exons 2 to 7 (red probe, retained in wild-type EGFR) and exons 9 to 28 (green probe, used to enumerate copies of both wild-type EGFR and EGFRvIII). Hybridization was performed according to the manufacturer's protocol (Agilent Technologies) on 4-μm FFPE sections from BT325.

Wild-type EGFR and EGFRvIII p-BABE-puro vectors were prepared as previously described (38). QuikChange site-directed mutagenesis (Stratagene) was used for generating EGFR vII using wild-type EGFR as a template. All plasmids were confirmed by sequencing.

Ba/F3 cells were obtained, thawed, and maintained as previously described (38). Proliferation assays were performed as previously described (37). The Ba/F3 cells were confirmed to be the parental line through a cell viability test following removal of IL3.

Mouse neural stem cells were derived from E14.5 ganglionic eminencies of p16Ink4a:p19Arf-null embryos, genotyped by qPCR on genomic DNA as described (13, 39). Cells were transduced by EGFR wild-type, EGFRvII, or EGFRvIII p-BABE retrovirus and selected with puromycin after 3 days for two passages. Transduction of the EGFR variant was confirmed by immunoblot. For subcutaneous xenograft studies, EGFRvII- or EGFRvIII-transduced mouse neural stem cells (10⁶ cells) were injected in the left flank of nude mice with EGFR wild-type injected on the right as control. Five mice were injected for each EGFR variant. All animal experiments were performed in accordance with Dana-Farber Cancer Institute animal facility regulations and policy under the protocol #09-016.

For growth inhibition assays, Ba/F3 cells (10,000 cells) were plated in 180 μL of media in 96-well flat-bottom plates. Twenty-four hours after plating, drugs prepared in 20 μL of cell culture medium were added to cells. All drugs were purchased from LC Laboratories. The concentrations of erlotinib and lapatinib used for the assay ranged from 0.1 nmol/L to 10 μmol/L. The concentrations of afatinib and neratinib ranged from 0.001 to 100 nmol/L. The cells were incubated for another 72 hours and the viable cell numbers were measured using CellTiter Glo (Promega).

Cultures were serum-starved for 3 hours before EGF stimulation and harvesting. EGF (PeproTech) stimulations were performed using 25 ng/mL for 10 minutes. Anti-EGFR, phospho-EGFR (pY1068), AKT, phospho-AKT (pS473), ERK1/2, and phospho-ERK1/2 (pT202/Y204) antibodies were purchased from Cell Signaling Technology. Anti-vinculin was from Sigma-Aldrich.

M. Meyerson reports receiving a commercial research grant from Bayer and has ownership interest (including patents) in and is a consultant/advisory board member for Foundation Medicine. No potential conflicts of interest were disclosed by the other authors.

Conception and design: J.M. Francis, C.-Z. Zhang, C.L. Maire, J. Jung, V.A. Adalsteinsson, M. Meyerson, K.L. Ligon

Development of methodology: J.M. Francis, C.-Z. Zhang, C.L. Maire, J. Jung, B. Blumenstiel, K.L. Ligon

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J.M. Francis, C.-Z. Zhang, C.L. Maire, V.E. Manzo, S. Haidar, J.C. Love, K.L. Ligon

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J.M. Francis, C.-Z. Zhang, C.L. Maire, S. Haidar, C.S. Pedamallu, A.H. Ligon, K.L. Ligon

Writing, review, and/or revision of the manuscript: J.M. Francis, C.-Z. Zhang, C.L. Maire, J. Jung, A.H. Ligon, M. Meyerson, K.L. Ligon

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J.M. Francis, C.-Z. Zhang, H. Homer, S. Haidar, K.L. Ligon

Study supervision: J.M. Francis, C.-Z. Zhang, J.C. Love, M. Meyerson, K.L. Ligon

The authors thank Dr. Juliann Chmielecki for generating the EGFRvII expression construct and Drs. Scott Carter and Marcin Imielinski for helpful discussions. The authors also thank the Brigham and Women's Hospital Tissue Biorepository Core (Dr. William Richards, Director) and the Department of Pathology faculty and staff of the Dana-Farber Cancer Institute Tissue Bank for their help and dedication. The authors also thank Leslie Gaffney for help in preparing figures.

This work was supported by the Bridge Project, a collaboration between The David H. Koch Institute for Integrative Cancer Research at MIT and the Dana-Farber/Harvard Cancer Center (DF/HCC), and the National Brain Tumor Society (to J.C. Love, K.L. Ligon, and M. Meyerson). This work was supported in part by grants to K.L. Ligon from NCI R01CA170592, P01CA142536, and the Sontag Foundation, and to M. Meyerson from NCI R01CA116020.