GBM induces bone disruption and OC activity in several anatomical regions of the calvarium that delineate osteogenic edges adjacent to the skull sutures

To investigate whether GBM induces changes in the skull thickness and density, we established two murine glioma models, SB28 and GL261, in immune competent mice. SB28 exhibited a clear mesenchymal subtype signature while GL261 had a hybrid signature of proneural/mesenchymal subtype according to our 12-gene panel for subtyping8 (Extended Data Fig. 1). High-resolution micro-computed tomography (microCT) imaging revealed that both intracranial tumors induced bone resorption, reducing bone density and thickness (Fig. 1a and Extended Data Fig. 2a,b). Bone resorption was most prominent at the osteogenic edges, adjacent to the skull sutures, that are known to be fused at early childhood (at 1 week of age in mice, corresponding to 2 years of age in humans) and ossified during adulthood (at 4 weeks of age in mice, corresponding to 24 years of age in humans)9. These osteogenic edges adjacent to the skull sutures were recently reported to be connected to the underlying dura through ossified vascular channels that allow immune cell trafficking4,5,6,7. Interestingly, the occipital bone, far from the tumor site, had the highest degree of erosion (Fig. 1a,b). Calvarial bone resorption was specific to intracranial tumor and absent in the subcutaneous GBM model. In addition, sham injury, fibroblast injection (mouse embryonic fibroblasts) or stroke lesions had no impact on the skull. In contrast, when the 4T1 breast cancer cells were intracranially injected, calvarial erosion was induced but with a different pattern (Extended Data Fig. 2c), suggesting the importance of the tumor origin. In both intracranial models of GL261 and SB28, both bone density and thickness were significantly reduced at all measured points, except bone density at mid-occipital (Fig. 1a,f,g). Furthermore, SB28 tumors induced greater bone resorption activity than did GL261 tumors, with wider demineralization in the skull above the tumor, and, interestingly, mineralization within the tumor itself (Fig. 1c). Bone resorption started in the midline at early-stage tumor (half-time), being evident by day 10 after tumor injection in SB28 and day 14 in GL261 tumors (Fig. 1d,e and Extended Data Fig. 2a). Taken together, these data suggest that GBM, but not other brain lesions or systemic tumors, induces bone resorption and OC activity, even at anatomical regions distal to the tumor site.

Fig. 1: GBM induces calvarial bone abnormality and OC activity in several anatomical regions of the calvarium delineating osteogenic edges adjacent to the skull sutures.

a, microCT imaging showed that GBM induced bone erosion in the skull mainly at the osteogenic edges adjacent to the skull sutures. Tumor was sterotactically injected in frontal lobe, and the occipital bone, far from tumor site, showed the highest erosion. b, Established criteria to measure skull density and thickness at defined anatomical regions indicated by arrows (bregma, lambda, mid-occipital, sub-occipital). c, SB28 tumor had calcification within the tumor as visualized by microCT. d, Statistical analysis of bone density showed significant reduction in bregma and lambda at early-stage/half-time tumor as compared with sham. SB28 had higher erosion extent than GL261. e, Statistical analysis of bone thickness showed significant reduction of bone thickness in lambda, bregma and sub-occipital in SB28 tumors and to lesser extent in GL261 at early-stage/half-time tumor. f, Statistical analysis of bone density showing significant reduction of bone density in bregma, lambda and sub-occipital in last-stage tumor. g, Statistical analysis of bone thickness showing significant reduction of bone density in lambda and sub-occipital in both SB28 and GL261, in addition to SB28 mid-occipital. Number of animals in statistics: in d,e, early-stage/half-time (4 mice sham; 9 mice GL261; 7 mice SB28); in f,g, late-stage tumor (6 mice sham; 6 mice GL261; 8 mice SB28). Data are presented means ± s.e.m. In d–f, significance was assessed by a one-way ANOVA with Tukey’s multiple comparison test (*P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001). See Supplementary Fig. 2. Scale bars: 1,000 µm. HU, Hounsfield unit.

Source data

GBM induces specific bone morphometry changes in the calvarium and increases skull channels

