Introduction

Pediatric Diffuse Midline Gliomas (pDMG), including Diffuse Intrinsic Pontine Gliomas (DIPG) are rare and aggressive brain tumors that arise in the pons, thalamus or spinal cord of children, most commonly between the ages of 5 and 10^1–3. pDMG are almost uniformly fatal, with a median overall survival of 9 to 11 months^4,5, thereby representing the leading cause of mortality in pediatric neuro-oncology^1,3. Clinical management of pDMG and especially of DIPG is a major challenge given their location in vital nervous centers and their leptomeningeal dissemination, which prevent any prospect of surgical intervention⁶. Radiotherapy, the current standard-of-care, is at best only transiently effective^3,7. Moreover, pDMG are highly resistant to currently available chemotherapies although promising combinations of molecules, including ONC201, are currently in clinical trials⁸.

High-throughput sequencing demonstrated that pDMG are associated with a major disruption of the epigenetic landscape, resulting in 80% of cases from a lysine-to-methionine substitution at position 27 (K27M) in genes encoding histone variants H3.3 (H3F3A/H3-3A) or H3.1 (HIST1H3B/H3C2 or HIST1H3C/H3C3)^4,9–11. H3K27M variants, also known as oncohistones, are affine for EZH2, the methyltransferase subunit of Polycomb Repressive Complex 2 (PRC2) and inhibit its activity. As a result, pDMG show abnormal methyl group deposition on H3K27 and a global loss of trimethylated H3K27 (H3K27me3)^12–17, a histone mark associated with transcriptional repression.

In addition to TP53, other genetic alterations have been described, notably in the PDGFR, EGFR or PIK3CA genes¹, suggesting an oncogenic synergy between H3K27-based epigenetic remodeling and the activation of several transcriptional programs^18,19. Accordingly, 20 to 30% of pDMG cases are associated with mutations in the ACVR1 gene, encoding for the bone morphogenetic protein (BMP) type I receptor ALK2^20–23, which leads to the overactivation of intracellular BMP signaling pathway^24–32. Such increase in BMP signaling in pDMG epigenetic context has been suggested to promote tumorigenesis by maintaining cells in a proliferative, mesenchymal-like and undifferentiated state in vitro and in vivo^24,25,27,29. In 69% of cases, these ACVR1 mutations are associated with the H3.1K27M mutation, whereas less than 20% are observed in H3.3K27M tumors¹. The question is then whether BMP signaling pathway could also be involved in the tumorigenesis of H3.3K27M-positive tumors. Recent data support the idea that H3.1K27M tumors mutated for ACVR1 would emerge from a ventral pool of oligodendrocyte progenitor cells (OPC) characterised by the expression of the transcription factor NKX6-1 and dependent on the Sonic Hedgehog (SHH) signaling²⁹. In contrast, H3.3K27M are thought to preferentially derive from dorsal progenitors expressing the dorsal PAX3- and BMP-dependent transcription factor²⁹. In other words, H3.1K27M tumors would depend on the concomitant acquisition of the ACVR1 mutation for transformation, whereas the BMP pathway could be activated in a tumor-independent manner in H3.3K27M tumors, due to their BMP-rich microenvironment. However, while the pro-tumorigenic activity of BMP signaling is clear in ACVR1 mutant pDMG, recent data have unexpectedly reported that BMP ligands may exert a tumor-suppressive activity in H3.3K27M-ACVR1 wild-type (WT) pDMG cellular models³³.

Here, we have integrated bulk, single-cell and spatial transcriptomic data from patients with functional approaches in cellular models to characterize the impact of BMP activation in H3.3K27M-ACVR1 WT pDMG. Bioinformatic analyses indicate that BMP signaling pathway activation ground state is independent of ACVR1 status in pDMG, which likely results both from BMP2/7 tumor-autonomous and microenvironment-driven signals in ACVR1 WT tumors. Functional modeling on pediatric glioma cell lines and spatial transcriptomics unveil that H3.3K27M and BMP2/7 synergize to induce a transcriptomic switch leading to a quiescent but invasive cell state. These results shed new light on the complex phenotype resulting from the synergy between activation of the BMP pathway and epigenetic remodeling induced by the H3.3K27M mutation in pDMG, and the interest of a therapeutic approach targeting the downstream oncogenic pathways responsible for the invasive potential of tumor cells.

Results

BMP pathway activation level is independent of ACVR1 mutational status in pDMG, and likely supported by BMP2/BMP7 in H3.3K27M tumors

Oncogenic^24,25,27,29 and tumor-suppressive³³ functions of the BMP pathway have both been reported in pDMG. To clarify this point, we first compared tumor clustering according to their ACVR1 status, based on the transcriptional profiling of 193 BMP pathway target genes as a reliable readout of the pathway activation state, including notably the downstream transducers ID1, ID2, RUNX2 and GATA3 (see list in Supp. Table 1). Principal component analysis (PCA) performed on this BMP target genes subset demonstrates no segregation of ACVR1 mutant (in blue) from ACVR1 WT samples (in red) in two independent publicly available transcriptomic cohorts of 31³⁴ and 40¹ patients (Figure 1A and Figure 1—figure supplement A, respectively). Consistently, the expression pattern of these genes does not segregate pDMG according to their ACVR1 status in an unsupervised hierarchical clustering analysis in both cohorts (Figure 1B and Figure 1—figure supplement B). To have a more comprehensive overview of BMP pathway activation level, we used the PROGENy method to infer downstream TGF-β/BMP response footprint from perturbation-response genes indicative of the pathway activity³⁵. Interestingly, we observed no significant difference in the TGF-β/BMP activation PROGENy score between ACVR1 mutant and WT tumors (Figure 1C and Figure 1—figure supplement C). Of note, although variable, this score was even higher than the maximum observed in ACVR1 mutant tumors in 12.5 (cohort 2) to 18.5% (cohort 1) of ACVR1 WT pDMG.

Figure 1.

