Generation and functional characterization of chemotherapy-resistant cells

Due to the lack of suitable models for studying the mechanisms of drug resistance in MB, we setup chemotherapy tolerant MB models (MB Resistant, MB-R) by weekly treating (see Materials and Methods for more details on treatment schedule) naïve/sensitive MB (MB-S) cell lines (DAOY and HD-MB03) and primary cultures (HuTuP33) with a drug cocktail consisting of four of the most commonly used chemotherapeutics for MB patient therapy [2], including Vincristine, Etoposide, Cisplatin and Cyclophosphamide, namely VECC (Fig. 1A and Additional file 1: Fig. S1A, B). VECC cycles-exposed cells were routinely tested until acquisition of drug resistant traits (commonly after 6–9 treatment cycles) through a resazurin-based cell viability assay. MB-R cells displayed a significant greater resistance to a 72 h VECC exposition (~ 3 to 4 folds in terms of GI50 shift) when compared to their relative MB-S counterparts (Fig. 1B), nevertheless showing no evident morphological changes associated to the cycling treatment schedule (Additional file 1: Fig. S1C). Once established and validated, MB-R cells were then functionally characterized and, considering the relative selectivity of most chemotherapeutics toward highly dividing cells, primarily tested for their proliferative potential. To this end, MB cell proliferation was examined through both trypan blue exclusion and EdU incorporation assays which concordantly disclosed identical growth rates to MB-S cells, demonstrating that the observed reduced VECC response was not dependent on differential proliferation (Fig. 1C and Additional file 1: Fig. S1D, E).

Fig. 1

Functional characterization of chemotherapy resistant MB-R cells. A Drawing summarizing the VECC treatment schedule (one VECC exposure per week for 24 h) to which MB cell lines and primary cultures (MB-S) have been subjected in order to eventually acquire drug resistant traits (MB-R). B Dose–response curves displaying increased resistance to VECC treatment (72 h) of MB-R models relative to their naïve/sensitive counterparts (MB-S). GI50 ratio: MB-R GI50/MB-S GI50 (n = 6). C Growth curves showing that both MB-S and MB-R cells possess equivalent proliferation (within 96 h) as measured by trypan blue exclusion assays (n = 3). D Limiting dilution assays comparing self-renewal potential of MB-S and R cells. Initiating cell frequency F of cells is reported. p calculated by extra sum-of-squares F test (n = 4). E Representative cytofluorimetric plots displaying functional assessment of MB-S/R drug extrusion potential by staining with Rho 123. Control unstained cells are reported as solid grey distributions

According to the “cancer stem cells (CSC) paradigm” and the knowledge that they may display refractoriness to anti-cancer agents by a potential enhanced drug efflux capacity [7], we analyzed both MB-S and MB-R cells for their self-renewing potential through limiting dilution assays. Indeed, HD-MB03-R and HuTuP33-R cells were characterized by a significant increase of self-renewal, with DAOY-R cells displaying no deviations from MB-S cells, probably due to a remarkably high “initiating cell” frequency even in basal conditions (Fig. 1D). However, this increased CSC phenotype was not accompanied by a concomitant enhancement of MB-R cell drug extrusion potential, as verified by a comparable functional efflux of Rhodamine 123, a shared substrate of several multidrug resistance proteins (Fig. 1E and Additional file 1: Fig. S1F). Collectively, these results suggest that additional mechanisms requiring a differential regulation of multiple intracellular pathways, may be involved in sustaining MB-R chemotolerance.

Proteomic and transcriptomic profiles converge to the identification of relevant pathways fueling MB-R resistance

Biological processes, in both healthy cells and disease, involve a highly dynamic and interactive network of molecular players, including transcripts, proteins, and metabolites, which participate at different levels in defining cell phenotype and behavior. Based on our initial results which preliminarily excluded that the observed MB-R chemoresistance could be dependent on the most common non-mutational mechanisms [6], we hypothesized that treatment resistance in MB would rely on a complex coordination of several biological processes. For this reason, we implemented a multilevel approach based on the integration of both the proteomic and transcriptomic profiles derived from all MB-S/R cells considered in this study.

