Tumoral and immune heterogeneity in an anti-PD-1-responsive glioblastoma: a case study.

Clinical benefit of immune checkpoint blockade in glioblastoma (GBM) is rare, and we hypothesize that tumor clonal evolution and the immune microenvironment are key determinants of response. Here, we present a detailed molecular characterization of the intratumoral and immune heterogeneity in an IDH wild-type, MGMT-negative GBM patient who plausibly benefited from anti-PD-1 therapy with an unusually long 25-mo overall survival time. We leveraged multiplex immunohistochemistry, RNA-seq, and whole-exome data from the primary tumor and three resected regions of recurrent disease to survey regional tumor-immune interactions, genomic instability, mutation burden, and expression profiles. We found significant regional heterogeneity in the neoantigenic and immune landscape, with a differential T-cell signature among recurrent sectors, a uniform loss of focal amplifications in EGFR, and a novel subclonal EGFR mutation. Comparisons with recently reported correlates of checkpoint blockade in GBM and with TCGA-GBM revealed appreciable intratumoral heterogeneity that may have contributed to a differential PD-1 blockade response.

INTRODUCTION

Anti-PD-1 immune checkpoint blockade, which serves to bolster the patient's own antitumor immune response, represents a novel therapeutic strategy in many cancer types, in addition to traditional treatments. Although checkpoint inhibition has produced outstanding results in those patients who do respond to this treatment, response rates remain stubbornly low for many tumor types (Ribas and Wolchok 2018). Because tumor clonal evolution and the immune microenvironment may ultimately determine the potential benefit of this novel therapeutic strategy, there is a crucial need to understand the conditions under which PD-1 checkpoint blockade can produce a clinically meaningful antitumor response in cancers with poor clinical responses.

Glioblastoma (GBM) is a very aggressive and highly heterogeneous cancer type with a median patient survival time of ∼14 mo (Delgado-López and Corrales-García 2016; Nam and de Groot 2017). There is increasing interest in immunotherapeutic treatment options for GBM, although clinical trials have largely proven unsuccessful in improving survival outcomes to date (Thomas et al. 2012; Reardon et al. 2014; Chin et al. 2018; Cloughesy et al. 2019). Reports of successful checkpoint blockade in glioblastoma have been linked to hypermutation and mismatch repair deficiency (Erson-Omay et al. 2015; Bouffet et al. 2016; Johanns et al. 2016; Zhao et al. 2019). Although tumor mutation burden (TMB) is modestly correlated to immunotherapy response (Vogelstein et al. 2013; Campbell et al. 2017; Yarchoan et al. 2017), GBM tends to present with few somatic mutations relative to other cancer types (Alexandrov et al. 2013; Hodges et al. 2017). Recent efforts to characterize genomic correlates of checkpoint-blockade response in GBM include a large study of 66 nonhypermutated GBM patients who were treated with immune checkpoint inhibitors at recurrence. An extensive genomic characterization was carried out, with particular attention to 17 long-term responders (Zhao et al. 2019). A key conclusion was that responders were enriched in BRAF/PTPN11 mutations, MAPK alterations, an altered Treg signature, and a branched pattern of clonal evolution. Nonresponders, on the other hand, were characterized by PTEN mutations and a linear pattern of clonal evolution. Here we present a case of an anti-PD-1-responsive GBM patient with unusually long survival after removal of recurrent disease. We provide a deep molecular and immunophenotypic survey of the immune landscape and intratumoral heterogeneity of the tumor that is overall representative of a typical GBM, including genomic stability, unmethylated MGMT, wild-type IDH, and the absence of hypermutation.

Leveraging multiplex immunohistochemistry (mIHC), RNA-seq, and whole-exome sequencing (WES) from the primary tumor and several biopsies at recurrence, we characterize the intratumoral heterogeneity and immune landscape of an anti-PD-1-responsive patient. We carried out expression, neoantigen, copy-number variation (CNV), human leukocyte antigen (HLA), and T-cell receptor (TCR) profiling across the sectors. Furthermore, we projected our patient sector profiles onto 155 The Cancer Genome Atlas (TCGA)-glioblastoma (GBM) RNA-seq samples to scale and contextualize regional GBM expression variation. We also interrogated intratumoral heterogeneity of the patient's recurrent sectors with the anti-PD-1-response biomarkers identified by recent studies (Zhao et al. 2019).

RESULTS

Case Report

