Longitudinal single-cell multiomic atlas of high-risk neuroblastoma reveals chemotherapy-induced tumor microenvironment rewiring

Abstract

High-risk neuroblastoma, a leading cause of pediatric cancer mortality, exhibits substantial intratumoral heterogeneity, contributing to therapeutic resistance. To understand tumor microenvironment evolution during therapy, we longitudinally profiled 22 patients with high-risk neuroblastoma before and after induction chemotherapy using single-nucleus RNA and ATAC sequencing and whole-genome sequencing. This revealed profound shifts in tumor and immune cell subpopulations after therapy and identified enhancer-driven transcriptional regulators of neuroblastoma neoplastic states. Poor outcome correlated with proliferative and metabolically active neoplastic states, whereas more differentiated neuronal-like states predicted better prognosis. Proportions of mesenchymal neoplastic cells increased after therapy and a high proportion correlated with a poorer chemotherapy response. Macrophages significantly expanded towards pro-angiogenic, immunosuppressive and metabolic phenotypes. We identified paracrine signaling networks and validated the HB-EGF-ERBB4 axis between macrophage and neoplastic subsets, which promoted tumor growth through the induction of ERK signaling. These findings collectively reveal intrinsic and extrinsic regulators of therapy response in high-risk neuroblastoma.

Fig. 1. Longitudinal single-cell RNA and ATAC atlas of high-risk neuroblastoma.

a, Overview of the multiomics studies on patient-matched longitudinal neuroblastoma specimens. b,c, UMAPs of the snRNA-seq data (b; n = 372,619 cells) and snATAC-seq data (c; n = 144,366 cells) annotated by major cell type category. d, Dotplot showing the mean expression of marker genes and the percentage of cells expressing them for each annotated cell type. e, Stacked barplots of cell type proportions across the snRNA-seq (left) and snATAC-seq (right) datasets. The cell types are colored as in b. f, Shifts in cell type proportions for each patient between initial diagnosis and post-therapy time points in the snRNA-seq data. Central lines indicate median values, the box edges mark the 25th and 75th percentiles and the whiskers extend 1.5 times the interquartile range. Samples from the same patient between time points are connected by a gray line (n = 22 pairs). Statistical significance was assessed using a one-sided Wilcoxon signed-rank test. Avg., average; DX, diagnosis; PTX, post-induction chemotherapy. Source data

Fig. 2. Therapy-induced neoplastic cell state shifts.

a, UMAP of inferred neoplastic cells from snRNA-seq data after integration and annotation of cell states. b, UMAPs showing the MES, ADRN and cell cycle S-phase signature scores of neoplastic cells. c, Violin plots of MES − ADRN signature score difference (top) and cell cycle S-phase score (bottom) across different cell states. The short black bars represent the median value in each group and the red dashed lines indicate y = 0. d, Heatmap of the top 15 enriched Kyoto Encyclopedia of Genes and Genomes pathways for each neoplastic cell state. The enrichment was conducted based on Fisher’s exact test using enrichR without multiple comparison adjustment. We highlighted in red some labels that are closely related to the naming of each cell state. e, Shifts in cell state proportion between diagnosis and post-therapy samples. A one-sided Wilcoxon signed-rank test was used to calculate statistical significance. Samples from the same patient between time points are connected by a gray line (n = 22 pairs). f, Differences in cell state frequencies between paired post-therapy and diagnostic samples stratified by ALK mutation status. A one-sided Wilcoxon rank-sum test was used to calculate significance (n = 18 ALK wild type (WT) and 4 ALK mutated). In e and f, central lines indicate median values, the box edges mark the 25th and 75th percentiles and the whiskers extend 1.5 times the interquartile range. g, Kaplan–Meier curves of overall survival based on neoplastic cell state using the Sequencing Quality Control project dataset. Patients were stratified into high and low groups based on the median value of the cell state signature score. P values were calculated based on the Cox proportional hazards model and adjusted by age, sex and MYCN amplification status. h, Kaplan–Meier curves of overall survival, with patients grouped into different neoplastic cell states based on maximum cell state signature scores. The numbers of samples per group are indicated. P values were calculated based on the Cox proportional hazards model and adjusted by age, sex and MYCN amplification status. The ADRN-calcium state was chosen as the baseline state. i, Proportions of neoplastic cell states in the initial diagnostic samples. Patients are grouped according to their responses to induction chemotherapy (top) and clinical events (bottom). A one-sided t-test was performed to compare the MES state proportion between two patient groups, as indicated by the vertical dashed line (P = 0.02 (top) and 0.05 (bottom)). cAMP, cyclic AMP; ECM, extracellular matrix. Source data