Starting from proteomics, matched MB-S/R cells were analyzed through a label-free LC–MS/MS approach, identifying, on average, 1409 ± 370 proteins per sample (for a total of 3506 proteins across all the analyzed MS runs) with a consistent protein expression across biological replicates (Additional file 1: Fig. S2A, B). Moreover, Principal Component Analysis (PCA) disclosed an evident differential distribution of MB-R cell samples relative to their sensitive counterparts (Fig. 2A). Accordingly, we identified 268 shared proteins (Differentially Expressed Proteins, DEPs) displaying a significant differential abundancy, based on their label-free quantification (LFQs) values, between MB-S and MB-R cells (Fig. 2B and Additional file 2: Table S1). With the aim to uncover any potential differentially regulated pathway, we performed String analysis on significant DEPs that allowed their categorization into four major protein–protein interaction clusters through k-means String Network clustering algorithm (Additional file 1: Fig S2C and Additional file 2: Table S1). Pathway enrichment analysis (Fig. 2C) using KEGG and Reactome databases then revealed that chemotherapy resistant cells share an altered expression of proteins involved in multiple processes including: central carbon and nucleoside metabolism (Cluster 1), protein synthesis and post-translational modifications (Cluster 2), mRNA metabolism and transport (Cluster 3), and cytoskeletal organization and intracellular trafficking (Cluster 4), thus suggesting a multi-faced regulation of pro-survival processes upon acquisition of resistant traits in MB cells.

Fig. 2

Proteomic and transcriptomic characterization of MB-R/S models. A PCA displaying distribution of MB-S/R samples from DAOY (■, n = 5) and HD-MB03 (●, n = 3) cell lines, and HuTuP33 (▲, n = 3) primary cells based on label-free MS-based protein expression data. B Heatmap reporting normalized protein abundancy of differentially expressed proteins (DEP, p < 0.05, n = 268 by t test) between MB-S and MB-R samples as in A; FC: fold change. C Pathway enrichment analysis (within the Gene Ontology (GO) terms) of DEP clustered into four distinct networks identified through STRING protein–protein interaction analysis (see also Additional file 1: Fig. S3C). Proteins significantly contributing to the reported GO enrichments are highlighted by colored dots. Cluster 1 (shades of blue): n = 86, Cluster 2 (shades of green): n = 85, Cluster 3 (shades of yellow): n = 51 and Cluster 4 (shades of red): n = 39 (Supplementary Table S1). D PCA displaying distribution of MB-S/R samples from DAOY (■, n = 4) and HD-MB03 (●, n = 4) cell lines, and HuTuP33 (▲, n = 4) primary cells based on Clariom S™ arrays-based whole transcriptome data. E Enrichment map summarizing GSEA results performed in the Hallmarks, C2cp, and C5bp MsigDB gene set collections (q < 0.1 and p < 0.05); NES: Normalized Enrichment Score. F Chord plot highlighting a high concordance between the enriched cellular processes/pathways identified through proteomic and transcriptomic analyses

In order to integrate these results into a more comprehensive intracellular network, a similar set of samples from MB-S/R cells was subjected to whole transcriptome analysis through Clariom™ S Affymetrix arrays. This analysis confirmed a clear-cut differential transcriptional behavior across MB-S and R samples (Fig. 2D) and identified 962 differentially expressed genes (DEGs) in MB-R relative to MB-S cells (445 up- and 517 down-regulated; Additional file 3: Table S2). This result corroborates the hypothesis that the acquisition of drug-resistance may be sustained, at least in part, by a peculiar transcriptional shift occurring in MB cells once they are exposed to a recurring chemotherapeutic stimulus. Mirroring the above-described approach of proteomics-driven interpretation of drug resistance, we performed a Gene Set Enrichment Analysis (GSEA) comparing MB-R with MB-S in each cell type. Common enriched pathways from the Hallmarks, C2cp and C5bp MSig databases highlighted a MB-R associated transcriptional dysregulation of several metabolic processes, comprising an increase of nucleoside and fatty acid metabolisms, protein translation, inflammation and immune system-related signaling, and cytoskeleton-dependent processes, including intracellular trafficking and cell adhesion/migration (Fig. 2E). Moreover, MB-R cells were characterized by a strong downregulation of neural differentiation genes (Fig. 2E), in accordance with the increased self-renewal of MB-R cells, already reported in Fig. 1D. Trying to integrate proteomic and transcriptomic information, a comparative analysis demonstrated a dramatic convergence between the identified altered pathways, that portrays MB resistant cells as endowed with a strong dysregulation of several metabolic functions (from biosynthetic pathways to cytoskeleton-dependent processes), remodeled gene transcription and protein translation activities, and increased activation of pro-inflammatory and stress responses (Fig. 2F).