To perform bone morphometry analysis and to quantify the changes in skull channels, we performed high-resolution (9 µm) microCT at early and late stages of tumor progression. The analysis was focused on the occipital bone to avoid the impact of surgery site or mechanical pressure related to tumor mass, bearing in mind that bone disruption in this region initiated in early-stage tumors, suggesting an early crosstalk between tumor and occipital region (Fig. 2a). The region of interest in the occipital bone is demarked by the lambdoid sutures and external occipital protuberance to enable exact comparison between animals from different groups. The inner table and outer plate represent the cortical bone in the occipital region, while the diploe represents the porous trabecular bone (see the Methods for details).

Fig. 2: GBM induces specific bone morphometry changes in the calvarium and increases skull channel openings.

a, High-resolution microCT (9 µm) of dissected skull from two GBM models (GL261 and SB28), during early-stage (day 8, day 15) and late-stage tumor progression. Aged mice (78–104 weeks) were used for late-stage only. Occipital bone distant from tumor site was visualized with skull channels highlighted and used for bone morphometry analysis. b–e, Statistical analysis of bone morphometry showing significant reduction in BV and Tb.Th in GL261 at days 8 and 15, and in late-stage SB28. Number of animals in statistics: b, day 8 (5 mice sham; 4 mice GL261; 6 mice SB28); c, day 15 (5 mice sham; 4 mice GL261; 4 mice SB28); d, late-stage tumor (5 mice sham; 4 mice GL261; 5 mice SB28); e, aged late-stage tumor (3 mice sham; 3 mice GL261; 3 mice SB28). Data are presented means ± s.e.m. In b–e, significance was assessed by a one-way ANOVA with Tukey’s multiple comparison test (*P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001). f, Statistical analysis quantification of the skull channel diameters at early-stage (day 8 and 15) and late-stage tumor. Three animals per cohort (except GL261 day 15) were used for each of the experimental groups (day 8, day 15 and late-stage). Data are presented means ± s.e.m. One-way ANOVA with Tukey’s multiple comparison test (*P < 0.05; **P < 0.01). g, Statistical analysis of skull channel counts showing significant increase for early-stage and late-stage tumors, GL261 and SB28. Three animals per experimental group (except GL261 day 15 were 4 animals). Data are presented means ± s.e.m. One-way ANOVA with Tukey’s multiple comparison test (*P < 0.05; **P < 0.01). Scale bars: 500 µm. NS, not significant.

Source data

Adult (8–12 weeks) and aged (78–104 weeks) mice were used. Statistical analysis showed that GBM already reduced bone volume (BV) by day 8 and this remained significant at day 15 in GL261, while this last one rebounded in late-stage tumor (Fig. 2b–d). In SB28, BV tended to decrease by day 8 and day 15 and reached a significant level in late-stage. However, the aged mice exhibited significant BV reduction in both late-stage tumor models (Fig. 2e). The bone volume/total volume (BV/TV) also tended to decrease at all stages but reached a significant level in late-stage SB28 only. For trabecular bone thickness (Tb.Th), the reduction pattern was evident in both tumor models, reaching a significant level in both early-stage and late-stage SB28 (Fig. 2b–d), and tended to be significant in SB28 aged mice (Fig. 2e). The trabecular bone number (Tb.N) tended to be reduced at early stages in both tumors, but it rebounded in SB28 by day 15, then continued to reduce in late-stage SB28, although not to a significant level. In GL261, the rebound occurred in late-stage tumor, in both adult and aged mice, without reaching a significant difference (Fig. 2b–f).

To investigate whether these changes were specific to the skull or whether GBM affected other bones similarly, similar analyses were performed on femurs dissected from the same mice included in the skull imaging (Extended Data Fig. 3a). Contrary to the skull, femur BV significantly increased by day 8 in SB28, and in late-stage GL261. In addition, BV/TV and Tb.N were increased by day 15 in GL261 and in late-stage SB28, but only Tb.N was significantly higher at both day 8 and day 15 in GL261 (Extended Data Fig. 3b–d). This suggests that GBM induces differential changes in bone structure and density in the skull compared with femurs. Finally, quantification of skull channels highlighted in Fig. 2a revealed a significant increase in the channel diameter by day 15 in SB28 and in late-stage GL261 (Fig. 2e), while their numbers demonstrated a significant increase in all tumor stages (Fig. 2f). These observations suggest that GBM affects the skull-dura channels in the occipital bone far from the tumor location by increasing their number and diameter.