The patient was a 67-yr-old woman who presented with a generalized tonic-clonic seizure and underwent a gross total resection of a right temporal GBM, IDH wild-type, MGMT promoter methylation not detected. IDH1 and IDH2 wild-type status was determined via WES somatic mutation calling and a cancer hotspot genotyping panel. MGMT methylation status was determined via methylation specific real-time polymerase chain reaction (PCR). She received standard-of-care treatment consisting of radiotherapy with concurrent temozolomide, and after two cycles of adjuvant temozolomide, her tumor recurred. She underwent a reresection of her right temporal tumor with persistent medial tumor after surgery (Fig. 1A,B) and morphological characteristics of a gliosarcoma. Four specimens were collected for sequencing, mIHC, and further analysis: the primary tumor, and three distinct portions of recurrent tumor prior to nivolumab treatment, A (lateral), B (inferior), and C (medial) (Fig. 1B,C). Within 2 wk after reresection, she started PD-1 immune checkpoint blockade and received 26 cycles of nivolumab spanning for 12 mo of treatment until tumor progression. After seven cycles of nivolumab (3 mo after reresection) she was treated with bevacizumab, a VEGF-A inhibitor, for symptoms including unsteady gait, a partial right third nerve palsy, and a right upper quadrantanopsia, which worsened since surgery. She received 20 cycles of bevacizumab (9 mo of treatment) until progression, maintaining an ECOG performance status of 2 until then (Fig. 1A). According to MRI volumetric analysis, sectors A and B maintained no tumor growth after resection. The tumor volume of sector C was 0.33 at 5 mo, 0.65 mL at 10 mo, and 0.86 mL at 12 mo. Along with longitudinal imaging, this suggests that the slow progression of disease arose from residual tumor near the location of sector C, followed by expansion toward the inferior regions of the brain near sector B (Fig. 1B). She survived 25 mo after the initial diagnosis including 6 mo after the treatment was discontinued (Fig. 1A). Based on the empirical cumulative density of survival in 155 TCGA-GBM patients, the probability that our patient survived this long by chance alone is 15.89% (Fig. 1D).

Figure 1.

(A) Patient clinical timeline showing treatment, survival, and points of sample collection. (B) MRI imaging of tumors resected for sequencing and analysis, followed by longitudinal scans at 5, 10, and 12 mo post-resection of the recurrent tumors A–C. Tumor volumes are shown directly on the MRI snapshots. Slow progression arises more clearly at the 10 mo near the medial region and expanding toward the inferior regions of the brain. (C) Summary of the downstream experimental analysis procedures. (D) Survival of the patient projected onto the survival distribution of 155 TCGA-GBM samples.

Genomic Analyses

Somatic Driver Genes and Known Correlates with Response in Recurrent Sectors

We sought to identify functional contributors to tumor progression, as well as the clonal heterogeneity among recurrent sectors and compared to the primary resection. We searched for mutations occurring in genes that have been previously reported to be somatic cancer drivers in the COSMIC Cancer Gene Census (Tate et al. 2019). Additionally, we searched for somatic mutations (Fig. 2B) and CNVs (Fig. 2C) in the PTEN, BRAF, and the MAPK families of genes to contextualize the patient sectors within known correlates of checkpoint blockade response (Zhao et al. 2019). We detected a nonsynonymous PTEN mutation that was conserved throughout the primary and recurrent sectors. Further, sectors A and B shared a TP53 mutation (Table 1; Fig. 2B). The primary tumor and all recurrent sectors had a conserved mutation in EGFR (dbSNP ID: rs121913428, COSMIC ID: COSM18425) that overlaps with a copy-number amplification in the primary tumor, with a uniform loss of intensity in the recurrent sectors (Table 1; Fig. 2B,C). However, sector C had an additional EGFR mutation at the Chr 7:55201281_G > C locus that is not reported in COSMIC (v70) or dbSNP (Fig. 2B; Supplemental Fig. 1). Selected variants are detailed further in Supplemental Table 1.

Figure 2.

(A) (Top to bottom) Tumor mutation burden (TMB), tumor purity, and genomic instability measures for each sample. (B) Notable somatic mutations, and (C) copy-number alterations. (D) Variant allele frequencies (VAFs) of altered loci identified in each detected cluster. The gray line indicates the mean VAF for each cluster. Each cluster is an inferred tumor subclone. (E) Bell plot showing the predicted cellular fraction for each cluster per sample overlaid upon the clonal evolution model and inferred clonal ancestry. (F) Predicted cellular fraction of each clone per sample depicted as a cloud of cells. (G) Heatmap depicting somatic mutation overlap among tumor biopsies. (H) Phylogenetic tree of inferred tumor subclone ancestry. Nodes represent inferred tumor subclones. Each node is labeled with the samples in which that clone was detected.

