Acquisition and stratification of short- and long-term GBM survivors
To establish a cohort of primary GBM samples for downstream interrogation, GBM tumour samples (N = 128) were procured from the GLIOTRAIN (GT) biobank [10]. An additional N = 5 samples from the LIH biobank were subsequently added to the GT cohort to form the expanded GT-cohort (N = 133 patients samples total). Within this expanded GT-cohort, patients were stratified based on OS, identifying N = 18 STS, N = 82 Intermediate Term Survivors (ITS), and N = 33 LTS (Figure S1). As a first step, we assessed the distribution of PN, CL and Mes gene expression subtypes [20] within the expanded GT cohort. Classification into molecular subtypes showed 10.7% PN, 42.8% CL and 38.9% Mes tumours (Figure S2). Further analysis of OS based on molecular subtypes showed no significant differences (p = 0.4) (Figure S3).
(Phospho-)protein Analysis Of Gbm Samples Using Rppa Demonstrates Heterogeneity In Key Signaling Pathways In Gbm
RPPA analysis quantified signaling proteins in 126 samples of the expanded GT cohort [21], which first underwent patient-to-patient clustering to identify potential (phospho)proteomic subtypes. 4 distinct clusters were identified: cluster 1–4, consisting of 30, 58, 21, and 17 samples (Fig. 1; 1A and 1B). Cluster 1 had high expression of apoptotic signaling proteins of apoptotic signaling proteins BCLXL, SMAC/DIABLO (Fig. 1C, 1D) and BAX (Figure S4), as well as PARP (Fig. 1E), PDK1 and FAK (Figure S4). Cluster 2 had high levels of HIF1α, AMPKα, cIAP (Fig. 1F-H), cleaved Caspase-9, Caspase 9(D1315 and D330), and APAF1 (Figure S5 A-D). Cluster 3 had significantly higher levels of VEGFR2 (Fig. 1I), while cluster 4 showed increased expression of mTOR (Fig. 1J). Pairwise t-test with ANOVA, Tukey HSD p-adjust. < 0.05 was applied to compare all the clusters.
We repeated clustering to group proteins with similar expression across samples. The heatmap (Figure S6) showed normalized protein expression and clinical parameters of all patients. We investigated how clinical factors affected protein levels. Patients of different sex, age, and GBM subtypes were evenly distributed across clusters. LTS patients were found in clusters 1, 2, and 3 (26.66%, 24.13%, 38.09%), while STS samples were mainly in clusters 2 and 3 (17.24%, 28.57%), but not in clusters 1 and 4.
Cluster-specific survival was analysed using patients' OS time in months (Fig. 2). Although cluster 4 had a trend towards shorter OS compared to clusters 1–3, no significant differences were found between clusters (pairwise t-test, p = 0.08) (Fig. 2A). Silhouette analysis indicated that clusters 1 and 4 were the most clearly defined clusters, with the highest separation (Fig. 2B). Cluster 1 had a silhouette coefficient of 0.75, and cluster 4 had a coefficient of 0.85, while clusters 2 and 3 had coefficients of 0.43 and 0.62, respectively.
Elevated levels of GAB1 (Y627), SRC (Y527), BCL2 (S70) and RAF (S338) phosphoproteins associate with OS and are differentially expressed in LTS and STS
To investigate the association of individual proteins with OS, we fitted 72 univariate Cox regression models across the expanded-GT cohort (N = 126). Six proteins significantly correlated with OS [p27, GAB1 (Y627), SRC (Y527), BCLXL, BCL2 (S70), and RAF (S338)] (Table S1; Likelihood ratio p-value < 0.05). We did not find any significant protein when adjusted for multiple comparisons.
As expected, we also identified significant differences in median levels of all six proteins between STS and LTS samples when analysed via pairwise t-test (Fig. 3A-F). Median levels of phosphorylated SRC(Y527) (p = 0.010; -0.15 and 0.49 median level in LTS and STS), GAB1 (p = 0.005; 0.45 and 0.01), BCL2 (p = 0.015; -0.49 and 0.12), RAF (p = 0.03; -0.5 and 0.01) and BCLXL (p = 0.010; -0.49 and − 0.46) expression were significantly lower in samples of patients with LTS compared to STS samples. In contrast, we found significantly greater median levels of P27 expression (p = 0.002) (-0.153 median level) in LTS samples compared to STS.
Rna-seq Transcriptomics Analysis In Sts And Lts
We analyzed gene expression variations between STS and LTS patient tumors to create transcriptomic signatures for defining survival groups. Our analysis identified 1577 differentially expressed genes (DEGs) (N = 737 downregulated and N = 99 upregulated) between the two groups (Figure. 4A). Among these, C6, OTX2, and DAZL were the most differentially expressed down-regulated genes in LTS, while RAX and ISL1 were the most differentially expressed up-regulated genes in STS.
GO enrichment analysis was used to identify enriched pathways for DEGs in LTS and STS samples. GO Biological process (BP) terms were mainly associated with cilium gene ontologies such as cilium movement, microtubule bundle formation, cilium assembly, and cilium organization for extreme responders (p.adjust < 0.01) (Fig. 4B). Upregulated DEGs in STS were enriched with terms such as animal organ formation, regulation of blood coagulation, and regulation of homeostasis (p.adjust < 0.006), among other developmental terms (Fig. 4C). Conversely, downregulated DEGs in STS were highly enriched in microtubule bundle formation, cilium movement, organization and assembly, as well as cilium and flagellum dependent cell motility (p.adjust < 2e-04) (Fig. 4D).
A Cilium Gene Signature Is Prognostic Within The Gt Cohort
As GO analysis revealed most genes related to cilium annotations, we assessed the relationship between patient survival and cilium gene expression identifying N = 44 genes involved in cilium gene ontologies. Survival analysis using on these 44 genes of the-expanded GT cohort, revealed an improved OS in tumours with higher cilium gene signatures expression (p < 0.001) (Fig. 5; 5A).
Since cilium gene signatures appear to be a positive prognostic maker in the expanded GT cohort, we performed IHC to assess ciliae presence in a representative cohort of LTS patients. We analysed three STS samples (GT.02.01, GT.02.02, GT.02.05) (Fig. 5B) and three LTS samples (GT.02.27, GT.02.39, GT.02.42) (Fig. 5C) for ciliae expression using anti-acetyl-alpha Tubulin as marker of ciliae [19]. Indeed, two of the three LTS samples analysed showed a strong presence of ciliae, which was not observed in any of the STS tumour samples.
Master Regulator Analysis Revealed Signalling Proteins And Tfs Associated With Upregulated Sts Samples
We used a Master Regulator analysis to identify new targetable signaling pathways by examining 1577 DEGs and locating clusters of TF binding sites in upstream regulatory regions. We identified 263 TFs and enhancers targeted by them, and then used CMA to identify two modules containing 13 TFs regulating regions of our genes of interest[17]. Using the TRANSPATH® database [18], we reconstructed the signaling network, identifying 25 distinct master regulators and 13 TFs in our upregulated gene network in STS samples in comparison to LTS samples (Figure S7). These could be novel targets for inhibiting overactivated signaling pathways in GBM. Table S2 shows the lists of identified master regulators and their associations with pathways, such as integrin signaling and cell cycle regulation. We did not find any significantly enriched TFs in promoters of downregulated genes in the STS versus LTS comparison, in accordance with previous workflow publications [22].