Integrated proteogenomic characterization of glioblastoma evolution

Abstract

The evolutionary trajectory of glioblastoma (GBM) is a multifaceted biological process that extends beyond genetic alterations alone. Here, we perform an integrative proteogenomic analysis of 123 longitudinal glioblastoma pairs and identify a highly proliferative cellular state at diagnosis and replacement by activation of neuronal transition and synaptogenic pathways in recurrent tumors. Proteomic and phosphoproteomic analyses reveal that the molecular transition to neuronal state at recurrence is marked by post-translational activation of the wingless-related integration site (WNT)/ planar cell polarity (PCP) signaling pathway and BRAF protein kinase. Consistently, multi-omic analysis of patient-derived xenograft (PDX) models mirror similar patterns of evolutionary trajectory. Inhibition of B-raf proto-oncogene (BRAF) kinase impairs both neuronal transition and migration capability of recurrent tumor cells, phenotypic hallmarks of post-therapy progression. Combinatorial treatment of temozolomide (TMZ) with BRAF inhibitor, vemurafenib, significantly extends the survival of PDX models. This study provides comprehensive insights into the biological mechanisms of glioblastoma evolution and treatment resistance, highlighting promising therapeutic strategies for clinical intervention.

Keywords: BRAF; longitudinal glioblastoma; neuronal; proteogenomics; recurrence; synapse.

Figure 1.. The molecular trajectory of longitudinal glioblastoma

(A) Summary of profiling platforms and number of matched longitudinal GBM patient samples analyzed by each molecular platform. The color white indicates matched samples not profiled by the indicated platform. SNU: Seoul National University, PS: Hôpital de la Pitié-Salpêtrière, SMC: Samsung Medical Center, CNU: Chonnam National University, and NCC: National Cancer Center. (B) Comparison of tumor mutational burden and aneuploidy fraction between IDH mutant and IDH wild-type primary and matched recurrent GBMs. Box plots span from the first to third quartiles, middle line represents median, and whiskers represent the 1.5x interquartile range. The P-values were calculated by the Wilcoxon rank-sum test. (C) Somatic genomic landscape of GBM driver genes grouped by oncogenic pathways. IDH mutant samples are grouped on the right side of the panel. Patients were ordered according to progression-free survival (PFS, months). Genomic alterations shared between primary and recurrent tumors are in green, private alterations are in yellow (primary specific) or purple (recurrent specific); in red are shared alterations with mutation replacement). The bar graph on the right represents the overall frequency rate of each shared or private genomic alteration. (D) Bubble plot illustrating the frequency of fCNVs and non-synonymous mutations of GBM driver genes in exclusively primary (yellow, left axis), exclusively recurrent (purple, right axis), and shared (green, upper axis) between primary and recurrent GBM. GBM driver genes with genetic alterations in at least 4 primary recurrent pairs are represented. Red and blue gene labels denote oncogenes and tumor suppressor genes, respectively. The size of each node represents the number of patients harboring the corresponding genetic alterations. (E) Bubble plot illustrating the frequency of fCNVs and non-synonymous mutations of all the genes in exclusively primary (yellow, left axis), exclusively recurrent (purple, right axis), and shared between primary and recurrent GBM (green, upper axis) GBM tumors (Fisher’s exact test, p < 0.10). Representative genes involved in neuronal differentiation, cell cycle, mitosis, and DNA repair are indicated. See also Figure S1 and Table S1.

Figure 2.. Impact of EGFR genomic alteration on downstream signaling upon recurrence

(A) Analysis of up-regulated protein in key oncogenic GBM pathways including RTK-RAS, p53, cell cycle, and chromatin modification. The pie charts show the percentage of primary (yellow) and recurrent (purple) samples showing up-regulated proteins. Selected phosphosites are indicated and color-coded according to elevation in primary or recurrent tumors. (B) Somatic genomic landscape of EGFR alterations, including amplification (Amp), mutation (Mut), EGFRvIII, and Fusion (Fus). Patients have been categorized as “Shared” if EGFR alterations are shared between primary and recurrent tumors (green), “Loss” (blue) if EGFR alterations are lost at recurrence, or “Gain” (Red) if it was gained in the recurrent tumor. (C) Comparison of EGFR protein abundance between matched primary (P) and recurrent (R) tumors in the corresponding patient groups. Box plots span from the first to third quartiles and middle line represents median. The P-values were calculated by the Wilcoxon rank-sum test. (D) Violin plots of the phosphorylation abundance in EGFR sites between primary and recurrent GBM patients in the “Shared” group. The P-values were calculated by the Wilcoxon rank-sum test. (E) Correlation plot of EGFR protein abundance with key regulators of PI3K-AKT-PTEN-mTOR and the RAS-RAF-MEK-ERK downstream pathways. Pearson Correlation Coefficient (PCC) is indicated once each protein abundance was normalized across samples. (F) Dots plot showing the changes of phosphorylation abundance in EGFR substrates when comparing recurrent and primary GBMs. The P-values were calculated by the Wilcoxon rank-sum test. See also Figure S2.

