Proteomics, Post-translational Modifications, and Integrative Analyses Reveal Molecular Heterogeneity within Medulloblastoma Subgroups

Abstract

There is a pressing need to identify therapeutic targets in tumors with low mutation rates such as the malignant pediatric brain tumor medulloblastoma. To address this challenge, we quantitatively profiled global proteomes and phospho-proteomes of 45 medulloblastoma samples. Integrated analyses revealed that tumors with similar RNA expression vary extensively at the post-transcriptional and post-translational levels. We identified distinct pathways associated with two subsets of SHH tumors, and found post-translational modifications of MYC that are associated with poor outcomes in group 3 tumors. We found kinases associated with subtypes and showed that inhibiting PRKDC sensitizes MYC-driven cells to radiation. Our study shows that proteomics enables a more comprehensive, functional readout, providing a foundation for future therapeutic strategies.

Keywords: MYC; NU-7441; SHH; mass spectrometry; medulloblastoma; multi-omics; network integration; phospho-proteomics; proteo-genomics; radio sensitization.

Figure 1:. Summary of data types included in this study, depth of proteomic data types, and cohort composition.

The extent of the data from proteomics, including post-translational modifications, is summarized at the top. pSTY: phosphorylation on serine, threonine, or tyrosine detected after immobilized metal affinity chromatography; pY: phosphorylated tyrosine detected after antibody purification; Total: the total number of features identified; All Samples: the number of features measured across all samples, i.e. without any missing values. The number of samples covered by each data type and their split by subgroups are summarized at the bottom. G3: Group 3; G4: Group 4; WGS: whole genome sequencing. Proteomics includes proteome, pSTY, pY, and acetylomics data sets. See also Table S1.

Figure 2:. Comparison of clustering results.

(A) The optimal clustering of DNA methylation data, RNAseq, and proteome, as determined using Pearson correlation as distance metric. k, number of clusters. Consensus scores are indicated using a color scale from white (samples never cluster together) to blue (samples always cluster together). (B) Comparison of the assignment of samples using the four “consensus subgroups” (Taylor et al., 2012), the DNA Methylation-based subtype-calls assigned in Northcott et al. (2017), which included most samples used in our study, and RNA expression assignments based on application of the classifier described in Cho et al. (2011). NA: no assignment available. See also Figures S1-S3 and Tables S1 and S2.

Figure 3:. Molecular and functional differences between SHHa and SHHb.

(A) Hierarchical clustering of 510 proteomic features that differ significantly between SHHa and SHHb (FDR < 0.005; ANOVA). Summaries of gene sets that significantly differ (FDR < 0.01; hypergeometric overlap tests) between SHHa and SHHb are shown alongside the heatmap. (B) A schematic of the pathway for the glutamatergic synapse. Each box represents a gene or gene set taken directly from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. The left and right boxes summarize SHHa and SHHb protein levels, respectively. Grey boxes indicate a lack of proteomic data for the gene product. Small white circles represent metabolites, solid directed arrows activation, and lines ending with a bar inhibition. Dotted lines indicate transport of metabolites. Solid lines and connected boxes represent physical interactions in the cell. Figure adapted from KEGG pathway hsa04724 (Kanehisa et al., 2017). (C) Clinical annotations of and recurrent genomic alterations in SHHa and SHHb samples. (D) Box-and-whisker plots of Spearman correlations between the mRNA and protein levels for approximately 8,700 genes for each tumor sample. Boxes range from first to third quartile in each group, with line indicating the group-median. The whiskers extend to the lowest and highest data points within 1.5 times the interquartile range from the box (p value shown from Mann Whitney U test). See also Figure S4 and Tables S3 and S4.

Figure 4:. G3a medulloblastomas have post-translational modifications of MYC protein that strongly correlate with MYC downstream activity and patient outcome.