Fig. 3. Transcriptional regulation of neoplastic cell states.

a, Stacked barplots of neoplastic cell state proportions in the snATAC-seq data at the diagnosis and post-therapy time points. The colors are as in b. b, Dotplot showing the top 15 transcription factors of the transcriptional regulatory network for each neoplastic cell state. The size of each dot represents the fraction of gene targets in the transcriptional regulatory network regulated by each transcription factor. The color represents the chromVAR deviation z score. c, Transcriptional regulatory networks for the ADRN-calcium, ADRN-proliferating and MES cell states. Diamonds represent transcription factors and circles represent target genes. The size of a transcription factor node is proportional to the average difference in the motif chromatin accessibility z score between a given cell state and the rest of the cell states. The size of a target gene node is proportional to the average fold-change of gene expression between a given cell state and the rest of the cell states. The node color indicates the direction of gene expression change between diagnosis and post-therapy samples in each cell state. The edge weight is proportional to the linear regression coefficient for the predicted enhancer–promoter interaction and the fraction of cells that are accessible at the enhancer peak. d, Coverage plot showing normalized chromatin accessibility for neoplastic cell states at the EZH2 locus. The E–P link track represents the predicted enhancer–promoter links colored by the regression coefficient. The transcription factor (TF) motifs present at the enhancer peaks are indicated. Differentially accessible peaks for the ADRN-proliferating state are highlighted. e, Normalized chromatin accessibility of putative EZH2 enhancers across neoplastic cell states. P values were calculated using edgeR on pseudo-bulk data without multiple comparison adjustment (Supplementary Methods). kb, kilobases. Source data

Fig. 4. TAMs in the neuroblastoma microenvironment.

a, UMAP of macrophage subsets (14,866 cells) from snRNA-seq data after integration and annotation. The colors are as in b. b, Dotplot of the average expression of marker genes and the percentage of cells expressing them for each annotated macrophage subset. c, Violin plots of signature scores for immunosuppressive, pro-inflammatory, angiogenesis and phagocytosis macrophages in our macrophage subsets. The short black bars represent the median value in each group and the red dashed lines indicate y = 0. d, Shifts in macrophage subset proportions between diagnosis and post-therapy samples. A one-sided Wilcoxon signed-rank test was used to calculate significance. Samples from the same patient between time points are connected by a gray line (n = 22 pairs). e, Dotplot showing predicted ligand–receptor interactions between neoplastic cells and macrophage subsets. The ligands are from macrophage subsets (top labels) and are listed first in each pair. The receptors are from neoplastic populations (bottom labels) and are listed second in each pair. Both y axes have been used for labeling due to space constraints. Important interactions involving ERBB4 are highlighted in red; the red dashed line indicates recurrent interactions. f, Comparison of the density of HB-EGF protein quantified by CODEX between neighbors of ADRN-like-2 (ERBB4^hi) neuroblasts and neighbors of other neuroblasts (top), as well as between diagnosis and post-therapy samples (bottom). The density was defined as the mean expression of HB-EGF on cells within a 40-µm window, excluding the marker within the center cell. Significance was assessed using a two-sided Wilcoxon rank-sum test. The numbers of cells are: n = 655,573 (other neuroblasts), 17,532 (ERBB4^hi neuroblasts), 221,677 (DX) and 451,428 (PTX). g,h, Representative cell type mask (g) and CODEX images (h). Arrows indicate macrophages (top) and neuroblasts (bottom). i, Distance from each neuroblast cell to the nearest CD163⁺CDCD68^hi macrophage across samples, stratified by neuroblast population. Numbers of cells in each group (from left to right): n = 396,277, 17,532, 241,995 and 17,301. Significance was assessed using a two-sided Wilcoxon rank-sum test. Outliers were truncated for visualization purposes. In d, f and i, central lines indicate median values, the box edges mark the 25th and 75th percentiles and the whiskers extend 1.5 times the interquartile range. Source data

Fig. 5. Spatial transcriptomic analysis of murine neuroblastoma.