To define what alternative mechanism could lead to the activation of the BMP pathway in ACVR1 WT tumors, we profiled the expression of all major BMP pathway components. Compatible with BMP signaling being active, several ligands, receptors, antagonists and intracellular transducers are expressed in ACVR1 WT pDMG, albeit at expression levels very similar to those observed in ACVR1 mutant tumors (Figures 1B, D and Figure 1-figure supplement B, D). By differential analysis, we found that only CHRDL1 is significantly overexpressed in ACVR1 WT tumors compared to mutant ones in two independent transcriptomic cohorts, as previously reported in pDMG cell lines³³ (Figure 1D and Figure 1-figure supplement D). Although its role in regulating the BMP pathway remains to be fully elucidated^36–40, the gain of CHRDL1 expression is unlikely promoting the activation of the BMP pathway in ACVR1 WT tumors. Since no difference is observed in receptor expression consistently between the two cohorts (Figure 1D and Figure 1-figure supplement D), we hypothesized that activation of this pathway likely results from tumor-autonomous and/or microenvironment-driven production of ligands. Out of all BMP ligands present in pDMG tumors, BMP2 and BMP7 are the two most highly expressed ligands in H3.3K27M pDMG in both cohorts with no significant difference according to H3 and ACVR1 mutational status (Figure 1E and Figure 1-figure supplement E-F). Once induced, BMP ligand production can be maintained by a positive transcriptional regulatory loop^41,42. We then reasoned that in ACVR1 WT pDMG, the priming signal may come from the microenvironment. Accordingly, it has been recently suggested that H3.3K27M DIPG likely occur in cells derived from dorsal PAX3+ BMP-reliant progenitors, and that the oncogenic transformation may result from a crosstalk with BMP ligands present in the microenvironment at that time²⁹. Regulation of dorsal glial cell fate during development has been shown to rely mostly on BMP4 and BMP7^43–46. By analysing data from the ALLEN Developing Mouse Brain Atlas⁴⁷, we observed that BMP7 is only expressed at E13.5 during embryogenesis. It is then progressively re-expressed post-natally from P14 notably in the pontine hindbrain, to reach a maximum in most territories including prosomere 2 at P28, which is compatible with the spatio-temporal window of occurrence of H3.3K27M pDMG tumors²⁹ (Figure 1F). BMP2 is only expressed at the E18.5 embryonic stage and BMP4 expression peaks at P14 corresponding to infancy⁴⁸, before decreasing during the post-natal period in these territories (Figure 1-figure supplement G). Similarly, by integrating transcriptomic data from the HBT program⁴⁹, we observed that BMP7 expression but not BMP2 and BMP4 gradually increases to reach its maximum in the midline structure in the mid-childhood period (green curve, Figure 1G and Figure 1-figure supplement G), then coinciding with the peak incidence of pDMG tumors¹.

Overall, the integration of transcriptomic data reveals that the induction of BMP signaling in ACVR1 WT pDMG, at a level equivalent to that observed in mutant tumors, could notably result from the production of BMP2 and/or 7, initiated by the expression of BMP7 by the microenvironment.

BMP7 synergizes with K27M to induce a transcriptomic program leading to quiescence and invasiveness in a low-grade glioma model

To functionally dissect the impact of H3.3K27M/BMP crosstalk and define whether it may have a rather oncogenic or tumor suppressive value, we first used the two previously described pediatric glioma Res259 and SF188 cell lines, which have been genetically modified to stably express and reproduce the epigenetic context either of WT or mutated forms of the variant H3.3⁵⁰. Interestingly, BMP7 expression significantly decreases after H3.3K27M induction in both SF188 and Res259 cells (Figure 2A), whereas BMP2 and BMP4 expressions are respectively increased or unmodified (Figure 2A and Figure 2-figure supplement A-B). To compensate for that decrease and mimic the high-level of BMP7 expression observed in pDMG tumors, we then assessed the impact of recombinant BMP7 depending on the H3.3 context.

Figure 2.

Using qRT-PCR as a first hint, we observed a significant increase in BMP targets expression following BMP7 treatment in both H3.3K27M-Res259 and -SF188 cells compared to H3.3 WT cell lines, with at least a 2-fold increase at 3 and 24 hours for both ID1 and ID2 (Figure 2A and Figure 2-figure supplement C). Interestingly, both levels and kinetics of SMAD1/5/8 phosphorylation in H3.3K27M versus WT overexpressing cells remain unchanged, indicating that BMP canonical pathway activation is not modified by histone mutational status (Figure 2-figure supplement D). To further investigate and characterize the specificity of the response of H3.3 WT and mutant cells to BMP7, we performed an RNA-sequencing (RNA-seq) analysis after 3h and 24h of treatment with recombinant BMP7. Using PROGENy analysis, we observed that the increase in TGF-β/BMP activation score induced by BMP7 is potentiated and remains significantly higher at 24 hours in H3.3K27M-Res259 compared to their wild-type counterparts (Figure 2B). While differentially expressed genes at 3 hours mostly correspond to the expression of the K27M mutation (Figure 2-figure supplement E-F), a subset of 900 genes appears differentially expressed (DE) specifically between treated and untreated H3.3K27M-Res259 cells, but not in H3.3 WT ones (Figure 2C). Enrichment analyses revealed that DE genes are notably associated with alteration in cell cycle regulation (Figure 2D-E), and that the downregulated ones in H3.3K27M BMP7-treated cells correspond to E2F targets (Figure 2-figure supplement G). Consistently, BMP7 treatment leads to a significant 12.2% and 16.7% increase of cells in the G0/G1 phase respectively in H3.3K27M-SF188 and -Res259 cells, while its effect is limited to a 3% and 10.3% increase in their WT counterparts (Figure 2F). Similarly, the number of H3.3K27M-Res259 cells is significantly reduced by 1.7-fold compared to non-treated cells upon BMP7 treatment, while the decrease is limited to 1.2-fold in WT ones (Figure 2-figure supplement H). Of note, this K27M-dependent BMP7 effect is associated with a significant 1.8-fold increase in cyclin dependent kinase inhibitor 1A (CDKN1A, encoding P21; Figure 2G) expression and a reciprocal decrease in RB1 phosphorylation (Figure 2H). This cell cycle blockade is unlikely to result from the entry of cells in senescence since there is no difference in beta-galactosidase activity (SA-β-gal) between wild-type and mutant cells upon BMP7 treatment (Figure 2-figure supplement I).