Table 1.

Variant table

Gene	Chromo some	HGVS DNA reference	HGVS protein reference	Variant type	dbSNP ID	VAF	ClinVar ID	Sample	Comments
EGFR	Chr 7	NC 000023.7: g.55174015 G > A	NM 001346941: p.G452D	Nonsynonymous SNV	rs121913428	0.137	362954	Recurrent A
EGFR	Chr 7	NC 000023.7: g.55174015 G > A	NM 001346941: p.G452D	Nonsynonymous SNV	rs121913428	0.117	362954	Recurrent B
EGFR	Chr 7	NC 000023.7: g.55174015 G > A	NM 001346941: p.G452D	Nonsynonymous SNV	rs121913428	0.106	362954	Recurrent C
EGFR	Chr 7	NC 000023.7: g.55201281 G > C	NM 001346941: p.D747H	Nonsynonymous SNV	NA	0.073	NA	Recurrent C	Validated with Sanger
EGFR	Chr 7	NC 000023.7: g.55174015 G > A	NM 001346941: p.G452D	Nonsynonymous SNV	rs121913428	0.837	362954	Primary
PTEN	Chr 10	NC 000023.10: g.87933214 T > TA	NM 000314: p.D153Rfs*26	Frameshift insertion	NA	0.28	NA	Recurrent A
PTEN	Chr 10	NC 000023.10: g.87933214 T > TA	NM 000314: p.D153Rfs*26	Frameshift insertion	NA	0.216	NA	Recurrent B
PTEN	Chr 10	NC 000023.10: g.87933214 T > TA	NM 000314: p.D153Rfs*26	Frameshift insertion	NA	0.056	NA	Recurrent C
PTEN	Chr 10	NC 000023.10: g.87933214 T > TA	NM 000314: p.D153Rfs*26	Frameshift insertion	NA	0.492	NA	Primary
TP53	Chr 17	NC 000023.17: g.7675088 C > T	NM 001126115: p.R43H	Nonsynonymous SNV	rs28934578	0.258	27413	Recurrent A
TP53	Chr 17	NC 000023.17: g.7676041 G > A	NM 001126118: p.R71C	Nonsynonymous SNV	NA	0.137	151920	Recurrent A
TP53	Chr 17	NC 000023.17: g.7675088 C > T	NM 001126115: p.R43H	Nonsynonymous SNV	rs28934578	0.254	27413	Recurrent B
TP53	Chr 17	NC 000023.17: g.7676041 G > A	NM 001126118: p.R71C	Nonsynonymous SNV	NA	0.103	151920	Recurrent B

Because mutation burden is a reported correlate of checkpoint blockade response (Hodges et al. 2017; Yarchoan et al. 2017), we measured the mutation burden of each recurrent sector using somatic mutations from WES. None of the sectors were hypermutated, with recurrent sector A showing the highest mutation burden and recurrent sector C showing the lowest (Fig. 2A).

Patient Sectors Show a Branched Pattern of Clonal Evolution with an Outlier Recurrence

We computed the percentage of overlapping somatic variants among sectors and observed relatively higher overlap between sectors A and B and very low overlap with sector C (Fig. 2G). These results were additionally corroborated using the SciClone clonal inference package on the somatic mutations to determine clonal structure, followed by clonal ordering and tree reconstruction using ClonEvol (Miller et al. 2014; Dang et al. 2017). Four distinct clones were detected, although clone 4 is a predicted descendant of clone 3 and was detected with a very small, residual cellular fraction among the recurrent sectors only. Clones 1 and 3 are dominant in the primary tumor and recurrent sector C. Clone 2, a predicted descendant of clone 1, was detected only in sectors A and B. Clone 3 was detected in all samples, although with less cellular fraction in the primary tumor. Furthermore, the phylogenetic tree shows a branched pattern of clonal evolution (Fig. 2C,H). The increased cellular fraction of clone 3 in the recurrent sectors, along with the residual presence of its descendant, clone 4, and the overall branched evolution pattern suggests a persistent tumoral clonal expansion in the recurrent sectors (Fig. 2D–H). This also suggests that sector C is an outlier with a different clonal structure compared to the other recurrent sectors.

Regional Sectors Show Variable Treg and T-Helper Infiltrate and a Conserved PD-L1 Signature

