Neuroblastoma is a common pediatric cancer, where preclinical studies suggest that a mesenchymal-like gene expression program contributes to chemotherapy resistance. However, clinical outcomes remain poor, implying we need a better understanding of the relationship between patient tumor heterogeneity and preclinical models. Here, we generated single-cell RNA-seq maps of neuroblastoma cell lines, patient-derived xenograft models (PDX), and a genetically engineered mouse model (GEMM). We developed an unsupervised machine learning approach ('automatic consensus nonnegative matrix factorization' (acNMF)) to compare the gene expression programs found in preclinical models to a large cohort of patient tumors. We confirmed a weakly expressed, mesenchymal-like program in otherwise adrenergic cancer cells in some pre-treated high-risk patient tumors, but this appears distinct from the presumptive drug-resistance mesenchymal programs evident in cell lines. Surprisingly however, this weak-mesenchymal-like program was maintained in PDX and could be chemotherapy-induced in our GEMM after only 24 hours, suggesting an uncharacterized therapy-escape mechanism. Collectively, our findings improve the understanding of how neuroblastoma patient tumor heterogeneity is reflected in preclinical models, provides a comprehensive integrated resource, and a generalizable set of computational methodologies for the joint analysis of clinical and pre-clinical single-cell RNA-seq datasets.

ABSTRACT Spatial transcriptomics technologies have recently emerged as a powerful tool for measuring spatially resolved gene expression directly in tissues sections, revealing cell types and their dysfunction in unprecedented detail. However, spatial transcriptomics technologies are limited in their ability to separate transcriptionally similar cell types and can suffer further difficulties identifying cell types in slide regions where transcript capture is low. Here, we describe a conceptually novel methodology that can computationally integrate spatial transcriptomics data with cell-type-informative paired tissue images, obtained from, for example, the reverse side of the same tissue section, to improve inferences of tissue cell type composition in spatial transcriptomics data. The underlying statistical approach is generalizable to any spatial transcriptomics protocol where informative paired tissue images can be obtained. We demonstrate a use case leveraging cell-type-specific immunofluorescence markers obtained on mouse brain tissue sections and a use case for leveraging the output of AI annotated H&E tissue images, which we used to markedly improve the identification of clinically relevant immune cell infiltration in breast cancer tissue. Thus, combining spatial transcriptomics data with paired tissue images has the potential to improve the identification of cell types and hence to improve the applications of spatial transcriptomics that rely on accurate cell type identification.

Retinoic acid (RA) is a standard-of-care neuroblastoma drug thought to be effective by inducing differentiation. Curiously, RA has little effect on primary human tumors during upfront treatment but can eliminate neuroblastoma cells from the bone marrow during post-chemo consolidation therapy-a discrepancy that has never been explained. To investigate this, we treated a large cohort of neuroblastoma cell lines with RA and observed that the most RA-sensitive cells predominantly undergo apoptosis or senescence, rather than differentiation. We conducted genome-wide CRISPR knockout screens under RA treatment, which identified BMP signaling as controlling the apoptosis/senescence vs differentiation cell fate decision and determining RA's overall potency. We then discovered that BMP signaling activity is markedly higher in neuroblastoma patient samples at bone marrow metastatic sites, providing a plausible explanation for RA's ability to clear neuroblastoma cells specifically from the bone marrow, seemingly mimicking interactions between BMP and RA during normal development.