The calvarial suture lines in patients with GBM had reduced thickness compared with age- and gender-matched controls without brain tumors and such changes were not affected by tumor size

To investigate whether GBM affects calvarial bone in humans, we used a restricted publicly available computed tomography (CT) dataset of The Cancer Imaging Archive (TCIA). Age- and gender-matched controls were mainly patients with strokes, hematoma, aneurysm, epilepsy and subarachnoid hemorrhage from our institute. The TCIA dataset (TCGA-GBM, TCGA-LGG, CPTAC-GBM, RTOG 0625/ACRIN 6677) included 26 patients with GBM with available CT scans as less than 10% of patients with GBM undergo CT imaging (accessed 27 September 2023). Image analysis was performed on these 26 patients with GBM (Extended Data Table 1). In addition, 22 age- and gender-matched controls were used after excluding any cancer history. However, age and gender had no significant correlation with the skull thickness in the control set except that the mid-occipital was thicker in females compared with in males (P < 0.0001). As the human skull is very different from that of a mouse, we included extra points to precisely compare the anatomical areas between subjects with tumor and controls (Fig. 3a,b). We used the crossing point of right and left midorbital axes with coronal sutures, in addition to bregma point, lambda and mid-occipital which we had used in the mice. The reduction in skull thickness was most significant at lambda and mid-occipital (P < 0.0001), followed by bregma, left and right midorbital–coronal axes intersection (P < 0.001, P < 0.01, P < 0.5, respectively) (Fig. 3c). The increased thickness seen at mid-occipital in females versus males disappeared in response to GBM, although the difference at right midorbital–coronal axes intersection was significant in the presence of GBM (P < 0.05) (Extended Data Fig. 4). However, no significant change was evident in bone density (Fig. 3d). Analyzing CT bone thickness with respect to tumor location, no significant differences were noticed, except the lambda in left frontal tumors was thicker than in right frontal tumors (Extended Data Fig. 5). To investigate the impact of tumor size on calvarial bone thickness, T1-weighted gadolinium-enhanced magnetic resonance imaging (MRI) images were volumetrically analyzed (axial, sagittal or coronal axis view) with maximum number of image slices per volume. The contrast-enhanced regions of the tumor were delineated, segmented and correlated with bone thickness. Of 26 patient scans used in CT analysis, 21 were available for MRI analysis and showed no significant correlation between tumor volume and calvarial thickness (Fig. 3e). We thus found that patients with GBM exhibit reduced skull thickness at specific anatomical locations of the calvarial sutures and these changes were not correlated with tumor size.

Fig. 3: Patients with GBM showed a reduced skull thickness compared with age- and gender-matched controls.

CT of patients with GBM and control patients without any tumor history (majority were stroke, subdural and subarachnoid hemorrhage, hematoma and epilepsy) were used for the analysis. a, Representative 3D visualization of the five anatomic areas that were analyzed for bone thickness and density. b, Sagittal CT section pointing to anatomic areas used in the analysis. c, Reduction in skull thickness was highly significant in lambda and mid-occipital, followed by crossing point of midorbital axes and bregma. d, The change in bone density was not significant. In c and d, sample size was 26 patients with GBM and 22 age- and gender-matched controls. Data are presented as means ± s.e.m. Significance was assessed by Mann–Whitney U test (*P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001). e, Linear fit model for tumor volume (cm3) and bone thickness on average measurements of bregma thickness, bregma to midorbital left side, bregma to midorbital right side, lambda and mid-occipital (3 measurements in mm). See also Supplementary Fig. 3. Student’s t-test was used for analyzing volume against bone thickness. A P value < 0.05 was set as statistical significance (*P < 0.05; **P < 0.01; ***P < 0.001). Scale bars: 10 mm.

Source data

GBM induces dynamic changes in the calvarial TRAP/Acp5+ OCs