To further dissect the crosstalk between BMP7 and the H3.3K27M mutation, we established a decision tree algorithm to specifically isolate which of the 900 DE genes post-treatment BMP7 in H3.3K27M cells correspond to a potentiation of the effect of the K27M mutation by BMP7 (Figure 2-figure supplement J, left panel), or to a specific effect of BMP7 in the K27M context (Figure 2-figure supplement J, right panel). We then pinpointed DE genes that specifically correspond to cooperative effects of K27M epigenetic alterations and BMP7-mediated transcriptional regulation. Interestingly, enrichment analyses revealed that these genes are involved in processes related to invasion/migration, including extracellular matrix organization, regulation of cell/neuron projection and adhesion (Figure 2I). Some of these genes such as ITGA2, ROBO2, and MMP28 are already induced following H3.3K27M expression, and their expression is further amplified by the BMP7 treatment (Figure 2J). Conversely, others, such as COL28A1, are specifically induced or repressed by the K27M+BMP7 context (Figure 2J), consistently with our decision tree algorithm. We then assessed the combined impact of H3.3K27M expression and exposure to BMP7 on the invasive properties of glioma cells using a Matrigel™-coated transwell assay (Figure 2K). Expression of the H3.3K27M mutation is sufficient to drive a moderate increase in invasion compared to the WT context. However, this phenomenon is largely amplified by BMP7, with a global 5.6-fold in the H3.3 mutant context, compared to 2-fold in the H3.3 WT one.

Altogether, these data support the fact that BMP7 is sufficient to induce a transcriptomic reprogramming specific to the H3.3K27M epigenetic context, which leads to the emergence of a quiescent but invasive cell state.

Combined BMP2/BMP7 expression drives a quiescent-invasive tumor cell state in pDMG

Considering the data obtained in the Res259/SF188 mechanistic models, we then sought to define whether this crosstalk between BMP and H3.3K27M was preserved with similar effects in H3.3K27M pDMG models. First, we observed that BMP7 is the most expressed ligand in two ACVR1 WT/H3.3K27M DIPG cell lines, but that few if any other BMP ligands, including BMP4, are (Figure 3-figure supplement A). Unlike in tumors, BMP2 expression is notably low in these cell lines (Figure 3-figure supplement A). BMP2 has been shown to be induced by hypoxia⁵¹ or reactive oxygen species (ROS)⁵². Accordingly, exposure of ACVR1 WT H3.3K27M DIPG cells to hypoxia or ONC201, which is known to significantly increase ROS production^53,54, are both sufficient to induce a significant increase in BMP2 expression (Figure 3A). Thus, BMP2 can be produced autonomously by tumor cells in response to specific stresses. To model the concomitant impact of BMP7 and stress-induced BMP2 in the H3.3K27M mutant background, we analysed the impact of BMP2 addition on DIPG 3D-spheroids. As shown in Figure 3B and Figure 3-figures supplement B-C, treatment of DIPG spheroids with increasing doses of recombinant human BMP2 triggers SMAD1/5/8 phosphorylation and leads to a strong dose-dependent decrease in growth rate. Consistently, Ki67 staining, a marker of proliferation, is significantly decreased upon BMP2 treatment (Figure 3C). Reciprocally, treatment of DIPG-spheroids with the BMP inhibitor LDN-193189 (LDN) leads to a slight increase in KI67-positive cells (Figure 3C). This effect was largely mitigated in BT245 and DIPGXIII cell lines in which the K27M mutant allele was removed¹⁷, indicating that it may depend on the K27M-specific epigenetic context (Figure 3-figures supplement D). Of particular interest, the combined knock-out of K27M and BMP inhibition with LDN treatment resulted in the failure of pDMG cells to proliferate and form gliospheres (Figure 3D). In parallel, we explored the impact of BMP activation/inhibition on DIPG cells migration/invasion propensity. BMP2 significantly increases the migration of tumor cells from Matrigel™-embedded 3D-DIPG spheroids, while an antagonistic effect was observed upon LDN treatment (Figure 3D and Figure 3-figure supplement E).

Figure 3.

To define whether such a BMP-induced quiescent-invasive cell state exists in ACVR1 WT H3.3K27M pDMG tumors, we first used publicly available scRNA-seq data from 10 patients’ biopsies published by Jessa and colleagues²⁹. Integration of all these data unveiled that the BMP-responsive (i.e. High PROGENy TGF-β/BMP score) pool of tumor cells correlate significantly with BMP receptors expression, in particular to BMPR2, BMPR1B, BMPR1A, and ACVR1 levels, respectively in 8, 5, 5 and 4 out of 10 samples (Figure 3E-F and Figure 3-figure supplement F). Interestingly, this BMP-responsive pool of cells is significantly enriched in genes involved in positive regulation of cell migration and in cell-matrix adhesion (Figure 3G), while showing a specific quiescent-compatible decrease in genes involved in transcription/translation processes (Figure 3-figure supplement G). To further define the extent and organization of this pool of cells, we performed spatial transcriptomic analysis of 3 H3.3K27M pDMG tumors. Data integration was performed using pDMG-derived Visium and scRNA-seq signatures recently published by Ren et al. and Filbin et al.^55,56 (Figure 3-figure supplement H). Consistently with observations in bulk and single-cell analyses of patients/models, the most highly expressed BMP ligands are BMP2 and BMP7 along with BMP8B (Suppl. Table 2 and Figure 3-figure supplement H), whose expression is not spatially delimited or associated with a specific gene signature. In line with single-cell analyses, we observed on the contrary that BMP receptors expression correlate significantly with the invasive niche score defined by Ren et al.⁵⁵ (Figure 3I), which can be spatially restricted to histologically delineated areas (samples 1-3, Figure 3I and Figure 3-figure supplement H) or more dispersed within the tumor (sample 2, Figure 3-figure supplement H). Moreover, the invasive niche score strongly correlates to the PROGENy TGF-β/BMP score (Figure 3J and Figure 3-figure supplement I). Coherently, this BMP-responsive invasive niche is characterised by a high expression of key markers of BMP activity (ID1), as well as markers of stemness (CD44, HEY1) and migration/invasion (TNC, IGFBP7, ITGAV) (Figure 3K).

Altogether, these data indicate that tumor-autonomous production of BMP2 and BMP7 synergize to maintain a quiescent-invasive niche in H3.3K27M DIPG.

Discussion

Major advances have been made in understanding the molecular bases of pDMG, among which DIPG, with the identification of major epigenetic remodeling processes induced by histone H3 mutations. However, the oncogenic mechanisms cooperating with these mutations to induce transformation and tumor escape remain largely undefined.