Figure 3.. Proteogenomic characterization of longitudinal GBM reveals increased neuronal activities in recurrent GBM

(A) Volcano plot representation of proteins (upper panel) and phosphoproteins (bottom panel) with differential abundance. x-axis indicates the difference of median expression between primary and recurrent GBMs; y-axis indicates the statistical significance derived from Wilcoxon rank-sum test. The P-values for each protein and phospho-protein were calculated by the Wilcoxon rank-sum test. (B) single sample Gene Set Enrichment Analysis (ssGSEA) analysis of proteins and phosphoproteins between primary and recurrent GBMs using all genesets available in the MSigDB. The P-values were calculated by the Wilcoxon rank-sum test for each geneset between primary and recurrent tumors. (C) Comparison of major synaptic pathway activities between matched primary and recurrent tumors. The lines within each violin plot represent the 25^th, 50^th, and 75^th quantiles. The P-value was calculated by the Wilcoxon rank-sum test. (D) Scatter plot showing the ratio of mRNAs, proteins, and phosphoproteins expression/abundance significantly changed between recurrent and primary GBMs. Red and blue dotted-line quadrants indicate upregulated and downregulated proteins, and phosphoproteins, respectively without any mRNA changes. The P-values were calculated by the Wilcoxon rank-sum test. (E) GO enrichment analysis of proteins and phosphoproteins from the red and blue quadrants in (D). The q values were calculated using the Benjamini-Hochberge method. (F-G) Cis and trans effects of major genomic alterations including mutations and CNVs on protein (F) and phosphoprotein (G) levels. Cis and trans effects of each genomic alteration were categorized into primary-specific (yellow), recurrent-specific (purple), or shared (green); see the track on the left of the panel. The P-values were calculated by the Wilcoxon rank-sum test. See also Figure S3 and Table S2.

Figure 4.. Integrative multi-omics clustering reveals the transition of GBM at diagnosis toward neuronal state at recurrence

(A) Heatmap of multi-omic features significantly associated with each subtype according to the proteomic-based functional classification of primary and recurrent GBM. Upper track, functional classification; bottom track, tumor type. Columns are individual tumors and rows are features. First panel, mutation calls are in green (p < 0.10, Fisher’s exact test). Second panel, fCNV gain or amplification calls are in red and orange (p < 0.05, Fisher’s exact test), respectively. Third panel, fCNV heterozygous or homozygous deletion calls (p < 0.05, Fisher’s exact test) are in blue and cyan, respectively. Fourth and fifth panels, 150 highest scoring genes/proteins in the gene/protein expression/abundance ranked list of each of the four GBM subtypes (MWW-GST test). Sixth panel, significant outlier phosphorylated proteins in each functional GBM subtype (p < 0.05, BlackSheep). Biological pathways significantly enriched (Fisher exact test, p < 0.05) and representative genetic alterations specific to each GBM subtype are indicated on the right. (B) Proteomic-based functional subtyping of primary and recurrent GBM. The transition plot of functional subtypes in primary and recurrent GBM shows an increased frequency of the NEU subtype at recurrence (p = 0.0099, χ2 test). (C) Bubble plot frequency of primary non-neuronal to recurrent neuronal matched pairs of GBM harboring fCNV in genes significantly associated with primary or recurrent tumors (Fisher’s exact test, p < 0.05). (D) Comparison of MGMT, MSH2, MLH1, and PMS2 protein abundance between matched primary and recurrent tumors that transition into C2-type. Box plots span from the first to third quartiles and middle line represents median. The P-values were calculated by the Wilcoxon rank-sum test. (E) Graphical illustration of the down-regulation of mismatch repair (MMR) encoding proteins leading to TMZ resistance. (PPR : proliferative/progenitor, NEU : neuronal, GPM : glycolytic/plurimetabolic, MTC : mitochondrial, CL : classical, MES : mesenchymal, PN : proneural) See also Figures S3-5.

Figure 5.. Single-cell analysis reveals increase of tumor-intrinsic neuronal state and neurons and oligodendrocytes in the TME of recurrent GBM