Bone resorption is a biological process mediated by OCs, a specialized immune cell of macrophage/monocyte origin. To investigate whether OCs in the skull undergo numerical or anatomical distribution changes, we used a fluorescent reporter mouse expressing tdTomato under the endogenous TRAP/Acp5 promotor10. These mice were injected with either GL261 or SB28. Light sheet imaging of the entire skull was used to give an overview of OC distribution during tumor progression. A dense OC cluster was evident in the frontal calvarial bone at day 8, while parietal bone above tumor showed concentrated clusters in late-stage tumor (Fig. 4a, Extended Data Fig. 6 and Supplementary Videos 1–3). To exclude any potential contribution of postmortem artifacts, we used intravital multiphoton imaging. To follow the dynamic changes in OC numbers, we performed a longitudinal study in which we killed mice at an early stage of tumor development (8 d), half-time (15 d) or late-stage, defined by exhibiting intracranial tumor symptoms (3–5 weeks). To avoid the impact of intracranial surgery and tumor growth, only parietal and occipital bones were analyzed. The imaging window was focused in the contralateral side from bregma to lambda. A significant reduction in OC numbers was noted in late-stage brain tumors (Fig. 4b). The distance between neighboring cells was also increased in both GL261 and SB28 at late-stage compared with sham, being wider in SB28 compared with GL261 (Fig. 4b). We classified the OCs into three types according to their size (Fig. 4c). The automatically extracted cell volume divided by three-dimensional (3D)-imaging area revealed a significant reduction in OC cell density in both tumors compared with in sham (Fig. 4d,f). The quantification of cell size and pattern of distribution showed that the small size and big size cells significantly decreased 8 d and 15 d after tumor initiation, then increased in late-stage GL261 tumor. The small size cells became even higher than in the sham control group. The number of big cell clusters decreased after 8 d and 15 d, then rebounded, but this increase was not significant compared with sham (Fig. 4e). In SB28 tumors, a significant decrease was already apparent by day 8 in all three cell types. Although a resurgence was recorded in the intermediate and big size cells by day 15, this increase was not significant compared with day 8 (Fig. 4g). These data suggest that GBM induced dynamic changes in OC cell density and distribution in the skull.

Fig. 4: GBM induces dynamic changes in skull OC number during tumor progression and disturbs OC cluster formation and distribution in the skull.

a, A representative light sheet imaging of SB28 showing segmented OC volume during tumor progression (day 8, day 15 and late-stage tumor). Red signal, OCs expressing TRAP-tdTomato; green, tumor cells. Scale bar: 1,000 µm. b, A representative maximum intensity projection of intravital multiphoton microscopic tiling of skull in TRAP-tdTomato transgenic mice. The skull of the same contralateral side was imaged in sham animal, GL261 and SB28 tumor. Red signal, OCs expressing TRAP-tdTomato; blue, bone tissues; green, tumor cells (crossed to contralateral hemisphere in SB28). Scale bar: 400 µm. Experiment was repeated four times. c, Distance between neighboring OCs was increased in both GL261 and SB28, and to a higher extent in SB28 than GL261 suggesting major loss of OC number in SB28. The distance between every pair of segmented OCs was calculated and the distribution was plotted. The OCs were divided into three groups according to their size: 900–3,000 µm3, small cells; 3,000–2,0000 µm3, big cells; >20,000 µm3, cell cluster. d, Representative images of segmented OCs in the skull classified according to their size. The OCs were divided into three groups according to their size: 900–3,000 µm3, small cells; 3,000–20,000 µm3, big cells; >20,000 µm3, cell cluster. e, GL261 tumors induce a decrease in skull OC numbers at days 8 and 15, and a resurgence in late-stage tumor. A longitudinal quantification in GL261 skull showed that the OC numbers declined significantly at early-stage (day 8) and half-time (day 15) of tumor progression compared with sham (P < 0.0001), then rebounded and increased in numbers significantly in late-stage tumor compared with days 8 and 15 (P < 0.0001), but still less than sham control (P < 0.05). OCs were calculated as areas of OCs per imaging area (mean = 3.72 sham; 2.737 GL261 day 8; 2.503 GL261 day 15; 3.102 GL261 late-stage). f, Quantification of segmented cells in GL261 skull according to their size showed that all three types were significantly decreased at days 8 and 15, compared with sham skull, then rebounded and increased in number during the late-stage of tumor progression to be higher than day 15 in all three subsets (P < 0.0001 for small and big size cells, P < 0.05 for cell cluster), but it was significant only for small size cells compared with sham skull. g, SB28 tumors induce a decrease in skull OC numbers during tumor progression, with a small resurgence at day 15. A longitudinal quantification in SB28 skull showed that the OC numbers declined significantly at early-stage (day 8) and half-time (day 15) of tumor progression compared with sham (P < 0.0001). All three types of tdTomato+ cells were significantly decreased in number compared with sham skull (mean = 3.720 sham; 0.6348 SB28 day 8; 0.5491 SB28 day 15; 0.5068 SB28 late-stage). h, Quantification of segmented cells in SB28 skull according to their size showed that all three types were significantly decreased at days 8 and 15 and late-stage compared with sham skull. A resurgence was recorded in the big size and cell cluster by day 15, but all three types of tdTomato+ cells were significantly decreased in number compared with sham controls (P < 0.0001). Multiphoton images used in c–h were collected from four times repeated experiments, 12 tiling images from: day 8, 3 mice GL261; 3 mice SB28; day 15, 3 mice GL261; 3 mice SB28; late-stage, 6 mice GL261; 4 mice SB28; and 3 mice sham. Data are presented as means ± s.e.m. Significance was assessed using a one-way ANOVA with Tukey’s multiple comparison test (*P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001). See also Supplementary Fig. 6. SHG, second harmonic generation.