Finally, in order to verify if the observed proteomic and transcriptional changes could be dependent on any recurrent rough genomic aberration accumulating as a consequence of the VECC-dependent DNA-damaging environment, we investigated if genomic copy number evolution (CNE) had occurred in MB-R cells. Indeed, CNE arising in cancer cells during treatments has been widely reported and peculiar therapy-selected cellular clones already recognized in several recent studies [22,23,24,25]. Nevertheless, through a methylation array-based identification of MB-S/R DNA copy number alterations, we did not retrieve any recurrent VECC-induced chromosomal aberration across the MB-R models analyzed. In particular, we could detect some additionally deleted regions in DAOY-R chromosomes 2 and a restricted 15q deletion in HuTuP33-R cells (Additional file 1: Fig. S3A, C). Conversely, HD-MB03-R were characterized by a negative selection of peculiar clones bearing discrete 2q and 6q deletions (Additional file 1: Fig. S3B). However, the “classical” pediatric brain tumors-associated genomic alterations (such as MYC amplifications or TP53 deletions [20, 26]) were maintained in MB-R cells and we could not identify any evident correlation between retrieved CNEs and dysregulated genes (not shown). These observations suggest that acquisition of drug resistance and the concomitant observed transcriptional shift is not dependent (or very partially) on CNE in our MB-R models.

MB patients correlating with a chemotolerant transcriptional phenotype are characterized by worse prognosis

Data generated so far, highlight a peculiar rewiring of the transcriptional and proteomic profile of MB-R cells compared to their MB-S fellows and reveal several dysregulated processes potentially contributing to the response to chemotherapy and progressive acquisition of resistance. Based on this information, we wondered if these models could be representative of a subset of highly aggressive and possibly chemotherapy-resistant MB tumors. To verify this hypothesis, we correlated the expression values of the above identified DEGs in MB-R with their corresponding expression levels in tumors belonging to a large pediatric MB patient gene expression cohort (GSE85217) [20]. Interestingly, correlation analysis allowed the identification of three main patient clusters, endowed with a differential grade of correlation with a MB-R dependent chemotolerant transcriptional profile (Fig. 3A). More importantly, MB patients displaying the highest transcriptional correlation with the MB-R chemotolerant state (Cluster 3), were also characterized by a significant worse prognosis in terms of overall survival, if compared to low-correlating tumors (Cluster 2; Fig. 3B). Moreover, when subgrouping MB patients according to their transcriptional subtype classification [15, 26], this differential outcome became even more clear, with Cluster 3 patient displaying a significant more unfavorable outcome in the SHH and Group 4 subgroups (MB patients belonging to the WNT subgroup were necessarily excluded due to their small number and extremely favorable prognosis). Of note, despite not statistically significant, also Cluster 3 patients belonging to Group 3 showed a worsened outcome (Additional file 1: Fig. S4A). These results highlight a strict correlation between our MB-R models and the most aggressive, and possibly therapy resistant, MB tumors retrieved from patients, confirming the relevance of our in vitro setting. Of note, the robustness of these models was supported by the highly consistent data obtained across diverse MB cell lines and primary cultures, and the relative independence from the molecular classification assigned to MB tumors.

Fig. 3

MB patients displaying high correlation with a MB-R associated transcriptional phenotype are characterized by worst prognosis. A Heatmap displaying the existing correlation (Spearman r) between MB-R models used within our study and MB patient samples from the GSE85217 dataset [20], based on the expression of previously identified DEGs between MB-S and MB-R models. MB patients have been sub-grouped into 3 different clusters by K means clustering (Euclidean distance). B Kaplan Meyer curves comparing the available overall survival (OS) data of MB patients belonging to each Cluster. y: years

Kinome activation profiling reveals a complex multi-level regulatory network sustaining resistance

Once verified the applicability of our MB-R models in the characterization of MB resistance to treatments, we evaluated the kinome activation status of MB-R cells and characterized any potential kinase activity deviation relative to sensitive cells. To this end, cell lysates from both MB-S and MB-R cells were subjected to phospho-tyrosine (PTK) and serine-threonine kinome (STK) activity evaluation through the PamGene® technology. Of the 196 PTK and 144 STK phosphorylation sites included into the PamChip® arrays, 17 PTK and 60 STK target peptides resulted as differentially phosphorylated between MB-R and S cells (not shown). According to these data, an upstream kinase analysis-based identification of potential differentially activated kinases in MB-R cells suggested an increased activity of SRC family (FGR, LCK, LYN, SRM), ErbB family (Her2, Her3, Her4), and TAM family (TYRO3, AXL, MER) kinases, and a significant over-activation of PDGFR, IGFR, FAK, Ryk and RON. On the contrary, MB-R cells displayed a reduced activity of kinases involved in calcium signaling (CAMK2A, CAMK4, CK2α1), AMPK/AKT/mTOR axis (PKAα, AMPKα1, AKT1, P70S6K), NF-kB pathway (IKKα and IKKβ), PKGs, and kinases related to cell cycle control including some CDK family members (CDK10, CDK14, CDK15) and CHK2 (Fig. 4A, left panel).