Along this line, a key challenge is to establish the precise role of the transcriptomic reprogramming induced by the BMP pathway, which remains controversial in these pediatric brain tumors^{24,25,27,29,33}. Herein, we performed a comprehensive integration of transcriptomic data, which first supports the view of BMP signaling being also clearly activated in ACVR1 WT/H3.3K27M pDMG. Further analyses will be necessary to define whether the lowest levels of activation correspond to samples with specific locations along the midline. Effect of BMP activation is potentiated by the epigenetic context of these tumors, then leading to a global transcriptional reprogramming. Although it has been previously described that ACVR1-mutant tumors exhibit a higher expression of ID1 and ID2²⁴, the extrapolation of BMP activity from a larger signature of BMP targets and on a score calculated from the inference of gene expression perturbations in response to the TGF-β/BMP pathway indicates the existence of a compensatory mechanism driving BMP activation in non-ACVR1 mutant tumors. By analysing the expression of BMP effectors in pDMG tumors and patient-derived DIPG models, we propose that the activation of this pathway in ACVR1 WT H3.3K27M tumors could be at least partially mediated by two complementary tumor-autonomous and microenvironment-dependent mechanisms. Upon initiation, the expression of BMP, and notably BMP7, could synergize with or even trigger the autocrine production of BMP ligands by the tumors, among which BMP2 and 7. Indeed, such positive feedback loops maintaining the expression of BMP ligands have already been described notably during development^41,42. The fact that the expression of BMP ligands is similar in ACVR1 mutant and WT tumors and independent of H3 status (Figure 1-figures supplement F) suggests that this regulatory loop is involved in most tumors, but probably by different mechanisms depending on the mutational context. Thus, BMP secretion by the microenvironment may prime the BMP activation in tumor cells and be required for oncogenic transformation. Once established, the dynamic modulation of BMP2 expression in response to stresses, such as hypoxia or treatments, could synergize with constitutive production of BMP7 to drive the emergence of an aggressive cell state. Of note, the expression of BMP-target genes ID1 to 4 was previously reported to be strongly decreased in DMG cell lines in which the K27M mutation was removed by CRISPR-Cas9¹⁷, suggesting that even the maintenance of BMP activation in a pool of tumor cells relies on the specific K27M-mediated epigenetic context. In addition, blocking the BMP pathway with LDN in these K27M-KO cells induces tumor cell death, supporting the view that the K27M/BMP oncogenic synergy plays a major role in maintaining oncogenic potential (Figure 3-figures supplement D).

Second, our data are in favor of a rather global oncogenic role of BMP pathway in pDMG gliomagenesis. The tumor suppressor activity of the BMP pathway had been largely extrapolated because of its positive impact on tumor cell quiescence³³. Indeed, using a genetically engineered glioma model, we confirmed that BMP7 is sufficient to potentiate the entry of cells in a quiescent state in a H3.3K27M-dependent manner, via a transcriptomic switch largely relying on the downregulation of E2F-targets cell-cycle regulating genes. Nevertheless, this quiescent cellular phenotype could paradoxically constitute an aggressive treatment-resistant state, as previously observed in adult glioblastoma^57–59 and thus explaining its increase in response to treatment³³. Along the same line, the impact of CHRDL1 increase in H3.3K27M tumors on BMP pathway inhibition may probably need to be qualified³³. Indeed, if CHRDL1 was first classified as a member of the chordin family of secreted BMP antagonists due to sequence homology, it has been shown that it exerts BMP-independent functions in synapses plasticity and maturation³⁹. Moreover, this protein has also been described as an activator of the BMP pathway⁴⁰, suggesting that its high level of expression may not be a robust marker of BMP activation state.

However, it cannot be ruled out also that different ligands of the BMP pathway may have different impacts on the cellular phenotype induced by the H3.3K27M mutation. Accordingly, it was shown that BMP4 treatment promotes differentiation of DIPG tumor cells, in line with its putative tumor suppressor activity³³. However, beyond quiescence, we observed that BMP activation by BMP2/7 in a H3.3K27M epigenetic context induces a transcriptomic switch rather conferring enhanced invasion potential to pDMG tumor cells. Because the level of BMP2/7 is particularly important i) in bulk, single-cell and spatial transcriptomic tumors, ii) in patient-derived cell models, and because the dynamics of BMP7 induction in the post-natal period coincide with the spatio-temporal window of pDMG onset, we believe that the role of the BMP2/7 couple is non-negligible in the pathogenesis of pDMG. Interestingly, this pair has already been shown to trigger BMP signaling as a heterodimer notably during embryogenesis, with a higher efficiency than homodimers^60–62. Keeping in mind that BMP2 appears to be dynamically regulated by tumor cells upon stress, these data suggest that activation of the BMP pathway may be finely regulated in pDMG to be maintained at optimal pro-oncogenic levels, without triggering the tumor-suppressive effects³³ or negative regulatory loops associated with its over-activation^42,63–65.

The next question is whether a therapeutic perspective can be defined by targeting the crosstalk between epigenetic modifications induced by the H3.3K27M mutation and transcriptomic reprogramming induced by BMP2/7. The BMP2/7-driven cell state that we described here fits with previous results obtained in other glioma models, in which a subpopulation of quiescent cells was identified as partially responsible of tumor invasiveness^57,66,67, a hallmark of pDMG aggressiveness⁶⁸. Given the complexity of the phenotype induced by H3.3K27M/BMP crosstalk and the pleiotropic role of BMP proteins in central nervous system, the therapeutic strategy to be developed should be based on the targeting of the downstream effectors, such as ID1 as recently described⁶⁹, or TNC, which may be responsible for the invasive phenotype. Upcoming challenges will be to precisely define the identity of these pro-invasive BMP effectors, to set up a combinatorial therapeutic approach, simultaneously targeting the proliferative compartment and the BMP-responsive H3.3K27M invasive cell state.

Materials and methods

Gene expression analyses of publicly available transcriptomic datasets

Gene-Expression analyses of H3.3K27M-pDMG transcriptomic cohorts

For cohort 1³⁴, HTSeq gene counts and somatic vcf of St Jude’s pDMG samples were downloaded from https://platform.stjude.cloud/. DESeq2 (v 1.36)⁷⁰ was used to normalise the data with the variant stabilisation transformation (vst)⁷¹. ACVR1 mutation status was assessed using tabix⁷² (v 1.15.1) to query the region of the ACVR1 gene chr2:157,736,251-157,876,330 (hg38). Identified variants were manually curated. PROGENy (v 1.18.0)³⁵ was used to infer TGF-β/BMP pathway activity with the “top” parameter set to 100 (default value).