Source data

GBM alters the immune landscape of SM

To perform deep molecular profiling of immune cells in the SM and to investigate the impact of GBM on SM and BM, we optimized a protocol to isolate cells from the skull and long bones and performed integrated scRNA-seq analysis. Our scRNA-seq datasets of sham SM (representing SM of healthy animal in homeostatic condition), sham BM, and SM and BM from SB28 and GL261 tumor-bearing mice were integrated. Cells originating from sham SM, GL261 SM, SB28 SM, sham BM, GL261 BM and SB28 BM were visualized in two dimensions by uniform manifold approximation and projection (UMAP) (Fig. 5a,b). CD45+ cells constituted 75.6% of the total cells in sham SM, 91.3% in GL261 SM and 90% in SB28 SM, while for BM, CD45+ made 76.4%, 89.5% and 91.5% of sham, GL261 and SB28, respectively (Supplementary Table 1).

Fig. 5: Global disturbance in myeloid lineage of the SM in response to GBM shown by scRNA-seq and confirmed by FACS.

a,b, 2D UMAPs and stacked bar charts visualizing the main myeloid cell populations within CD45+ cells in the SM and BM of sham control compared with GL261 and SB28 derived from scRNA-seq. a, In homeostatic SM, neutrophil subsets made up the largest myeloid populations, with Mature Neutro. and Transit. Neutro. 1 being the majority, and they almost doubled in response to GBM, whereas Pro. Neutro. increased twice to reach 8% in SB28 and four times in GL261 (16%). Prol. Monocytes, macrophages, Acp5+ macrophage, MEPs and GMPs made up a small subset of 0.2–4% which also increased in response to GBM, except for macrophages in SB28 SM, and DCs stayed almost the same (~2%) in SM. b, In sham BM, similar to SM, neutrophil subsets made the largest cell population, but with some differences including: Pre. Neutro. was the largest population in sham BM, slightly increased in GL261 BM and decreased in SB28 BM, whereas Mature Neutro increased three times, and Transit. Neutro. 1 doubled in GL261, but was only slightly elevated in SB28 BM. MEPs and GMPs as well ST-HSCs and Hemato. Prog. in BM were reduced in response to GBM, opposite to their population in SM, whereas Erythroid. Prog. increased in BM and decreased in SM. c, A dot plot depicting normalized expression of different marker genes used to annotate individual cell clusters within the myeloid lineage. Gene markers are shown on the x axis and cell types on the y axis. d, Flow cytometry analysis confirming that GBM alters the SM and BM homeostasis, pushing it towards myelopoiesis. Spectral flow cytometric analysis of CD45+ population in sham SM, GL261 SM and SB28 SM (top) as well as in sham BM, GL261 BM and SB28 BM (bottom). FACS analysis further recapitulated and validated the increase in percentage of neutrophils in the SM and the concurrent decrease in percentage of the lymphoid compartment. Lin., lineage.

Source data