To quantify the tumor infiltrating immune populations, we performed mIHC imaging and quantification. Histologically, the primary tumor exhibits the classical features of a GBM, including palisading necrosis, mitotic figures, and microvascular proliferations. The recurrent tumor was diagnosed as gliosarcoma histologically. Recurrent tumor A had mostly sarcomatous characteristics, with fascicles throughout and a necrotic zone along an edge. Recurrent B had a mixed, heterogeneous sarcomatous aspect and a glial appearance morphologically. Recurrent C appeared largely gliomatous with extensive necrotic areas. We observed intense, homogenous CD68+ macrophage infiltrate in all samples. Overall, we observed low adaptive immune infiltrates across all samples. In the primary tumor and recurrent sector C, CD3+ T-cell infiltration was scarce and mostly observed in pockets of stromal areas, whereas it was rarely seen in between individual tumor cells. Recurrent tumor sector A had strong CD20+ B- and CD3+ T-cell infiltration, including lymphoid follicle formations in some areas of the sample, whereas recurrent tumor B had relatively homogenous CD3+ T-cell infiltrate and sparse, rare individual CD20+ B cells in between tumor cells. In the primary tumor and recurrent sector B, membranous and cytoplasmic expression of PD-L1 was observed in individual tumor cells. Recurrent tumor sectors A and C had intense PD-L1 staining surrounding necrotic regions, which was mostly interpreted as nonspecific staining (Fig. 3A; Supplemental Fig. 2). We were able to quantify FOXP3+ regulatory T-cells and DC-LAMP+ mature dendritic cells in recurrent sectors A and C and found higher density of both markers in A compared to C (see Supplemental Table 2).

Figure 3.

(A) Multiplex immunohistochemistry (IHC)-stained images for the primary tumor, and recurrent sectors A–C (top to bottom) placed with the respective per-marker quantification. (B) Estimated percentage of immune cell types based on bulk RNA-seq xCELL deconvolution. (C) Proportion of TMB attributed to neoantigens. The total number of predicted immunogenic (ic50 < 500) neoantigens for each sample is labeled on top of the barplots. (D) Selected immunogenic neoantigens for loci that have undergone immune editing relative to the primary tumor (top) and relative to the recurrences (bottom). Text inside each cell indicates the minimum ic50 for predicted neoantigens in that gene.

We additionally estimated immune cell fractions using xCELL transcriptomic deconvolution (Aran et al. 2017). Although such methods are known to have strong limitations, especially compared to the mIHC assay, they can help offer some insight into how well RNA-seq alone estimates tumor immune quorum. Overall, transcriptomic deconvolution estimated similarly intense monocyte infiltrate levels across samples and relatively high macrophage and CD4+ T-memory cell infiltrate. There was low naive B-cell infiltrate and almost no B-memory cell infiltrate estimated across samples. xCell detected variable levels of Th1 and Th2 T-helper cell infiltrate in which A and B have the highest estimated cell fractions and primary and recurrent C have low to no estimated T-helper infiltrate. Low, but somewhat variable Treg infiltrate was detected in the samples using xCell (Fig. 3B).

Neoantigen Profiling Reveals Immune Editing

To estimate the extent of immune editing between recurrent sectors and the primary lesion, we sampled the neoantigen repertoire via in silico predictions of HLA binding affinity to neopeptides arising from somatic mutations (see Methods). Although suboptimal compared to mass spectroscopy, such predictions have proven to yield useful priors for identifying the subset of immunogenic mutations in a scalable fashion (Desrichard et al. 2016; Nogueira et al. 2018; Shembrey et al. 2019). The total immunogenic neoantigen burden was highest in recurrent sectors A and B and lowest in recurrent C. The immunogenic neoantigen burden in primary was only slightly higher than that of recurrent C (Fig. 3C). Indeed, we found that 110 somatic neoantigens were not likely immune-edited, 96 were later gained in the recurrent sectors, and only 14 were originally in the primary tumor but lost in recurrence (Supplemental Table 1). We note that the conserved PTEN mutation results in a potentially immunogenic neoantigen, as well as the conserved EGFR mutation and the novel, more immunogenic EGFR mutation that is specific to recurrent sector C (Fig. 3D).

Transcriptomic Heterogeneity of Regional Sectors Contains a Strong Immune Component

