Single-cell transcriptomics and epigenomics unravel the role of monocytes in neuroblastoma bone marrow metastasis

Abstract

Metastasis is the major cause of cancer-related deaths. Neuroblastoma (NB), a childhood tumor has been molecularly defined at the primary cancer site, however, the bone marrow (BM) as the metastatic niche of NB is poorly characterized. Here we perform single-cell transcriptomic and epigenomic profiling of BM aspirates from 11 subjects spanning three major NB subtypes and compare these to five age-matched and metastasis-free BM, followed by in-depth single cell analyses of tissue diversity and cell-cell interactions, as well as functional validation. We show that cellular plasticity of NB tumor cells is conserved upon metastasis and tumor cell type composition is NB subtype-dependent. NB cells signal to the BM microenvironment, rewiring via macrophage mgration inhibitory factor and midkine signaling specifically monocytes, which exhibit M1 and M2 features, are marked by activation of pro- and anti-inflammatory programs, and express tumor-promoting factors, reminiscent of tumor-associated macrophages. The interactions and pathways characterized in our study provide the basis for therapeutic approaches that target tumor-to-microenvironment interactions.

Fig. 1. Integrated analysis of NB-infiltrated BM by single-cell RNA sequencing.

a Overview of our experimental and analysis approach. DTC disseminated tumor cells, ME microenvironment b UMAP of the overall scRNA-seq dataset (n = 80,789 cells) after removal of patient effects. Numbers denote clusters, while colors distinguish cell types. c Relative abundances of inferred cell types in each cluster. The annotation to the right indicates the final cell type classification assigned to each cluster. d Expression of canonical cell type marker genes in clusters shown in b. e Comparison of copy number variation regions (colored areas) predicted from scRNA-seq data (bottom rows) to regions deduced from log R ratios (dark lines) in SNP array data (top rows). One patient from each respective NB subtype is shown. Source data are provided as a Source Data file.

Fig. 2. Characterization of the metastatic NB cells in the BM.

a Expression of NB marker genes in NB cells across patients (left) and in other cell types (right). b Proportion of NB cells assigned to each adrenal medulla cell type in primary and metastatic tumor cells. NB subtypes are grouped into MNA and non-MNA (A + S) tumors. c Correlation analysis of pseudobulk gene expression demarcates the M tumors from A and S tumors. Source data are provided as a Source Data file.

Fig. 3. Cell-cell communication analysis of the NB tumors with the BM microenvironment.

a Number of ligand-receptor interactions of metastatic NB cells with the BM microenvironment. b Relative contribution of each ligand-receptor pair to the overall communication network. c Ligand-receptor pairs and their relative communication score in each cell type interaction. d, e Relative importance of each cell type based on the computed four network centrality measures of MIF and MK signaling networks, receptively. f, g RNA expression levels of ligands and receptors involved in the MIF and MK pathway, respectively, in individual patients. h–i Label-free protein quantitation of MIF and MK in ganglioneuroma (GNM, n = 6 biologically independent samples), primary NB tumors (NB-TU, n = 3 biologically independent samples), human mesenchymal stem cells (hMSC, n = 3 independent experiments), and NB cell lines (n = 3 biologically independent samples), and of CD44, CD74, LRP1, and NCL in primary monocytes stimulated with LPS (n = 3 independent experiments) and corresponding controls (n = 3 independent experiments). Data were subjected to two-tailed unpaired Student’s t test for MIF (GNM vs. NB-TU, p = 0.46 and hMSC vs. NB-CL, p = 0.0025) and MK (GNM vs. NB-TU, p = 0.0005 and hMSC vs. NB-CL, p = 0.0022) or two-tailed paired Student’s t test (CD44, p = 0.022; CD74, p = 0.027; LRP1, p = 0.017; NCL, p = 0.23). j Secreted levels of MIF and MK as determined by ELISA in cell culture supernatants of CLB-Ma (n = 4 independent experiments for PBMCs, n = 5 independent experiments for all the other conditions) and CHLA90 (n = 2 independent experiments for PBMCs, n = 3 independent experiments for all the other conditions). k Mean fluorescence intensity (MFI) of cell surface proteins: CD44, CD74, CXCR4, LRP1, and intracellular receptor: NCL, as determined by flow cytometry in NB cells, co-cultured either directly with PBMCs or through a trans-well for 3 days. Cells were gated for GD2+ CD45−, CD14+ CD16+, and CD14− CD16+ populations (n = 3 independent experiments). Data (j–k) were subjected to one-way ANOVA and corrected using Dunnett’s post hoc test for multiple comparisons. Asterisks indicate statically significant changes compared to control (PBMC or CLB-Ma): *p < 0.05, **p < 0.01, ***p < 0.001. Data are presented as mean ± standard error of the mean (h–k). Source data are provided as a Source Data file.