a,b, UMAP projection of Xenium transcriptomic data (993,070 cells) annotated by cell cluster (a) and major cell type (b). c, Dotplot showing the normalized expression levels of marker genes and the percentages of cells expressing them for each annotated cell type. Each row represents a cell cluster, with the average gene expression for each cluster normalized to a range between 0 and 1 across clusters. d,e, Dotplots showing the normalized signature scores for each neoplastic cell state (d) and macrophage subset (e). Each row represents a predicted cell subpopulation, with the average signature score for each cell state normalized to a range between 0 and 1 across cell types. f,g, Barplots displaying the proportions of projected neoplastic cell states (f) and macrophage subsets (g) for treated mice and controls. h, Spatial co-localization analysis of ligand–receptor interactions predicted between neoplastic cells and macrophage subsets based on snRNA-seq data. The Hbegf–Erbb4 interaction is highlighted in red. i, Representative image illustrating the spatial co-localization of Erbb4⁺ neuroblasts and Hbegf⁺ macrophages. The dots represent individual transcripts. j, Comparison of spatial distances between Hbegf⁺ macrophages and Erbb4⁺ neuroblasts versus other neuroblasts. Significance was assessed by two-sided (left) and one-sided (right) Wilcoxon rank-sum test. Source data

Fig. 6. Macrophage-secreted HB-EGF activates ERK signaling and promotes proliferation.

a, Representative western blot of cell-surface HB-EGF (pro-HB-EGF) from THP-1 macrophages in monoculture and co-culture with neuroblastoma cell lines. b, Quantification of the results from a across all replicates, normalized to β-actin (n = 3 for COG-N-297 and n = 4 for the others). Central lines indicate median values, the box edges mark the 25th and 75th percentiles and the whiskers extend 1.5 times the interquartile range. c, Ligand concentrations in the media of THP-1 macrophage monoculture and co-culture with neuroblastoma cell lines, measured by ELISA (n = 3). d, Phosphorylated ERBB4 (pERBB4) levels in neuroblastoma cells after co-culture with THP-1 macrophages, measured by ELISA (n = 4 for CHLA15 and n = 3 for the others). Significance in b and d was calculated using a Welch’s two-sided t-test. e, Representative western blots showing ERK activation in neuroblastoma cells with and without macrophage co-culture and the HB-EGF inhibitor CRM197. f, Quantification of the results from e across all replicates, normalized to total ERK (n = 5, 3 and 4, from left to right). g, Representative western blots showing AKT phosphorylation in neuroblastoma cells with and without macrophage co-culture and the HB-EGF inhibitor CRM197. h, Quantification of the results from g across replicates, normalized to total AKT (n = 3 for the treatment group and n = 5 for the others). i, Representative images (left) and quantification (right) of the area of colony formation for neuroblastoma cells co-cultured with THP-1 macrophages with or without treatment with CRM197 (from left to right: n = 3, 3 and 2 per condition). j, Representative images (left) and quantification (right) of the area of colony formation for neuroblastoma cells co-cultured with THP-1 macrophages with or without treatment with the pan-ERBB inhibitor afatinib (n = 4). The experiments in i and j were repeated two to four times with two to four samples per condition and the averaged values across samples are shown. Significance was calculated using a one-sided paired t-test (f, i and j) or two-sided paired t-test (h). The error bars in a–c, f and h–j represent means ± s.d. k, Neuroblastoma cells stimulate HB-EGF secretion from THP-1-derived macrophages, which reciprocally induce the phosphorylation of ERBB4 on neuroblastoma cells. Activation of ERBB4 stimulates proliferation via the ERK pathway. Panel k was created with BioRender.com. Source data

Fig. 7. Transcriptomic analysis of mono- and co-cultured macrophages and neuroblasts.

a, Left, UMAP plots showing cells from monoculture and co-culture with THP-1 macrophages, colored by the normalized expression levels of PHOX2B and CD68 for the CHLA15 cell line. Right, barplots depicting the proportions of projected neoplastic cell states across experimental conditions. b, UMAP and neoplastic cell state inference of co-culture experiments with the CHLA20 cell line as in a. The P values in a and b were calculated using a two-sided proportion test, comparing either co-culture versus monoculture or co-culture with treatment versus co-culture without treatment. No multiple comparison adjustment was made. c, Differential pathway analysis of CHLA15 cells using the hallmark pathways from the Molecular Signatures Database. Pathways in the co-culture condition were compared versus the monoculture condition and pathways in the co-culture with treatment condition were compared versus the co-culture without treatment condition. d, Pathway analysis of co-cultured CHLA20 cells as in c. The P values in c and d were calculated based on Fisher’s exact test using enrichR and corrected using the Benjamini–Hochberg procedure. *P < 0.1; **P < 0.01; ***P < 0.001. FDR, false discovery rate; IFN, interferon; IL-2, interleukin 2; mTORC1, mammalian target of rapamycin complex 1; TNF, tumor necrosis factor. Source data