For cohort 2¹, gene expression (z-score) and tumor variants were downloaded on https://pedcbioportal.kidsfirstdrc.org/ from the dataset “phgg_jones_meta_2017”. Samples were filtered based on the following location: brainstem or midline, and the following histone mutation: WT, H3.1K27M and H3.3K27M. Only the 4 datasets with ACVR1 mutant samples were kept for the analysis to avoid biases: “PMID:21931021|PMID:24705251”, “PMID:21931021|PMID:22286216”, “PMID:22389665|PMID:24705252”, “PMID:24705251”. TGF-β/BMP pathway activity was also inferred using PROGENy, but the “top” parameter was set to 178 genes: as half of the top 100 genes of the PROGENy model for the TGF-β/BMP pathway are not covered in the cohort 2, we identified the top 100 genes with the highest absolute coefficient in PROGENy’s model to compute pathway activity with the same number of genes in both cohorts.

For both cohorts, Principal Component Analysis (PCA) was performed using FactoMineR (v.2.4)⁷³ and plotted with factoextra (v1.0.7)⁷⁴.

Patterning of BMP expression in the developing brain

Spatio-temporal gene expression data of BMP2, BMP4 and BMP7 from human developing and adult brain samples were obtained from the Human Brain Transcriptome project (hbatlas.org)⁴⁹. Heatmaps showing relative BMP2, BMP4 and BMP7 expression by in situ hybridisation (ISH) in murine brain across development were obtained from the ALLEN Developing Mouse Brain Atlas (developingmouse.brain-map.org)⁴⁷.

scRNA-seq analyses of H3.3K27M ACVR1 WT DMG

Using the scRNA-seq data published by [Jessa et al., 2022] (GSE210568)²⁹, we selected H3.3K27M-ACVR1 WT pDMGs localized either in the pons or in the thalamus (10 samples). The same processing steps as the authors were performed except for the regression of mitochondrial proportion and the number of UMIs. Normalization was performed with the LogNormalize function of Seurat⁷⁵. All samples were filtered to keep tumoral cells as identified by the authors. As a readout of BMP pathway activity, we used PROGENy to infer the TGF-β/BMP score. We then stratified cells into TGF-β-High and TGF-β-Low groups using the 95th quantile as cutoff. FindMarkers with default parameters was used to identify differentially expressed (DE) genes (FDR <= 0.05) in each group. Enrichment analyses were then conducted using the enrichR package⁷⁶ and the Gene Ontology Biological Process gene signatures⁷⁷, separately on DE genes upregulated in TGF-β-High and TGF-β-Low for each sample. Only the top 10 significant enrichments (FDR <= 0.05) are shown on the plots. Identification of DE BMP ligands and receptors was run with FindMarkers using a logfc.threshold of 0.15 (instead of the default 0.25). For the visualization of all samples on a shared UMAP, Harmony⁷⁸ (v 0.1.1) was used to integrate the 10 samples on the 50 principal components computed. The first 30 harmony-corrected principal components were then used to compute the SNN graph and UMAP.

RNA-sequencing of inhouse pDMG samples

As part as the Share4Kids program, a third cohort was constituted from leftovers DMG samples, obtained through biopsies performed at the Pediatric Hematology and Oncology Institute (iHOPE, Lyon) and the Hôpital Femme Mère Enfant (HFME, Lyon). Tissue banking and research were conducted according to national ethics guidelines, after obtaining the written informed consent of patients. This study was approved by the ethical review board of the BRC of the Centre Léon Bérard (n°BB-0033-00050, N° 2020-02). This BRC quality is certified according to AFNOR NFS96900 (N° 2009/35884.2) and ISO 9001 (Certification N° 2013/56348.2). Biological material collection and retention activity are declared to the Ministry of Research (DC-2008-99 and AC-2019-3426). For RNA-seq library construction of cohort 3, total RNAs from tissues were isolated using the AllPrep® DNA/RNA FFPE kit (Qiagen, 80224) following manufacturer’s instructions. Libraries were prepared with Illumina Stranded mRNA Prep (Illumina, 20040534) following recommendations. Quality was further assessed using the TapeStation 4200 automated electrophoresis system (Agilent) with High Sensitivity D1000 ScreenTape (Agilent). All libraries were sequenced (2×100 bp) using Agilent SureSelect RNA XTHS2 All Exon V8 (Agilent, G9991A)according to the standard agilent protocol.

Quality control of reads was performed using FastQC (v.0.11.9)⁷⁹, followed by trimming of Illumina adapter sequences with Cutadapt⁸⁰ (v.3.4) using the -a CTGTCTCTTATACACATCT and -A CTGTCTCTTATACACATCT parameters. Reads were mapped to the GRCh38 human genome using “two-pass” mode STAR (v.2.7.9)⁸¹, with Ensembl v104 annotations. Gene counts were then computed using HTseq-count (v.0.13.5)⁸² with the following parameters: “--order pos” and “--stranded reverse”. HTseq count files were then loaded in R (v 4.2.0) and two filtration steps were applied using the annotations of org.Hs.eg.db (v 3.15.0)⁸³: genes with low counts (less than 10 reads across samples) were removed; non-protein-coding genes were removed. Filtered gene counts were converted into a DESeq2 object with the design parameter set to account for the tumor histone mutation (HGG/DIPG) and normalized with vst using “blind = FALSE”. The R package ggplot2 (v 3.3.5)⁸⁴ was used to plot the expression of BMP ligands.

Spatial transcriptomics analyses

FFPE tissue sections were placed on Visium slides and prepared according to 10X Genomics protocols. After H&E staining, imaging and decrosslinking steps, tissue sections were incubated with human-specific probes targeting 17,943 genes (10X genomics, Visium Human Transcriptome Probe Set v.1.0). Probes hybridised on mRNA were captured on Visium slides and a gene expression library prepared following 10X genomics dedicated protocol and sequenced on Illumina NovaSeq 6000 targeting 50,000 reads per spot.