Fig. 4. Cellular composition of NB-infiltrated BM.

a Cell type abundance in NB subtypes (n = 4, 2 and 5 patients for MYCN amplified, ATRX^mut, and sporadic, respectively) compared to control (n = 5 patients). b Correlation analysis of gene expression changes between NB patients and controls for BM microenvironment cell types in NB subtypes. c Gene set enrichment analysis using MSigDB Hallmark 2020 database with genes sorted according to their log-fold change of expression relative to the control group. p-values are based on an adaptive multi-level split Monte-Carlo scheme as implemented in the R package fgsea. d Expression of exemplary genes in our patient cohort for the top two enriched/depleted pathways identified in c. e Label-free protein quantitation of targets identified in d in U937 cells (n = 3 independent experiments). f, g Protein expression in monocytes and macrophages derived from multiplex images of neuroblastoma BM samples with no (control, C, n = 3 patients) and high (>200 cells, n = 5 patients) NB cell infiltration. Scale bar, 45 µm. Wilcoxon–Mann–Whitney with FDR-corrected p-values: ns not significant, **p ≤ 0.01, ***p ≤ 0.001. h Secreted levels of M1 and M2 markers in co-cultured NB cells with PBMCs (n = 2 independent experiments). i Representative FACS plots of NB cells and CD45+ populations in co-cultured NB cells and PBMCs, and controls, along with PBMCs induced to acquire an M1 or M2 phenotype at day 3. j Percentage (%) of cells expressing the M1 marker, CD86 or the M2 marker, CD163 in CD45+, CD14+ CD16+, and CD14dim CD16dim populations as well as % and MFI of MHC class I and II markers (n = 3 independent experiments). Dots show the raw data (a, g) and boxes display the median value and 25 and 75% quartiles; the whiskers are extended to the most extreme value inside the 1.5-fold interquartile range. Data were subjected to two-tailed paired Student’s t test (e) or one-way ANOVA (j, compared to control: PBMC or CLB-Ma), and corrected using Dunnett’s post hoc test for multiple comparisons. Asterisks indicate statically significant changes, *p < 0.05, **p < 0.01, ***p  <  0.001. Data are presented as mean ± standard error of the mean (j). Source data are provided as a Source Data file.

Fig. 5. Single cell regulatory landscape of NB-infiltrated BM.

a UMAP projection of scATAC-seq profiles of NB and BM cell types. Dots indicate individual cells, whereas colors indicate cell type identity. b Scatterplot of scRNA-seq and scATAC-seq log fold changes for genes in TNFα pathways and E2F targets. c Representative sequencing tracks for the NFKB1 and KDM2B loci show distinct pseudo-bulk ATAC-seq peaks in the NB subtypes compared to controls. The ATAC-seq data have been normalized with Signac and the scale on the y-axis was chosen for optimal visualization of peaks for each sample. d Dot plot depicting odds ratio and adjusted p-value (calculated with hyper-geometric statistical test implemented in HOMER) of motif enrichment in promoter and distal regions between NB subtypes (n = 4, 2 and 5 patients for MYCN amplified, ATRX^mut, and sporadic, respectively) compared to control samples (n = 5 patients). e Differences in active binding of individual TFs between NB subtypes inferred with footprinting analysis, where the confidence interval represents variations between the patients of the group. Source data are provided as a Source Data file.

Fig. 6. Schematic illustration of the main findings and conclusions.

The figure depicts interactions between NB and myeloid cells in the bone marrow compartment, which are mediated through the MIF (Macrophage Migration Inhibitory Factor) and MK (Midkine) pathways. M, MYCN amplified; A, ATRX^mut; S, sporadic.