Extended Data Fig. 1. Quality assessment of snRNA-seq and snATAC-seq data.

a) Bar plots showing the number of cells sequenced (top) and the number of genes detected per cell (bottom) in each snRNA-seq sample after quality control filtering. b) Bar plots showing the number of cells sequenced (top) and the number of unique chromatin fragments detected per cell (bottom) in each snATAC-seq sample after quality control filtering. c) Dot plot showing the average gene activity and percentage of cells with gene activity for cell type marker genes in the snATAC-seq data. d) Cell type composition of each sample based on snRNA-seq data. e) Scatter plot of cell type proportions in snRNA-seq and snATAC-seq data. Each dot represents a sample. Pearson correlation coefficient was indicated in the plot. f) Shifts in cell type proportion between diagnosis and post-therapy samples in snATAC-seq data. Samples from the same patient between time points were connected by a grey line (n = 22 pairs). Statistical significance was assessed using a one-sided Wilcoxon signed-rank test. Source data

Extended Data Fig. 2. Malignant cell calling by using whole genome sequencing (WGS) and snRNA-seq data.

a) Genome segmentation for two representative patients using whole genome sequencing (WGS). Black points are equal sized bins. Red line shows the Hidden Markov Model (HMM) copy number states. b) Heatmaps of copy number alteration profiles from samples in (a) using snRNA-seq data derived from Clonalscope. c) UMAP representation of cells from snRNA-seq data from samples in (a) colored by predicted malignant or non-malignant prediction. d) Receiver Operating Characteristics (ROC) curve for supervised classification by artificial neural networks. Value of the Area Under the Curve (AUC) is shown inside the plot. e) Correlation between histologic estimation of neoplastic cell percentage by pathologists (x-axis) and predicted malignant cell proportions (y-axis). Each dot represents an individual sample. Linear regression line is shown in red. Pearson correlation coefficient (r) is shown. Significance was assessed using a two-sided test for association between paired samples. f) Heatmaps of inferCNV results of representative patient samples after inference of putative normal cells (top panels) and neoplastic cells (bottom panels). Known recurrent neuroblastoma copy number variations (CNVs) are highlighted with a magenta box. Source data

Extended Data Fig. 3. Malignant cell annotation and survival analysis.

a) Barplot of fraction of malignant cells in each cell type. b) Heatmap showing normalized expression of top 15 up-regulated genes in each neoplastic cell state, down sampled to 100 cells per state for illustration purpose. c) Heatmap of top 15 enriched Gene Ontology Biological Process terms for each neoplastic cell state. d) Stacked barplots of neoplastic cell state proportions between diagnosis and post-therapy timepoints. e) Malignant cell state proportions in snRNA-seq data of each patient across timepoints, ALK mutation, and MYCN amplification status. f) The changes in cell state proportions between diagnosis and post-therapy samples stratified by MYCN amplification status. (n = 11 MYCN-amplified, 11 MYCN-non-amplified). The one-sided Wilcoxon rank-sum test was used to calculate significance. Source data

Extended Data Fig. 4. Projection of neoplastic cells onto normal adrenal medullary developmental trajectories.

a) Stacked barplots of projected cell type proportions across initial diagnosis and post-therapy samples. b) Shifts in projected cell type proportions for each patient between initial diagnosis and post-therapy time points. Samples from the same patient at different time points were connected by a grey line (n = 22 pairs). Significance was assessed using a one-sided Wilcoxon signed-rank test. SCPs, Schwann cell precursors. c) Heatmap showing projected cell type fractions for each neoplastic cell state. d) Heatmap showing the enrichment of projected cell types for neoplastic cell states. Pearson residuals were first calculated based on the contingency table to test the independence of the projected cell type from the neoplastic cell states. The p-value for each cell in the heatmap was calculated based on the assumption that the Pearson residual follows a normal distribution. The enrichment score was defined as the -log10(p-value) multiplied by the sign of the Pearson residual. e) The reference UMAP of normal adrenal medullary developmental trajectories. f) Visualization of neoplastic cells on the reference UMAP stratified by cell state. Source data