To define the cell identity of obtained clusters, an immune cell marker panel was created by combining literature searches11,12,13. We identified 13 clusters in the myeloid lineage (Fig. 5a–c) and 12 clusters in the lymphoid lineage (Fig. 6a–c), in addition to hematopoietic stem cells (HSCs) and erythroid progenitors (Supplementary Table 1). The frequency differences and functional properties of major cell types were confirmed by flow cytometry, comparing sham SM versus SB28 SM/GL261 SM and sham BM versus SB28 BM/GL261 BM (Figs. 5d and 6d, Extended Data Figs. 7–9 and Supplementary Figs. 1 and 2). The gating strategy is presented in Supplementary Fig. 2. Surprisingly, a global distortion was observed in SM hematopoietic compartments in both lymphoid and myeloid lineages, shifting the niche toward myeloid predominance (Fig. 5a,b and Supplementary Table 1). To simplify the changes among SM hematopoietic cells in response to GBM, we studied myeloid and lymphoid subsets separately.

Fig. 6: GBM induces distinct DEGs and pathways across myeloid lineages in SM and BM.

a, Venn diagrams showing DEGs and differential pathways induced by two GBM models, GL261 and SB28, in SM as well as BM compared with sham SM and sham BM, respectively. GL261 and SB28 induced common as well as tumor-specific DEGs and pathways in both SM and BM. pyDESEQ2 was used to identify DEGs in all cell clusters. Venn diagram depicts DEGs with more than 1.2-fold change in upregulation/downregulation. For pathway enrichment, GSEA was used to identify enriched GO pathways and only pathways having FDR Q value less than 0.25 were retained for analysis. b, Stacked bar charts showing significantly upregulated and downregulated pathways induced in each cell type by both GL261 and SB28 in SM (left) and BM (right). There was overall upregulation of pathways in SM while BM had downregulation across most cell types. Prol. Monocytes in SM observed most downregulated pathways while Transit. Neutro. 2 in BM observed most upregulated pathways. c, A heatmap depicting NES of select pathways that were downregulated or upregulated across different cell types in SM in response to both GL261 and SB28. Average NES combined from GL261 and SB28 is plotted. d, A heatmap depicting NES of select pathways that were downregulated or upregulated across different cell types in BM in response to both GL261 and SB28. Average NES combined from GL261 and SB28 is plotted. e, 3D PCA plot showing pseudo bulk transcriptomic profile of different samples. Each sample was randomly split into three pseudo replicates and raw counts were summed for each gene to generate bulk RNA profiles with three pseudo replicates for each of the samples. pyDESEQ2 was used to analyze these bulk RNA profiles and normalized counts were then used to perform PCA. Subsequently, the first three principal components were plotted in 3D space. f, Venn diagrams showing DEGs and differential pathways between SM and BM in sham, GL261 and SB28. Sham SM had a significant number of DEGs as well as differential pathways compared with sham BM which remained differential even in the presence of tumor. SB28 induced more DEGs as well as differential pathways compared with GL261. pyDESEQ2 was used to identify DEGs in all cell clusters. Venn diagram depicts DEGs with more than 1.2-fold change in upregulation or downregulation. For pathway enrichment, GSEA was used to identify enriched GO pathways and only pathways having FDR Q value less than 0.25 were retained for analysis. g, A heatmap depicting NES of select pathways that were downregulated or upregulated across different cell types in SM versus BM. NESs for both GL261 and SB28 are plotted unless a pathway is specific only to GL261 or SB28. Gray boxes depict nonenrichment of a pathway in a given cell type and condition. NES, normalized enrichment score.

GBM triggers differential expansion of the myeloid lineage in SM versus BM

The myeloid lineage composed ~41% in sham SM and increased to 90% and 72% in SM of GL261 and SB28, respectively, whereas BM consisted of 67% myeloid cells which increased to 91% and 82% in GL261 and SB28, respectively (Fig. 5a,b).