(A) Heatmap of multi-omic data for MYC, showing protein levels, post-translationally modified peptides, mRNA expression, DNA copy number, and DNA methylation. Scale bar is of log2 values for all data sets except DNA methylation, which ranges from 0 to 1 on a linear scale. Methylation values, converted to percentages, are superimposed on the heatmap. The phosphorylation locations are provided both based on RefSeq (as used in our data sets) and Uniprot (frequently used in the literature). S77 or 79: phosphorylated at either S77 or S79 but not at T73; S359* / T365: phosphorylated at T365 as well as a proximal site that could not be uniquely determined. (B) Representative confocal microscopy images showing pS62 and pT58 MYC in G3a, G3b, and Group 4 tumors. (C) Scatter bar charts showing mean immunofluorescence densities quantified across 50–100 nuclei in 3–4 fields of view captured at 40x for MYC pS62 and pT58. Error bars represent standard error of means (SEM) and significance reported is from two-tailed unpaired t-tests. (D) Kaplan-Meier plots of overall survival and progression-free survival of Group 3 medulloblastoma from the current cohort (left) and the current cohort combined with Group 3 samples from Cho et al. (2011). The survival analysis on the combined cohort was performed using the classifier defined by the proteomic clusters G3a and G3b (“Proteome-based Classifier” in the center) and the original classifier from Cho et al. (2011) (c1/c5, right, gray scale). In the combined cohort with proteomebased group assignments (G3a/G3b), both outcome predictions reached statistical significance (p values as shown; Log-rank test). See also Figure S5.

Figure 5:. Medulloblastoma subtypes differ in kinase regulation and substrate levels.

Two different methods used to analyze the kinome are shown. On the left, upstream kinases are predicted from differential phosphopeptides between subgroups using PhosphoSitePlus database (Hornbeck et al., 2015). On the right, kinases are predicted from scoring differential phosphopeptides using sequence specificity motifs from Scansite (Obenauer et al., 2003). Heatmaps show the median levels of peptides matched to an upstream kinase, with the number of peptides matching each kinase shown in parentheses. Kinases found by both methods are annotated with an asterisk. Bars on the side of the heatmaps indicate whether target peptides correlate with protein or phosphorylation levels of upstream kinases; and if DrugBank (Law et al., 2014) lists any drugs targeting the kinases. See also STAR Methods and Table S5.

Figure 6:. MYC status correlates with PRKDC phosphorylation, and predicts increased sensitivity to PRKDC inhibition with irradiation.

(A) Representative Western blots of medulloblastoma cell lines performed in at least triplicate. Antibody against pT58/S62 MYC detects either site or both sites. (B) Quantification of Western blots represented as fold change compared to expression in DAOY cells. All proteins are normalized to β-actin. Significance was determined via one-way ANOVA with a Dunnett multiple comparison test to compare the normalized signal for each antibody across the five cells lines; **p < 0.001, ***p < 0.0001. Error bars represent mean normalized signal ± SEM. (C) Correlations between normalized means for specified antibodies determined by calculating a Pearson correlation coefficient. Error bars indicate ± SEM, and are depicted but are not visible for some data points because of scale. (D) Representative confocal images of indicated cell lines showing pS62 MYC and pS2056 PRKDC. (E) Experimental design of dose-response curve of PRKDC inhibitor NU7441 for 18 hours prior to irradiation. (F) Viability assay of medulloblastoma cell lines treated as indicated. Plotted is the mean of 6 biological replicates ± SEM for each dose; note, small error bars for D458 are depicted but obfuscated by trend lines. (G) Histogram of mean IC₅₀ values ± SEM for DAOY vs D458 treated with NU7441 ± irradiation. *p < 0.01; ***p < 0.0001; ns = not significant; IR = irradiation.

Figure 7:. Network methods relying on known protein-protein interactions identify pathways relating to SHHb and G3a tumors.

Omics Integrator output showing network views of proteins, posttranslational modifications, and genomic alterations associated with G3a (A) and SHHb (B). Node shapes indicate data type and colors indicate log2-based fold change between groups as described in the legends. Phosphopeptides are labeled with their phosphorylation sites (based on RefSeq) after the ‘@’ symbol. Nodes associated with selected pathways are highlighted with yellow background. Grey nodes were added by Omics Integrator and have no associated proteomic data for our samples. (A) mRNA levels of transcriptional targets of MYC are shown at the bottom. Thick borders highlight proteins that are also shown as direct transcriptional targets of MYC. SNVs, indicated by diamonds, are color-coded to show which subtype the genomic alteration was seen in: G3a, red; G3b, blue. (B) Kinases that were also found by our independent Kinome Analysis (Figure 5) are highlighted with thick borders. The color of genomic lesions (diamonds) indicates the subtype in which they occur: SHHb, red; SHHa, blue. The location of color in a diamond indicates the type of genomic lesion: upper triangle, SNV; lower triangle, indel.