(A) UMAP plot of single cells from matched primary and recurrent tumors (4 patients, 8 tumors), colored according to malignant or non-malignant cell subpopulations (left panel) and primary or recurrent status (right panel). (B) UMAP plot of single cells colored according to the median of the copy number of genes in chromosome 7 (left panel) and chromosome 10 (right panel) predicted by inferCNV as hallmarks of glioblastoma. (C) UMAP plot of single cells colored according to the malignant and non-malignant cell types (left panel) and averaged gene expression of corresponding marker genes (right panel). OPC, Oligodendrocyte progenitor cells. (D) Bar plots showing the cell type composition in each tumor analyzed. (E) Enrichment of non-malignant cell types between primary and recurrent GBM. Asterisks indicate chi-squared test derived standardized residuals above 2.5 as index of the statistical significance and strength of the association. (F) Tumor cell type distribution based on functional pathway (upper panel) and glioma cell state (lower panel). Asterisks indicate chi-squared test derived standardized residuals above 2.5 as index of the statistical significance and strength of the association. AC, Astrocyte. (G) Box plots showing neuronal and proliferative-progenitor state enrichment scores (NES) in each individual tumor cell from primary and recurrent GBM. Cells are colored according to the neuronal (blue) or proliferative-progenitor (cyan) state enrichment score. Box plots span from the first to third quartiles and whiskers show 1.5× interquartile range (Wilcoxon test). (H) Two-dimensional plot of functional tumor cell state enrichment scores. Each quadrant corresponds to one GBM subtype, and the position of dots (tumor cells) reflects the relative subtype-specific NES of each cell as indicated on the x and y axes; cells are colored according to primary and recurrent status.

Figure 6.. Functional validation of enhanced synaptic activity in recurrent GBM

(A) Co-culture of paired patient-derived GBM cells (primary and recurrent) and mouse cortical neuron cells. Expression of MAP2, synapsin-1 in paired GBM cells was measured by immunoblot analysis. The number of presynapses was analyzed. Data are represented as mean ± SD from triplicate wells. Statistical significance was assessed using Student’s t-test.Scale bars = 20 μm. (B) Immunofluorescence stains with MAP2 (red), synapsin-1 (green), and human-specific nuclei (white) antibodies in the co-culture condition of (A). The neurite length of GBM cells and the number of the synapse was analyzed by counting colocalized punta of synapsin-1 and MAP2. Data are represented as mean ± SD from triplicate wells. Statistical significance was assessed using Student’s t-test. Scale bars = 20 μm. (C-D) Representative images (C) and quantification (D) of immunostaining with DAAM1 (red) and SOX2 (green) antibodies in patient GBM tissue specimens. The ratio of DAAM1+ cells among SOX2+ cells is shown in primary-recurrent paired GBM tissues whole sections from 4 patients. Scale bars = 40 μm. See also Figure S6.

Figure 7.. In vivo modeling of neuronal transition at recurrence

(A) Overview of experimental design for proteogenomics analysis and modeling for therapeutic resistance by TMZ treatment in PDX models. (B) Survival curves of mice treated with 10 mg/kg TMZ (n=4) or vehicle (n=4) (p = 0.0288, log-rank Mantel-Cox test). T indicates the time point (in days) of tissue preparation for proteogenomic analysis. (C) Heat map of the transcriptome (upper panel), proteome (middle panel), and phospho-proteome (bottom panel) of PDX-GBM treated multiple times with TMZ. Differential expression/abundance at each time point was compared to the other times (log2(FC) TMZ/control > 0.58). Biological pathways significantly enriched for each platform analyzed are presented on the left (Fisher exact test, p < 0.05). (D) Co-regulation clustering of temporal proteomics data from PDX contrasting TMZ and vehicle treatment reveals two distinct clusters of proteins having the opposite expression pattern by day 35 after administration. (E) Pathway enrichment analysis of the two co-expressed clusters of proteins from Figure 7C. (F) Representative multiplex immunohistochemistry for MAP2 and DAAM1in TMZ- and vehicle-treated PDX model. Scale bars = 20 μm. See also Figure S7.

Figure 8.. Identification of BRAF as a key therapeutic target of GBM evolution

(A) Heat map of master kinases with differential activity in each unsupervised multi-omics-based (left) and functional (right) GBM subtypes. (B) Box plots showing BRAF activity in each individual tumor stratified according unsupervised multi-omics based (left) or functional GBM subtype (rigth panel). Tumors are colored according to primary or recurrent status. Box plots span from the first to third quartiles, middle line represents median, and whiskers show 1.5× interquartile range. (left panel, C2 versus all other subtypes: p = 1.71e-08; right panel, NEU versus all other subtypes: p = 2.21e-05; Wilcoxon test) (C) Representative images (left) and quantification (right) of spreading assays in paired GBM cells treated with vehicle, 100 nM vemurafenib, or 100 nM dabrafenib (PMZ). Scale bars = 10 μm. Data are represented as mean ± SD from triplicate wells. Statistical significance was assessed using Student’s t-test (**p < 0.01; ***p < 0.001). (D) Immunoblot analysis of MAPK signaling-related proteins in paired GBM cells treated with vemurafenib. (E) Survival analysis of in vivo PDX models that were treated with either vehicle (n=6), TMZ (n=6), vemurafenib (n=6), or the combination of TMZ with vemurafenib (n=6) (p = 7.2X10⁻³, log-rank Mantel-Cox test). (F) Representative MRI and H&E sections of the mouse brains from (E). Scale bars = 2 mm. See also Figures S7 and S8.