We examined differential gene expression among the regional sectors using a pooled-comparison strategy to help mitigate sample size bias. In the pooled comparison of recurrent sector C versus recurrent sectors A, B, and the primary tumor, we found evidence suggesting that sector C was enriched in aberrant expression of the H19 locus. We also found a global up-regulation in the pooled comparison of recurrent B versus A, C, and the primary tumor, including evidence of aberrant expression in the MAPK4 gene and several HOX genes (Supplemental Fig. 3). In cases in which the comparison involved a single sample in a group, variance estimates were estimated from all samples.

Using gene set enrichment analysis (GSEA) from KEGG, MSigDB Hallmark, Biocarta, and Reactome pathway signatures, we computed we found evidence for strong regionally relative immune pathway dysregulation among the sectors. For example, we find significant up-regulation of T-helper cell pathways, T-regulatory pathways, and the antigen processing and presentation pathway specific to primary and recurrent C. Further analysis revealed a MAPK signaling signature specific to recurrent sector C (P = 0.0386, normalized enrichment score = −1.40) and evidence toward global KRAS pathway dysregulation in recurrent sector B and an epithelial to mesenchymal transition (EMT) pathway down-regulation in sectors A and B. We also find a CTLA4 pathway up-regulation specific to the primary tumor and a PD-1 pathway down-regulation in sectors A and B (Fig. 4A; Supplemental Table 3). GSEA on the anti-PD-1-response gene sets from the Zhao et al. cohort revealed a pretreatment responder signature in sector C (Supplemental Fig. 4A).

Figure 4.

(A) Heatmap showing normalized enrichment scores for a curated set of pathways that is significantly enriched among regionally relative differentially expressed genes (DEGs). (B) Principal component analysis of patient sector gene expression profiles overlaid onto the profiles of 155 TCGA-GBM samples. Color indicates the log(Months Overall Survival) and shape indicates survival status at last follow-up. (C) Voom-normalized and batch-corrected expression of checkpoint gene expression across samples after integrating TCGA-GBM expression profiles. P-values describe unusual variance among the sectors compared to TCGA (see Methods).

To contextualize the regional expression profiles in terms of a “typical GBM,” we projected them onto a set of 155 TCGA-GBM samples using principal component analysis. This revealed that patient sectors achieve significant extremes of variation, implying a large relative scale of regulatory and cell-type heterogeneity at the level of gene expression (Fig. 4B). Upon integrating expression profiles for a defined set of immune marker genes, we found that sectors A and B overexpressed FOXP3, whereas primary and sector C showed considerably lower FOXP3 expression. Indeed, we find that FOXP3 expression varies more across regional sectors than expected given the TCGA cohort (P = 0.0118, via resampling). PD-L1, on the other hand, was shown to be significantly less varied among the regional sectors than expected in the TCGA cohort (P = 0.0409, via resampling) (Fig. 4C; Supplemental Fig. 4B).

Taken together these observations imply that primary and recurrent C cluster together at the transcriptomic level and not just the mutational level. Additionally, this provides evidence at the pathway level for a differential T-cell activation signature.

DISCUSSION

We characterized the intratumoral and immune heterogeneity of a primary GBM and multiple recurrent tumor sectors to provide insight into the GBM molecular landscape in the context of an unusually long survival following adjuvant anti-PD-1 treatment. Although recent studies have observed correlates of checkpoint blockade response in nonhypermutated GBMs, they typically rely on single-biopsy data and do not directly examine the characteristics of intratumoral heterogeneity in the context of anti-PD-1 responders. Here, we sought to extensively characterize regional heterogeneity in this recurrent GBM to better understand the molecular mechanisms behind the extended survival and evaluate whether the enrichment in response signature was uniform across recurrent sectors. Interestingly, we found an enrichment in both responder and nonresponder signatures across recurrent sectors.

This case study is representative of a typical GBM, including an unmethylated MGMT locus, a wild-type IDH genotype, and a low somatic mutation burden (Hegi et al. 2009; Delgado-López and Corrales-García 2016; Nam and de Groot 2017). Despite these classic markers suggesting a typical poor prognosis, the patient's overall survival was approximately double what would be expected. Moreover, the patient surpassed the median survival of GBM patients with MGMT promoter methylated tumors by 7 mo (Hegi et al. 2009). Historical data on GBM patients indicates that most patients, especially those with our patient's tumor profile, do not generally survive >6–8 mo after recurrence. Further, based on TCGA integration, the probability of our patient's long survival occurring by chance alone was relatively low. Finally, volumetric MRI analysis shows that two of the three recurrent tumors did not grow back following anti-PD-1 therapy. These observations suggest that checkpoint blockade conferred a clinical benefit in this case (Thomas et al. 2012; Reardon et al. 2014; Delgado-López and Corrales-García 2016; Nam and de Groot 2017).