Raw reads were pre-processed using the spaceranger count pipeline (v 2.0.0) and the human GRCh38 reference provided by 10X Genomics. Filtered H5 matrices were loaded using the Load10X_Spatial function from Seurat. The following pipeline was applied to process raw counts: counts were normalized with the LogNormalise method and a scale factor of 10000; the top 4000 variable features (“nfeatures” parameter) were identified using the vst method; data were then scaled for all genes, the first 50 components of the PCA were computed using the 4000 variable features identified; SNN graph was constructed using the first 30 dimensions of the PCA (dims parameter); clusters were identified with the Louvain algorithm along several resolutions ranging from 0.1 to 1; UMAP was also computed using the first 30 dimensions of the PCA (dims parameter). Signatures identified by [Ren et al., 2023]⁵⁵ and [Filbin et al., 2018]⁵⁶ were used in conjunction with sc-type⁸⁵ to annotate the clusters at a resolution of 0.5. Annotated clusters were validated with an anatomopathologist, leading to the subclusterisation of the “Invasive niche” of HFME-1 into 2 new clusters (“Invasive_niche_1” and “Invasive_niche_2”). Scoring of Ren’s signatures were made with the AddModuleScore function of Seurat. PROGENy was used to score TGF-β/BMP pathway activity. Markers of the “Invasive niche” were identified with the FindMarkers function using a logfc.threshold = 0.15.

Cell culture and treatments

Human pediatric low-grade glioma cell line Res259 (grade II, diffuse astrocytoma) and high-grade glioma cell line SF188 (grade IV, glioblastoma) were kindly provided by Dr Samuel Meignan. Both isogenic cell lines overexpressing H3.3 wild-type and H3.3K27M were generated as previously described by [Rakotomalala et al, 2021]⁵⁰. All cell lines were grown as a monolayer in DMEM medium with GlutaMAX-I, 4,5 g/L D-Glucose and 110 mg/L pyruvate (Gibco, 31966) supplemented with 10% foetal bovine serum (FBS), 1X MEM non-essential amino acid solution (Gibco, 11140050) and 1X Penicillin-Streptomycin (Gibco, 15140122). Cells were incubated in a humid atmosphere at 37 °C with 5% CO2.

To assess the impact of BMP7 treatment on wild-type or H3.3K27M-mutated Res259 or SF188 cell lines, 1.5×10⁵ cells were plated into 6-well plates in complete medium. After 24 hr, cells were rinsed with PBS and fresh medium and serum-deprived in 1% FBS medium for 3 hours. 50 ng/mL of human recombinant BMP7 (Peprotech, 120-03P) was then directly added to the medium, and cells were harvested at indicated time points for further experiments.

For pDMG cell lines grown as spheroids, HSJD-DIPG-007, HSJD-DIPG-012, HSJD-DIPG-013 and HSJD-DIPG-014 were kindly provided by Dr Angel Montero-Carcaboso. BT245 and SU-DIPGXIII, KO or not for H3.3K27M mutation, were kindly provided by Dr Nada Jabado¹⁷. pDMG cell lines were grown in complete culture medium as described by [Monje et al., 2011]⁸⁶, in 96-well ULA plates (Corning, 7007) or 25cm² low attachment culture flasks (Corning, 431463). Medium was changed twice a week and spheroids were splitted every 1 to 2 weeks when reaching 800 to 1000 µm of diameter using TrypLE Express Enzyme (Thermo Fisher Scientific, 12605010) preheated to 37 °C. Hypoxia and ONC201 treatments were induced on HSJD-DIPG-012 spheroids once a diameter of 600 µm was reached. To induce hypoxia, HSJD-DIPG-012 spheres were incubated for 3 hr in a dedicated incubator at 37 °C, 1 % O2, while the controls were incubated at 37 °C in normoxia conditions. Spheroids were similarly treated or not with 20 µM of ONC201 (Selleckchem, S7963) for 96 hr. Exploration of BMP2 effect on HSJD-DIPG-012 and HSJD-DIPG-014 spheroids was done using 10 ng/mL of human recombinant BMP2 (Peprotech, 120-02C) before protein quantification by western blot.

All cultures were tested every month for mycoplasma using the MycoAlert Mycoplasma Detection Kit (Lonza, LT07-318), in accordance with the manufacturer’s instructions.

RNA-seq data processing and analyses of Res259 cells

For RNA-seq library construction, 1000 ng of total RNAs were isolated using the Nucleospin RNA kit (Macherey-Nagel, 740955) following manufacturer’s instructions. Libraries were prepared with Illumina Stranded mRNA Prep (Illumina, 20040534) following recommendations. Quality was further assessed using the TapeStation 4200 automated electrophoresis system (Agilent) with High Sensitivity D1000 ScreenTape (Agilent). All libraries were sequenced (2×75 bp) using NovaSeq 6000 (Illumina) according to the standard Illumina protocol.

Raw sequence quality was assessed using FastQC (0.11.9)⁷⁹. The trimming step was omitted, as 5’ and 3’ read bases had a quality greater than Q30, and no adapter fragments were detected. Reads were then pseudo-aligned using Kallisto (0.46.2)⁸⁷ with Ensembl v96, human genome build GRCh38. The rest of the analyses were performed using R version (4.2.0). Differential expression (DE) analyses were conducted using the DESeq2 package (1.36.0)⁷⁰ using default parameters. Genes with corrected p-value (Benjamini–Hochberg) < 0.05 and |log2FoldChange (LFC)| > log2(1.5) were considered DE. EnrichR package (3.0)⁷⁶ was used for overrepresentation analysis using the following databases: Gene Ontology (GO) Biological Process 2021, GO Molecular Function 2021, GO Cellular Compartment 2021, KEGG 2021 Human, MSigDB Hallmarks 2020 and ChEA 2016. Enrichment significance was assessed by Fisher’s exact test, and p-values were corrected with the Bonjamini-Hochberg method. Overrepresentation analyses were run separately on upregulated and downregulated DE genes. Principal Component Analysis (PCA) was performed using FactoMineR and plotted with factoextra. Boxplots and heatmaps were respectively made with ggplot2 and ComplexHeatmap (v 2.12.1)⁸⁸ using the vst normalized expression. TGF-β pathway activity was inferred with PROGENy using the first 100 genes of the model (top=100).

qRT-PCR profiling