This myelopoiesis dominancy was driven by neutrophils which composed the largest myeloid population in homeostatic SM (28.2%) and almost doubled in response to GBM (68.2% in GL261 and 52.5% in SB28) (Fig. 5a,b and Supplementary Table 1). Neutrophils comprised four main subpopulations, with Mature Neutro. (defined as Camp-lo/NoNgp-lo/NoLy6g-hiRetnlg-hiCxcr2hi) being the largest cell population (13.3%) which almost doubled in response to GL261 and SB28 (28.7% and 21.8%, respectively). Pre. Neutro. (pre.neutrophils) (Camp-hi, Ngp-hi, Ly6g-loRetnlg-loCxcr2lo) and Transit. Neutro. 1 (transit neutrophils) (Camp-No, Ngp-No and Ly6g-hiRetnlg-loCxcr2lo) expanded twice in SB28, but to a much greater extent than in GL261 (2–4 times compared with sham SM). In sham BM, the Mature Neutro. and Transit. Neutro. 1 represented the biggest populations of CD45+ cells (14–16%) and underwent differential changes in response to tumor type, with the first one being increased to 28% in SB28 and 31% in Gl261, whereas Transit. Neutro. 1 increased to 17.4% in SB28 and to 23.7% in GL261. The Pre. Neutro. were ~6% in sham BM and expanded almost two times in response to GL261, but did not change in SB28 tumors (Fig. 5a–c and Supplementary Table 1). Furthermore, brain tumors led to an increase in Prol. Macrophages (proliferative macrophages), Hemato. Prog. (hematopoietic progenitors) and MEPs (megakaryocyte–erythroid progenitor cells) in the SM, but not in BM. Monocytes, which made up 0.5% in sham BM, increased to 1% in GL261 BM and 1.3% in SB28 BM. While Acp5+ macrophages increased from 1.7% to 3.8% in SB28 SM, their frequency did not change in GL261 SM, but they were highly reduced in the BM in both tumor models (dropped from 0.5% to 0.2%) (Fig. 5e). The ST-HSC (short term-HSC) population (~1.4% of total CD45+ in SM and 2% in BM) was defined as cKit+CD34+Flt3+CD27+Sca-1Lo (Fig. 5a–c and Supplementary Table 1). The high expression of CDK6 in this cluster suggested that these were short-term HSCs (not quiescent) rather than long-term HSCs14. ST-HSC numbers were increased in SB28 SM from 1.4% to 2%, but did not change in GL261 SM. Hemato. Prog. were defined as Ly6c2+Ccr2+Stmn1+Mki67+, which existed in sham SM (0.7%) and increased to 2% in response to SB28, and to a lesser degree in GL261, 0.9%. Both ST-HSCs and Hemato. Prog. were reduced by almost ~50% in BM of both tumors. The Erythro. Prog. (erythroid progenitor) cells expressed Gypa, Alas2 and hemoglobin genes (Hba-a1/2, Hbb-bs and Hbb-bt) and existed in sham SM (6.8%), and decreased to ~1.6% in GL261 SM and 1.2% in SB28 SM, whereas they increased in BM, from 0.2% in sham BM to 3.5% in GL261 BM and 5% in SB28 BM (Fig. 5a–c and Supplementary Table 1).

The differentially expressed gene (DEG) analysis comparing tumor SM with sham SM showed 5,552 shared genes between SM of both tumor types, GL261 and SB28, with 864 genes specific to GL261 and 1,158 specific to SB28 (Fig. 6a). The pathway analysis showed that the myeloid lineage expansion in the SM of both tumors had 544 shared pathways, which were dominated with upregulated genes related to immune response and proliferation in most cell clusters, including ST-HSCs, GMPs (granulocyte-monocyte progenitors), Hemato. Prog., Erythro. Pro., monocytes, DCs (dendritic cells), macrophages, Acp5+ macrophages and Pre. Neutro. (Fig. 6b,c and Supplementary Table 1), whereas the BM had 527 shared pathways dominated with downregulated genes in most cell clusters, except for Pro. Neutro. (pro.neutrophils) and Prol. Macrophages. (Fig. 6b,d). Semantic clustering of enriched gene ontology (GO) terms revealed inflammatory response, proliferation and angiogenesis, DNA replication and migration as commonly occurring GO themes (Supplementary Figs. 4 and 5). The inflammatory and immune response pathways in SM showed upregulation of several pathways, such as positive regulation of macrophages, monocytes, response to type II interferon, positive regulation of neutrophil migration and others (Fig. 6c). Conversely, downregulated pathways related to type II interferon response, IL-1 response, monocyte chemotaxis, B cell proliferation and others, and were evident in BM (Fig. 6d).