As noted, our patient has a mixed predicted response profile, with expressed markers found in both anti-PD-1 responders and nonresponders (Zhao et al. 2019). For example, the nonresponse PTEN alteration was found in all recurrent samples and the primary, whereas a TP53 mutation was found in recurrent tumors A and B, as well as a focal deletion in both TP53 and PTEN in the primary tumor. However, we also observed response markers, including the overall branched clonal evolution pattern, an observed silencing of neoantigens previously found in the primary tumor, BRAF and PIK3A copy-number alterations, and MAPK pathway alterations. Further, although sector C had a responder signature according to GSEA, this conflicts with volumetric evidence showing that C recurred, whereas A and B did not. This explicitly demonstrates the scale of intratumoral heterogeneity in this patient, as these biomarkers would have only been partially informative in predicting treatment response.

Projection onto the TCGA-GBM cohort revealed large-scale expression heterogeneity across sectors. In an immune-specific context, for given checkpoint-related genes, the expression variance was much larger than expected by chance. We hypothesize that this high variation within the primary tumor is reflected in the recurrences as evidence of the disparity of response markers to checkpoint blockade. Additionally, radiation therapy, which was part of the standard-of-care treatment, has been shown to elicit antitumor immune response and tumor clonal expansion (Rajani et al. 2019). Therefore, it may be possible that radiation therapy played a role in augmenting an antitumor immune response in some recurrent sectors and enhancing resistance in others. This hypothesis is corroborated by the observed immune editing denoted by a pronounced shift in gain of immunogenic neoantigens in the recurrence sectors. Further, the observed expansion of tumor subclones from clonal inference and tree reconstruction may also support this hypothesis. We additionally observed a marked increase in CD3+ T cells through IHC and T-cell populations via xCELL deconvolution in recurrent tumors A and B relative to the primary. Finally, this hypothesis is further substantiated by evidence of an alteration in the interferon-gamma pathway, which is a known mechanism of radiation-induced immune cell recruitment (Rajani et al. 2019).

One notable feature was a T-helper enrichment in recurrent sectors A and B both at the pathway level and via cell-type proportion estimates derived from xCELL (Aran et al. 2017). Similarly, the immune landscape observed in the IHC underscores differences found across recurrences with and highlights that CD3+ T cells are robustly distributed across sectors A and B, whereas they are scarce in the primary tumor and sector C. However, there are also observable patterns among them that may help explain treatment response. For example, the overall pattern of clonal evolution, with recurrent sectors A and B clustering together and similarly primary and recurrent C, was observed in both the immune and genomic surveys. Of all the recurrent tumors, sector A also had the highest PD-L1 expression according to IHC quantification. Additionally, sector A had the highest mutation burden and the most immunogenic neoantigen burden. These observations, which are all associated with an increased probability of anti-PD-1 response, suggest that recurrent sector A may have contributed the most to eliciting PD-1 blockade response, whereas sector C may have contributed the least. This hypothesis is additionally corroborated in the longitudinal MRI imaging that followed the resection of the recurrent tumor, in which the slow progression began near the medial (C-like) region, grew toward the inferior (B-like) region, but did not reach or arise from the lateral (A-like) region. Although we lack post-treatment molecular data to confirm this hypothesis, future studies may be able to further evaluate heterogeneity-driven mixed response profiles.

In summary, we provide a detailed, comprehensive, multiregional tumor profiling of a typical GBM that successfully responded to immune checkpoint blockade. Overall, we show key differences between the recurrent sectors that could have affected treatment outcomes and highlight the need for further studies that can provide a broader survey of intratumoral heterogeneity across anti-PD-1-responsive GBMs. These proposed studies may instruct us how to take advantage of tumor heterogeneity at the immunologic level in this patient population to successfully tailor treatment to GBM patients.

METHODS

Patient Recruitment and Sample Collection

All tumor samples used in this study were derived with consent from a patient enrolled in our institutional IRB-approved cancer biorepository tissue procurement protocol. Resected tumor tissue was processed for clinical diagnosis according to clinical standard operating procedures. A diagnosis of GBM was established according to the criteria of the WHO Classification of Tumors of the Central Nervous System.

Genomic Characterization

Whole-Exome Sequencing