Total RNA was extracted using Nucleospin RNA kit (Macherey-Nagel, 740955) or RNeasy Plus Micro (Qiagen, 74034) following manufacturer’s instructions. 500 to 1000 ng was reverse transcribed using the iScript cDNA Synthesis kit (Bio-Rad, 1708891) according to the manufacturer’s instructions. Expression of BMP2 [forward primer (5ʹ-TGCGGTCTCCTAAAGGTCG-3ʹ) and reverse primer (5ʹ-GAATTCAGAAGCCTGCAAGG-3ʹ)], BMP7 [forward primer (5ʹ-GGGTGGGTCTCTGTTTCAG-3ʹ) and reverse primer (5ʹ-CCTGGAGCACCTGATAAACG-3ʹ)], ID1 [forward primer (5ʹ-GGTGCGCTGTCTGTCTGAG-3ʹ) and reverse primer (5ʹ-TGTCGTAGAGCAGCACGTTT-3ʹ)] and ID2 [forward primer (5ʹ-CCCAGAACAAGAAGGTGAGC-3ʹ) and reverse primer (5ʹ-GAATTCAGAAGCCTGCAAGG-3ʹ)], were assessed by real-time quantitative QRT-PCR on a LightCycler® 480 instrument (Roche) using the LightCycler® 480 SYBR Green I Master Mix (Roche, 04707516001), according to the manufacturer’s instructions. HPRT1 [forward primer (5ʹ-AAGAGCTATTGTAATGACCAGT-3ʹ) and reverse primer (5ʹ-CAAAGTCTGCATTGTTTTGC-3ʹ)] expression was used as a housekeeping gene.

The expression of a panel of 84 genes of the TGF-β/BMP signaling pathway was assessed in HSJD-DIPG-007, HSJD-DIPG-012 and HSJD-DIPG-014 cell lines by real-time quantitative QRT-PCR on a LightCycler® 480 instrument (Roche) using a RT² Profiler PCR Array (Qiagen, 330231), according to the manufacturer’s instructions. RNA were extracted using the RNeasy mini kit (Qiagen, 74104) and 500 ng of RNA were reverse transcribed using the RT2 First Strand Kit (Qiagen, 330401) according to the manufacturer’s instructions.

Western Blot

Cells were lysed in RIPA Buffer (50 mM Tris-HCl pH 8, 150 mM NaCl, IGEPAL 1%, 0.5% sodium deoxycholate, 0.1% SDS) containing a protease and phosphatase inhibitor cocktail (Thermo Fisher Scientific, 78440). Protein contents were estimated using the Thermo Scientific™ Pierce™ BCA Protein Assay Kit (Fischer Scientific, 23225). Samples were diluted with distilled water to achieve equal concentrations and a loading buffer (4x Laemmli Sample Buffer, Biorad, 1610747) containing 100 mM DTT (Sigma-Aldrich, 11583786001) was added. Protein extracts were then analysed by immunoblot. Briefly, proteins were loaded into SDS-polyacrylamide gels for electrophoresis (Mini Protean TGX gels, Biorad, 4561034) and blotted onto polyvinylidene fluoride sheets (Trans-Blot Turbo Transfer PVDF pack, Biorad, 1704156) using the TransBlot technology (Bio-Rad). Membranes were blocked with 5% BSA in TBS for 1 hr and then incubated overnight at 4° C with anti-pSMAD1/5/8 (1:1000, Cell Signaling Technology, 13820), anti-SMAD1 (1:1000, Cell Signaling Technology, 6944), anti-SMAD1/5/8 (1:1000, Abcam, ab80255), anti-pRb (1:500, Abcam, ab47763), and anti-Rb (1:1000, Cell Signaling Technology, 9309). After three washes with TBS-Tween 0.1%, membranes were incubated with the appropriate HRP-conjugated secondary antibody (1:20.000, Jackson ImmunoResearch). HRP-conjugated β-Actin antibody (1:25.000, Sigma-Aldrich A3854) and HRP-conjugated GAPDH antibody (1:2000 Cell Signaling Technology, 8884), used as loading controls, were incubated for 1 h at room temperature. After three washes with TBS-Tween 0.1%, detection was performed using the Amersham ECL Prime Detection Reagent (Cytiva, RPN2232). Membranes were imaged on the ChemiDoc Touch Imaging System (Bio-Rad).

Proliferation and cell cycle characterisation assays

For proliferation assay, 5×10⁴ Res259 or 7.5×10⁴ SF188 cells were plated into 6-well plates in complete medium. After 72 hr of BMP7 treatment, total cell number and viability were quantified by image cytometry on a NucleoCounter NC-3000 (Chemometec) according to the procedure provided by the manufacturer, using a co-staining of Acridine Orange and DAPI (Chemometec, 910-3013).

For cell cycle analyses, 3×10⁵ Res259 or SF188 cells were plated into 10 cm-petri dishes in complete medium. After 24 hr of BMP7 treatment, 10 μM EdU was added for 1.5 hr. Cells were then harvested using Trypsin (Gibco, 25300054) and washed twice with PBS. Click-it reactions were performed using the Click-iT Plus EdU Alexa Fluor 647 Flow Cytometry Kit (Invitrogen, C10634), according to the manufacturer’s instructions. Cells were counterstained with DAPI and analysed on a BD FACSAria Fusion flow cytometer (BD Biosciences). Data analysis was performed with FlowJo v10.7.1 Software (BD Biosciences). Cells were identified on a Side Scatter (SSC) vs Forward Scatter (FSC) dot plot and cell debris and aggregates were excluded from analysis based on FSC signals.

For β-galactosidase activity determination, 5×10³ Res259 cells were plated into 12-well plates in complete medium. After 72 h of BMP7 treatment, cells were fixed for 5 min in 0.5% glutaraldehyde, rinsed twice with PBS, and incubated for 48 hr at 37 °C in senescence-associated beta-galactosidase (SA-β-Gal) staining solution as previously described⁸⁹. SA-β-Gal+ cells were counted manually, and total number of cells were counted using a DAPI-counterstaining and using the Fiji software⁹⁰.

Invasion assays