Extended Data Fig. 5. Survival analysis of neoplastic cell state signatures.

a) Kaplan-Meier curves of event free survival based on neoplastic cell states using the Sequencing Quality Control project (SEQC) neuroblastoma RNA-seq dataset. Patients were stratified into high and low groups based on the median value of the cell state signature score. b) Kaplan-Meier curves of event free survival based on the maximum cell state signature scores in the SEQC dataset. Patients were grouped based on the cell state with the maximum signature score. The number of samples and p-value for each group are indicated in the parentheses. The ADRN-Calcium state was chosen as the baseline. c-f) Kaplan-Meier curves of overall survival (c-d) and event free survival (e-f) using the Cangelosi et al. RNA-seq dataset based on patient stratification (c, e) and maximum cell state signature score (d, f) as in (a) and (b) respectively. P-values in a-f were calculated using the Cox proportional hazards model and adjusted by age, sex and MYCN amplification status with no multiple comparison adjustment. g-h) Comparison of proportions of neoplastic cell states based on the deconvolution of SEQC (g) and Cangelosi et al. (h) datasets between MYCN-amplified and MYCN non-amplified samples. MYCN-amplified, non-amplified: n = 96 and 401 (g); n = 84 and 333 (h). i-j) Comparison of proportions of neoplastic cell states based on the deconvolution of SEQC (i) and Cangelosi et al. (j) datasets between disease stages. Stage1, 2, 3, 4s, 4: n = 120, 78, 62, 52, 181 (i); n = 59, 70, 56, 49, 184 (j).The one-sided Wilcoxon rank-sum test was used to calculate statistical significance (g-j). Source data

Extended Data Fig. 6. Neoplastic cell states in snATAC-seq data.

a) UMAP of integrated snATAC-seq data (93,261 cells) annotated by neoplastic cell clusters. b) Patient proportions in each of the snATAC-seq clusters. c-e) Coverage plots showing normalized chromatic accessibility for each neoplastic cell state at the YAP1 (c), PHOX2B (d), and MYCN (e) loci. f) Dot plot of the average gene activity and percentage of accessible cells of the ADRN and MES genes for each predicted neoplastic population. g) Coverage plots for differentially expressed genes and associated chromatin peaks in each neoplastic cell state: SMC4 (ADRN-Proliferating), APOE (MES), ZIC2 (ADRN-Baseline), KCNQ3 (ADRN-Calcium), DBH (ADRN-Dopaminergic), RPL32 (Interm-OxPhos) loci. The E-P link track represents the predicted enhancer-promoter links colored by the regression coefficient, and the TF motifs present at the enhancer peaks are indicated. The differential accessible peaks (DAPs) for the corresponding cell state are highlighted in blue and yellow for promoter and enhancer peaks, respectively. Source data

Extended Data Fig. 7. Epigenetic regulation of neoplastic cell states.

a) Shifts in cell state proportions in snATAC-seq between diagnosis and post-therapy samples. Samples from the same patient at different time points were connected by a grey line (n = 13 pairs). Significance was assessed using a one-sided Wilcoxon signed-rank test. b) Heatmap of top 20 transcription factors (TFs) with differential motif chromatin accessibility in each cell state. c) Transcriptional regulatory networks (TRNs) for ADRN-Dopaminergic, ADRN-Baseline and Interm-OxPhos cell states. Diamond represents TF and circle represents target gene. The size of a TF node is proportional to the average difference in motif chromatin accessibility z-score between a given cell state and the rest of cell states. The size of a target gene node is proportional to the average fold change of gene expression between a given cell state and the rest of cell states. Node color indicates direction of gene expression change between diagnosis and post-therapy samples in each cell state. The edge weight is proportional to the linear regression coefficient for the predicted enhancer-promoter interaction and the fraction of cells that are accessible at the enhancer peak. d) Fractions of state-specific genes in each TRN that were upregulated, downregulated, or non-significantly changed post-therapy. e) Coverage plot showing normalized chromatin accessibility and gene expression for the MES state-specific gene NECTIN2. The E-P link track represents the predicted enhancer-promoter links colored by the regression coefficient, and the TF motifs present at the enhancer peaks are indicated. The differential accessible peaks (DAPs) for the MES state are highlighted in blue and yellow for promoter and enhancer peaks, respectively. Source data