Furthermore, to investigate whether SM and BM are different niches and respond differently to the same type of tumor, we compared side-by-side SM versus BM from GL261 and SB28. Principal component analysis (PCA) of pseudobulk scRNA-seq profiles showed that sham SM and BM are different and remain different in response to GL261 while clustering away from sham (Fig. 6e). In contrast, SB28 brought them closer in the PCA space. There were 2,691 DEGs between sham SM versus sham BM, 1,997 of which stayed among DEGs in response to different tumor types, GL261 or SB28, representing 133 differential pathways across all conditions and suggesting SM and BM keep these differences even with existing tumor condition (Fig. 6f). SB28 SM versus BM seemed to have higher DEGs and differential pathways shared with sham SM versus BM than GL261, suggesting GL261 affected SM and BM in a more similar way than SB28 (Fig. 5e–g). Comparison of SB28 SM versus BM showed higher DEGs and pathways specific to the tumor model than GL261, suggesting that SB28 induced a more specific disturbance in the myeloid lineage of SM and BM than did GL261 (Fig. 5e–g). In SB28 SM versus BM, MEPs, Erythro. Prog., macrophages and Acp5+ macrophage clusters seemed to be enriched, with upregulated immune response and activation, in addition to OC activation in Acp5+ macrophages. In contrast, in GL261 SM versus BM, most differential pathways were upregulated in GMPs, macrophages, Pre. Neutro. and Hemato. Prog (Fig. 5g). The functional characterization of the main myeloid populations via flow cytometry showed that macrophages and DCs from BM and SM had higher expression of MHC-II+, suggesting an increase in antigen presentation under tumor condition (both SB28 and GL261), whereas IL-10 was significantly increased in macrophages in the SB28 tumor, but not in GL261 (Extended Data Figs. 7 and 8a,b). Conversely, DCs from this last one showed higher expression of PD-L1 in both BM and SM compared with sham, and only higher IL-10 in SM DCs (Extended Data Fig. 8c,d). Proliferating neutrophils were significantly higher in the SM but not in the BM of GL261, while proliferating monocytes were higher in both BM and SM of GL261 (Extended Data Fig. 9a–d). In the SB28 model, the total number of neutrophils was significantly higher compared with sham SM, but KI67 did not reach a significant level, while both proliferating monocytes and proliferating DCs were significantly higher in SM, not in the BM (Extended Data Fig. 9e,f). These data together demonstrate that SM and BM myeloid compartments are different in the homeostatic condition and have a differential immune response specific to the tumor type.

GBM induces differential alterations in the lymphoid lineage of SM and BM

Opposite to myeloid lineages, a dramatic change in the lymphoid lineage within the SM of GBM animals was evident, with more than 55% reduction in total lymphoid compartments in SB28 and 87% in GL261 SM. In sham BM, the lymphoid lineage composed ~33% of CD45+ cells which decreased to 11% and 4% in BM of SB28 and GL261, respectively. A wide population of various B cell subsets that belong to different developmental stages occupied 49% of total CD45+ cells in sham SM, making B cells the most numerous compartments in the homeostatic SM niche, whereas they made up <29% of CD45+ in BM (Fig. 7a,b and Supplementary Table 1). The SM lymphoid reduction was mainly driven by major loss of almost all B cell subsets (60–94%) (Fig. 7c). This loss was due to a higher extent of GL261 than SB28 BM.

Fig. 7: scRNA-seq and flow cytometry analysis shows a dramatic change in the lymphoid lineage within the SM and BM in response to GBM.

a, 2D UMAPs and stacked bar charts depicting different cell populations and their percentages within the lymphoid subset present in sham SM compared with GL261 and SB28 SM. b, 2D UMAPs and stacked bar charts of corresponding BM derived from scRNA-seq. a,b, In homeostatic SM, B cells were the dominant population within the lymphoid subset, with some differences between SM and BM including higher Pre-B. cells in SM than BM, while this last one had higher Late Pro. B cells which almost disappeared in response to tumor. Early developmental B cell subsets were the most decreased cell types in both SM and BM. While immature B cells were much higher in SM (biggest B cell population) versus BM and reduced almost by 50% in SM and not affect in BM or increased in case of GL261 tumor. Mature B cells also seemed to expand in BM of both tumors while decreasing in SM, and NK cells, NKT cells and T cells were also present in small frequencies. In response to both tumors, T cells, regulatory T cells (Treg cells) and NKT cells increased in the SM while the percentages of NK cells, plasma B cells an