Genomic DNA extraction was conducted using Purelink Genomic DNA Mini kits (Thermo Fisher Scientific) to yield ample material for NGS. Whole-exome sequencing was conducted using a NextSeq 550 (Illumina) sequencer using the SureSelect WES v5 chemistry (Agilent) and 125-bp paired-end, high-output sequencing mode. Samples recurrent A, recurrent B, recurrent C, primary tumor, and matched blood yielded 256×, 241×, 218×, 271×, and 111× coverage across 246 million, 235 million, 215 million, 261 million, and 106 million reads, respectively (Supplemental Table 4).

Data processing for WES samples was performed using a custom Nextflow pipeline that is available on GitHub (https://github.com/losiclab/exoseq). Raw reads were trimmed and aligned to the hg38 reference genome using Trim Galore! (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/), bwa mem, and SAMtools, respectively (Li et al. 2009; Li 2013). Duplicate reads were then marked using Picard MarkDuplicates (https://broadinstitute.github.io/picard/), and bam quality scores were recalibrated for known technical bias using GATK4 base recalibration (McKenna et al. 2010). Quality control metrics were compiled using FASTQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) for raw reads and Picard CollectMultiMetrics for aligned reads.

Somatic Mutation Calling

Somatic variants were called using Mutect2 in WES tumor-blood matched data, with the intervals parameter set to all coding regions, and with 1000 Genomes as the germline resource. Variants were filtered for quality control using GATK FilterMutectCalls with default parameters (McKenna et al. 2010; Cibulskis et al. 2013; Poplin et al. 2018). Only variants with a VAF > 5% were retained for further analysis. Variants were annotated with ANNOVAR (Wang et al. 2010).

Neoantigen Prediction

Patient HLA haplotypes were determined using Optitype. Neoantigen prediction was performed on somatic mutation callsets using the netMHCpan version 4.0. Putative neoantigens with a predicted ic50 >500 nM and more than five tumor reads supporting the mutation were retained for further analysis (Lundegaard et al. 2008; Szolek et al. 2014).

CNV Calling

Copy-number status of whole-exome data was inferred using cnvkit (Talevich et al. 2016).

Clonal Architecture Modeling and Tumor Purity Estimation

Purity estimates were determined using THetA2. To determine cellular fraction and tumor subclone composition, read counts were first extracted from all mutated genomic positions using bam-readcount (https://github.com/genome/bam-readcount). Then, the kernel density distribution of VAFs was visualized for each chromosome. VAFs from Chromosomes 1, 7, 9, 12, and 14 were removed because of unusually high variance in an effort to reduce noise. Clonality was estimated with the SciClone package using the filtered VAFs, CNV calls, and a minimum depth parameter of 300 reads (Miller et al. 2014). The phylogenetic tree of the identified clones was modeled and visualized using the ClonEvol package in R (Dang et al. 2017).

RNA Sequencing

Poly(A)-selected RNA sequencing was performed with 100-bp reads using an Illumina HiSeq 2500 instrument. Raw reads were aligned to the hg38 reference genome, and gene-level read counts were extracted using STAR (Dobin et al. 2013). Aligned reads were further processed with Picard (https://broadinstitute.github.io/picard/) MarkDuplicates, followed by GATK SplitNCigarReads, IndelRealigner, BaseRecalibrator, and AnalyzeCovariates, respectively (McKenna et al. 2010).

Differential Gene Expression and Pathway Analysis

Using gene-level read counts, regionally relative DEGs were computed using DESeq2 with a design matrix that uses pooled sample comparisons in an effort to compensate for lack of statistical power (Love et al. 2014). Gene set enrichment analysis was applied to all DESeq2 results using the curated KEGG, Hallmark, Reactome, and Biocarta signatures publicly available from MSigDB (Liberzon et al. 2011, 2015).

TCGA Expression Integration

TCGA-GBM raw counts matrices (n = 155) were downloaded from the GDC data portal (https://portal.gdc.cancer.gov/), combined with the case resection samples, and corrected for batch effect, age, and gender. This was followed by principal component analysis.

For each checkpoint gene, bootstrap resampling was used to select the gene expression values for four random TCGA samples out of the 155, for a total of 2,081,158 resamplings, or 10% of the total possible unique combinations of TCGA resamplings. For each of these resamplings, the variance was computed, followed by visualization as a kernel density estimate. The variance among the four case biopsies was overlaid on this distribution, and the P-value was calculated as the area under the curve from the patient sample variance to the nearest extreme of the distribution.

Infiltrating Lymphocyte Characterization

Transcriptomic Immune Cell Deconvolution

Transcriptomic profiles for the samples were subject to cell-type deconvolution using xCell (Aran et al. 2017), using the normalized CPM expression values as input for the algorithm. Results were visualized using R.

Multiplex Immunohistochemistry

Four-micrometer tissue sections were cut from the primary and recurrent tumor blocks. The slides were stained with CD3 (Ventana, 2GV6), CD20 (Dako, L26), PD-L1 (Cell Signaling, E1L3N), CD68 (Dako, KP1), FOXP3 (Abcam, 236A/E7), and DC-LAMP (Novus Biologicals, 1010E1.01) by using Multiplex Immunohistochemical Consecutive Staining of Single Slide (MICSSS) methodology (Remark et al. 2016). Whole slide images were generated using a slide scanner (Hamamatsu NanoZoomer S60) for each of the multiplex staining cycles. QuPath image analysis software (version 0.2.0 m2) was used for the quantification of positive cells for each marker (Bankhead et al. 2017). To analyze marker coexpression, a pseudofluorescence composite image of all chromogenic markers was created as follows for each case.

The same regions of interest (ROIs) were selected from images of each marker and exported as lossless PNG formatted images without any downsampling. Images of different immunostains belonging to the same ROI were then transferred to the Fiji-ImageJ 1.52h image analysis platform and coregistered by using the TrakEM2 plug-in (Cardona et al. 2012; Schindelin et al. 2012). Color deconvolution was done by using H-AEC vectors for each image. This process split the red–green–blue (RGB) type MICSSS images into three 8-bit channels including hematoxylin (blue), AEC (red chromogen color), and residual (green) channel. The best hematoxylin channel was selected to be used as the nuclear channel, and AEC channels for each marker were saved. AEC channels representing staining of each marker were assigned to different colors by using the lookup tables (LUTs) function of Fiji. The hematoxylin channel was assigned to blue color in order to mimic fluorescent DAPI staining, whereas other channels were assigned a color per marker. Color inversion was done on all channels and then merged to achieve a multiplexed pseudofluorescent image. Brightness and contrast settings were optimized to facilitate visualization for each immunostain channel by comparison with original chromogen images by a pathologist, although the underlying image pixel values were not changed for quantification by QuPath. Composite multiplexed images were analyzed for positive cell quantification of single markers and/or combination of markers by using QuPath image analysis software (Bankhead et al. 2017).

ADDITIONAL INFORMATION

Data Deposition and Access

Raw sequencing data is deposited in Sequence Read Archive (SRA) (https://www.ncbi.nlm.nih.gov/sra/) under accession number PRJNA566392. The analysis source code is deposited in GitHub (https://github.com/LosicLab/gbm-case-study).

Ethics Statement

All tumor samples used in this study were derived with consent from a patient enrolled in the Icahn School of Medicine at Mount Sinai institutional IRB-approved (GCO#1: 18-2163) cancer biorepository tissue procurement protocol.

Author Contributions

A.H. and R.Y. treated the patient and performed the surgical biopsies, respectively. R.S. performed the library preparation and subsequent sequencing assays with assistance from N.T., I.L., and R.Y. S.G., G.A., and I.L. performed the mIHC assay and analysis. K.N. performed the MRI volumetric analysis. P.R. performed the raw large data and downstream bioinformatics analysis for the study with assistance from B.L. P.R. prepared the manuscript with assistance from B.L., A.H., S.G., R.S., G.A., and R.Y. B.L., A.H., and R.Y. designed the project and supervised the work.

Funding

This work was funded by the Genetics and Genomics Department, Icahn School of Medicine at Mount Sinai, as well as the Icahn Institute for Data Science and Genomic Technology.

Competing Interest Statement

S.G. is a named coinventor on an issued patent for multiplex immunohistochemistry to characterize tumors and treatment responses. The technology is filed through Icahn School of Medicine at Mount Sinai (ISMMS) and nonexclusively licensed to Caprion (provides integrated platforms in proteomics and immune monitoring services to the pharmaceutical and biotechnology industry). Mount Sinai has received payments associated with licensing this technology, and both Mount Sinai and S.G are entitled to future payments. S.G. reports receiving commercial research funding from Immune Design, Pfizer, Janssen R&D, Genentech, Regeneron, and Bristol-Myers Squibb and is a consultant/advisory board member for Merck, OncoMed, Neon Therapeutics, and B4CC.