Extended Data Fig. 8. Therapy-induced shifts in macrophage populations.

a) UMAP plots of macrophages from snRNA-seq data clustered at different resolutions. b) Stacked bar plots of macrophage subset proportions in snRNA-seq data at diagnosis and post-therapy timepoints. c) Representative signature genes for immunosuppressive (A2M, CD84, HLA-E, SPP1), angiogenic (VCAN, VEGFA, VAV2, CXCR4), phagocytic (MERTK, MRC1) and proinflammatory (IL18, CD80) macrophages. d) Proportions of THY1⁺ and proliferating macrophages in diagnosis and post-therapy samples. Significance was assessed using a one-sided Wilcoxon signed-rank test (n = 22 pairs). e) Macrophage subset proportions in snRNA-seq data of each patient across timepoints, ALK mutation, and MYCN amplification status. f-g) Proportions of macrophage subsets in initial diagnostic samples. Patients were grouped according to their responses to induction chemotherapy (f) and adverse clinical events (g). A one-sided t-test was performed to compare the proportion of IL18⁺ macrophages between two patient groups, as indicated by the vertical dashed line; p = 0.005 (f) and 0.02 (g). Source data

Extended Data Fig. 9. Clinical correlation and transcriptomic validation of macrophage populations.

a) Comparison of macrophage subset proportions between MYCN amplified and non-amplified neuroblastoma patients using CIBERSORTx deconvolution of bulk RNA-seq data (SECQ and Cangelosi et al. datasets). MYCN-amplified, non-amplified: n = 96 and 401 (SEQC); n = 84 and 333 (Cangelosi et al.). b) Comparison of macrophage subset proportions between disease stages in bulk RNA-seq data as in (a). Stage 1, 2, 3, 4s, 4: n = 120, 78, 62, 52, 181 (SEQC); n = 59, 70, 56, 49, 184 (Cangelosi et al.). Significance was assessed in (a-b) using a one-sided Wilcoxon signed-rank test. c) Dot plot showing predicted ligand-receptor interactions between neoplastic cell states and macrophage subsets. Ligands were from neoplastic cells and receptors were from macrophage subsets. d) Normalized gene expression of ERBB4 ligands (HB-EGF, EREG, TGFA) in different macrophage subsets from DX and PTX samples, upregulated post-therapy using logistic regression with no multiple comparison adjustment (Methods) e) Cell-cell interaction analysis using CellChat and LIANA. Left, HB-EGF-ERBB4-mediated interactions among neoplastic cell states and non-neoplastic cells by CellChat. Macrophage subsets are colored in blue and neoplastic populations are colored in red. Right, Top 10 ligand-receptor pairs identified between VCAN+ macrophages (ligand source) and neoplastic cell states (receptor source) by LIANA using three representative signaling pathway databases. HB-EGF-ERBB4 interaction is highlighted in red. Source data

Extended Data Fig. 10. Annotation and protein expression in CODEX spatial imaging.

a) List of CODEX panel markers for different cell types. b) UMAP visualization of the ~841,000 cells in CODEX data. c) Stacked barplot of cell type proportions across samples in the CODEX dataset. d) Heatmap of normalized protein expression across cell types in CODEX data. e) Representative images of CODEX immunofluorescence, cell phenotype mask (CPM), and hematoxylin and eosin (H&E) stain for selected cell types. f) Hematoxylin and eosin (H&E) image (left) and cell phenotype mask (right) of representative images for each CODEX sample. g) Comparison of the density of TGFA protein quantified by CODEX between neighbors of ADRN-like-2 (ERBB4^hi) neuroblasts and the neighbors of other neuroblasts (top), and between diagnosis and post-therapy samples (bottom). The density was defined as the mean expression of TGFA on cells within a 40 µm square, excluding the marker within the center cell. Significance was assessed using a two-sided Wilcoxon rank-sum test. The numbers of cells are: n=655,573 and 17,532 (top); n = 221,677 and 451,428 (bottom). h) Distance from each neuroblast cell to the nearest CD163⁺CD206^hi macrophage across samples, stratified by neuroblast population. Significance was assessed using a two-sided Wilcoxon rank-sum test. Outliers were truncated for visualization purposes. Source data