For 2D pediatric glioma cell lines (Res259 and SF188), 8.0 µm-Boyden chambers (Falcon, 353097) were coated with 100µL of 10% Geltrex™ LDEV-Free, hESC-Qualified, Reduced Growth Factor Basement Membrane Matrix (Gibco, A1413301) completed or not with 50 ng/mL human recombinant BMP7. During the 3 days preceding the experiment, Res259 WT or K27M cells were pre-conditioned or not with 50 ng/mL human recombinant BMP7. 2.5×10⁴ cells were plated over Geltrex™ in 400 µL 1% FBS-DMEM completed or not with 50 ng/mL human recombinant BMP7. Boyden chambers were deposited in 24-well plates filled with 700 µL of 10% FBS-DMEM. After 48 hr, Boyden chambers were washed with PBS, incubated for 20 minutes with methanol to fix the cells and Crystal violet 0,1% (Sigma-Aldrich, V5265) was used to color cells allowing counting. Pictures of the chambers were taken using an EVOS-7000 (Invitrogen) and analysed using the Fiji software.

For pDMG cell lines grown as spheroids, HSJD-DIPG-012 were plated in 10% Matrigel (Corning, 354277) in Nunc™ Lab-Tek™ Chamber Slide (Thermo Fisher Scientific, 177402). Media and matrigel were supplemented with either 10 ng/mL human recombinant BMP2 or 1 µM LDN-193189 (Selleckchem, S2618). Invasiveness was measured at 48 hr after seeding. Images were captured using an EVOS-7000 (Invitrogen) and analysed using the ‘ROI’ feature in the Fiji software.

Spheroids growth monitoring

To assess the effects of BMP2 on pDMG spheroids, 1×10⁴ cells of HSJD-DIPG-012 or HSJD-DIPG-014 were seeded in 96-well plates (Corning, 7007). Cells were treated at seeding with concentrations of 1, 5, 20, 50 or 200 ng/mL of human recombinant BMP2 (Peprotech, 120-02C). Media were renewed twice a week, with or without recombinant BMP2. Cell growth was monitored using an EVOS-7000 (Invitrogen) and pictures of control and treated 6-days or 9-days spheres were used to calculate sphere area of HSJD-DIPG-014 and HSJD-DIPG-012, respectively.

Immunohistochemistry

HSJD-DIPG-012 spheroids were treated when reaching a diameter of 600 µm either with 200 ng/mL of human recombinant BMP2 or with 1 µM of LDN-193189. Spheroids were then fixed in 4% PFA for 1 h at room temperature, dehydrated and embedded in paraffin. 8 µm-sections were then deparaffinized, rehydrated, and heated in a citrate buffer (0.01 M; pH 6.0) for 15 min. Sections were then incubated overnight at 4 °C with appropriate dilution of anti-Ki67 (1:100, Dako, M7240) in TBS containing Horse Serum (2%). After several washes in TBS-Tween 0.05 %, sections were incubated with the secondary antibody for 1 hr, and then washed again in TBS-Tween 0.05 %. Color was developed with 3,3’-diaminobenzidine (DAB, Vector Laboratories, SK-4100) incubation for 1 to 3 min and with Hematoxylin (Vector Laboratories, H-3401) incubation for 3 min. Images were captured using an EVOS-7000 (Invitrogen) and quantification of Ki67-positive cells was assessed using the “Cell Positive Detection” QuPath function⁹¹, with parameters: “Optical density sum”, “thresholdCompartment”: “Cell: DAB OD mean”, “thresholdPositive2”: 0.25.

Statistical analyses

Data are represented as means ± std. Sample size and replicates are stated in the corresponding figure legends. Using GraphPad Prism 9 software and R (4.2.0), the statistical significance between two groups were determined by one-tailed Mann-Whitney signed-rank tests, apart from the proportion of Ki67 stained cells, which was assessed using Fisher’s exact Test. The following symbols used to denote (<0.05: *; <0.01: **; <0.001: ***; <0.0001: ****).

Data availability

Raw bulk RNA-seq and Visium data generated in this study will be available via GEO. RNA-seq data of pDMG samples from cohort 1 were obtained from St. Jude Cloud⁴⁶ (https://www.stjude.cloud). Cohort 2 of glioma samples can be accessed on https://pedcbioportal.kidsfirstdrc.org/ with the following dataset id: “phgg_jones_meta_2017”. The scRNA-seq cohort of H3.3K27M ACVR1 WT pDMGs is available on GEO (GSE210568) and https://zenodo.org/record/6929428. Data from cohort 3 are described in Supp. Tables 3 and 4.

Code availability

All data analyses’ codes will be made available without restriction upon request.

Author contributions

Study concept and design were performed by P.H., S.Meyer., C.B., M.H., V.R., E.C., and M.C.. Acquisition of data was performed by P.H., S.Meyer., C.B., M.H., A.B.C., A.R., L.M., J.T., N.L., P.Lewandowski., M.C.B, L.B., A.D., I.R., J.Y.B., C.D., V.A., A.M.C., M.L.G., E.P., A.V., P.G.H., S.Meignan, P.Leblond., V.R., E.C., M.C.. Analysis and interpretation of data were performed by P.H., S.Meyer., C.B., M.H., A.B.C., V.R., E.C., M.C.. Drafting of the manuscript was performed by P.H. and M.C.. Writing of the manuscript was performed by P.H., S.Meyer., C.B., M.H and M.C.. Critical revisions of the manuscript for important intellectual content was performed by P.H., S.Meyer., C.B., M.H, A.B.C., V.R., E.C., and M.C.. M.C. obtained funding. V.R., E.C. and M.C. co-supervised the study.

Ethics declarations

Competing interests

The authors declare no competing interests.

Acknowledgements

We thank the patients and their families who consent to participate in this study, as well as charities that support this work, including the Ligue Nationale contre le Cancer, and the Ligue Départementale contre le Cancer de Haute-Savoie et de l’Ain. Our warmest thanks go to Augustine and its Wonder team (“Wonder Augustine” charity) that support this work since the initiation, to Léonie (“Au Pays de Léonie”) and to “Les Etoiles Filantes”, which has funded all the omics analyses. We also thank clinical teams from IHOPe and HFME for their support and contributions, as well as all the facilities from the CLB and CRCL. We are grateful to Séverine Tabone-Eglinger, Loïc Sebileau and Anne-Sophie Bonne for their help with the management of regulatory procedures. P.H. and M.H. received financial support from the Dev2Can LabEx and ITMO Cancer of Aviesan, respectively. C.B. received financial support from the University Claude Bernard Lyon 1 in the framework of the “Etoiles 2022” call for projects. S.Meyer. is supported by a private donation from T. family. The Share4Kids project is supported by INCa and “Enfants Cancers Santé” foundation.