Fig. 4

Kinase activation status of MB-S/R cell models. A Barplot summarizing the significant differential activation status of indicated kinases (both PTK and STK identified through the PamGene’s UKA) (left panel) and indication of their relative contribution (grey slots; right panel) to gene ontology term enrichments (green lines) reported in B. B Commonly enriched gene ontology processes between differentially activated kinases (together with their relative differentially phosphorylated peptides; green) reported in A and DEGs from gene expression analysis (blue). Positive/negative enrichments are indicated by a differential red/blue coloring of panel sections above

Enrichment analysis based on combined upstream identified kinase activation status and phosphorylation of target peptides was performed through GSEA (Additional file 1: Fig. S5) and integrated with the already available transcriptional enrichments to highlight which relevant pathways, collaboratively affected by both kinase-dependent signaling and transcriptional modulations, were significantly associated with the observed drug resistant phenotype. Overlapping enrichments summarized in Fig. 4B suggest that kinases differentially activated in resistant cells cooperatively affect cell adhesion and migration, control cell differentiation, regulate immune system processes, and enhance the cell response to stress, eventually influencing cell survival. These results further corroborate the hypothesis that the chemotolerant state displayed by MB-R models is sustained, at different levels, by a complex regulatory network involving the concurrent tuning of kinase activation, gene transcription, and protein production.

High-throughput screening identifies nucleotide metabolism as a functional vulnerability of MB-R cells

The applied multiomic approach defines a significant deregulation of several pathways in resistant MB cells, converging to cell metabolism, RNA/protein homeostasis, and immune response, thus significantly impacting on patient outcome. To identify the specific dysregulated pathways, that functionally sustain the chemotolerant phenotype displayed by MB-R cells, both sensitive and resistant cells were screened for their sensitivity/resistance to a large library of compounds, characterized by widespread targets and mechanisms of action. In particular, the implemented HTS approach consisted in a resazurin-based cell viability screening of 3533 compounds from different commercially available drug libraries, including 1280 FDA-approved drugs (Additional file 1: Fig. S6A). To balance the cohort of our MB-R models, we generated and characterized (according to the same treatment schedule and analytical tools) chemotherapy resistant counterparts of Med-411 primary MB cells, which displayed a consistent behavior with previously described MB-R models in terms of retained morphology, VECC GI50 shift, proliferation, self-renewal, drug extrusion capacity, and CNE (Additional file 1: Fig. S7).

A primary screening was performed by treating MB-S/R cells for 72 h at a fixed 5 µM dose for each compound. After data normalization, the robust Z score was used to rank all the tested compounds according to their capacity to inhibit cell growth or cell viability [21]. The top 5th percentile of ranked compounds derived from each screened MB model (cumulatively n = 303) was then additionally tested through 6-points fivefold dilution dose–response experiments in both MB-S and MB-R cells identifying 56 out of 303 compounds (18.48%) with a selective action against one or more MB-R cell models (Fig. 5A, red shaded dots). On the other hand, 23 out of 303 (7.59%) compounds displayed increased activity on at least one MB-S model (Fig. 5A, green dots). Classification of 5th percentile ranked compounds in definite subgroups according to their pharmacological action or target, allowed to easily recognize at least two main drug clusters displaying high selectivity towards MB-R cells, classifying into the “antimetabolites” and “glycosides” compound families. Given the high toxicity displayed by the selected glycosides when administered to normal human cells (Additional file 1: Fig. S8), we excluded this class of drugs and its related biological processes from further analysis. It is well known that most antimetabolites interfere with the metabolic processes responsible for the correct availability of purine or pyrimidine nucleotide precursors, by mildening their synthesis or competing with them during nucleic acid synthesis. Considering that, among the transcriptional changes related to the acquisition of resistant phenotype, different processes associated to the metabolism of nucleosides were specifically enriched (Fig. 2E), we performed an additional GSEA analysis aimed at a more detailed identification of specific enriched gene sets within a collection of antimetabolites-related processes. This confirmed that both purine and pyrimidine metabolic processes are significantly up-regulated in chemotherapy resistant cells (Fig. 5B), suggesting that nucleoside metabolism may represent a potential vulnerability of MB-R cells.

Fig. 5

Functional identification of MB-R vulnerabilities by HTS. A Distribution of 303 active compounds (5th percentile of all ranked compounds), selected from primary HTS, tested in dose–response experiments, and arranged according to their target/mechanism of action. Each compound is represented as a coloured dot, with increased potency against n MB-R models represented by progressively darker red and increased activity against MB-S cells shown by green color. B Enrichment map summarizing transcriptional GSEA results (p < 0.05) performed in a subset of Gene Set terms (across the C2cp, and C5bp MsigDB gene set collections) correlated to the antimetabolites function, grouped according to their similarity

Antimetabolites sensitize MB-R cells to chemotherapy

Our experimental approach allowed to functionally select nucleoside metabolism as a relevant cellular process characterizing the resistant phenotype of MB cells which could be then exploited as a specific susceptibility of these cells, with evident therapeutic purposes. Starting from these premises, we expanded the antimetabolite compound collection by including additional molecules possibly excluded from the previous analyses since (i) absent in the HTS library, or (ii) ranked below the 5th percentile. A total of 41 antimetabolites were then tested in dose–response experiments in the two available primary MB-R models. Among all tested compounds, 16/41 (39.02%) demonstrated a significantly increased (GI50 ratio > 2) activity against at least one primary MB-R culture (Fig. 6A), thus strengthening the results obtained from the primary HTS. Indeed, by comparing the normalized GI50 ratio between MB-S and MB-R pairs, it became particularly clear that the purine analogues class of antimetabolites displayed a remarkable efficacy against MB-R (Fig. 6A and Additional file 1: Fig. S9A). Interestingly, in accordance with the activity of antimetabolites, we found that purine analogues targets were overexpressed in MB-R cells, together with folate metabolism genes, instead of pyrimidine analogues target genes, which rather displayed reduced expression (Fig. 6B). In particular, the adenosine analogues clofarabine (GI50 ratio of 5.60 ± 3.38 in HuTuP33 and 18.95 ± 0.07 in Med-411), cladribine (GI50 ratio of 18.84 ± 2.44 and 3.54 ± 0.01), and fludarabine (GI50 ratio of 3.92 ± 1.84 and 5.27 ± 2.46) and the guanine analogue 8-azaguanine (GI50 ratio of 4.22 ± 0.8 and 14.75 ± 3.91) showed a remarkable increased activity against MB-R cells (Fig. 6C and Additional file 1: Fig. S9A, B). As expected, the prodrug nelarabine was inactive in both MB primary cultures, regardless of their resistance status [27]. Of note, almost half of antimetabolites were not toxic in MB cells (within the tested dose range). For these compounds, the GI50 values were not calculated due to lack of cell killing/antiproliferative effect (< 50% at the highest dose).

Fig. 6

Nucleoside analogues display increased activity in MB-R cells. A Distribution of 41 different antimetabolites according to their pharmacological classification and target. Each tested compound is represented as a progressively red dot when displaying increased activity against MB-R primary cultures. B Heatmap displaying the differential expression, in MB-R cells, of a subset of genes serving as targets for the antimetabolite drugs tested within the study. C Representative dose–response curves of the most active purine analogues against MB-R primary cultures relative to their MB-S counterparts. D HSA synergy matrixes calculated for the selected purine analogues as in C when combined with VECC in both HuTuP33-R (top panels) and Med-411-R (bottom panels) primary cultures. Compounds have been considered synergistic when HSA ≥ 10, antagonistic when HSA ≤ 10, and non-interactive for 10 > HSA > -10. Relative dose–response matrixes are reported in Supplementary Figure S10

To further increase the potential translational relevance of our functional findings, we conducted drug synergism studies by combining the previously identified most effective purine analogues fludarabine, cladribine, clofarabine, and 8-azaguanine with the VECC chemotherapeutic cocktail in a two-fold dilution matrix layout and evaluated their potential synergism by calculating HSA synergism scores [28]. All tested compounds displayed a significant synergistic action when combined with VECC (HSA score ≥ 10; Fig. 6D and Additional file 1: Fig. S10), suggesting that the combined use of antimetabolites together with the classical MB chemotherapeutics could represent a successful strategy to sensitize MB resistant cells to conventional chemotherapy.