To investigate the tumor heterogeneity of FNB, we sequenced a total of six samples ( Figure 1A ) including two cases of FNB(F1, F2). The six samples were divided into two groups. The P1 group included F1 and two cases of stage I/II sporadic neuroblastoma (A8 and A53), while the P2 group included F2 and two cases of stage IVs sporadic neuroblastoma (A5 and A24) ( Figure 1C ). All six samples were examined histopathologically, including one case of ganglioneuroblastoma and five cases of neuroblastoma, and N-MYC was not amplified in all cases. After stringent quality filtering, we obtained a total of 35,369 cells from the six samples, and two cases of familial neuroblastoma were identified with 9482 and 8912 cells. A total of ten cell types were identified in the six samples, including neuroendocrine cells(NEs), Schwann cells, T cells, B cells, mononuclear phagocytes(MPs), fibroblasts, endothelial cells(ECs),mural cells, plasmacytoid dendritic cells(pDCs), and hepatocytes ( Figure 1B ). The cell types shared by the six samples included neuroendocrine cells (NEs), T cells, ECs, and MPs. Uniform manifold approximation and projection (UMAP) was used to visualize cell clusters and violin plots were used to show the canonical markers in each cell type ( Figure 1E ). Consistent with the literature, the T cells were distinguished by high expression of CD2,CD3D,TRAC, and TRBC2 (20, 21), while the B cells were characterized by high expression of MS4A1,CD79A, and CD79B (20, 22). The pDCs specifically expressed IL3RA, CLEC4C, and GZMB (23, 24); the MPs had high expression of CD14, VCAN, CD1C, and C1QA (20, 25, 26); and the mural cells had high expression of RGS5,ACTA2,PDGFRB, and NOTCH3 (27, 28).We also identified fibroblasts by high expression of DCN,COL1A2, and COL1A1 (29, 30), ECs by high expression of CDH5,PECAM1,VWF, and CLDN5 (31), and NEs by high expression of CHGB,CHGA,NPY, and SCG2 (32, 33).The schwann cells were characterized by the expression of S100B,CRYAB,and MPZ (34, 35), the hepatocytes were characterized by the expression of ALB,APOA1,HP (26, 36).

We further compared the proportion of cell types in FNB versus sporadic neuroblastoma, and among the six shared cell types, we found that FNB had fewer NE cells and more immune cells in terms of number share ( Figure 1D ), which also represented a more complex tumor immune microenvironment and closer immune cell-to-cell communication in FNB compared to sporadic neuroblastoma. In addition, we identified partial hepatocytes in F2, associated with tissue samples that could carry partial liver tissue.

To classify the cells into malignant and non-malignant cells, we used the inferred CNV algorithm ( Figure 2A ) to demonstrate the clonal structure of the cells. The results showed that NEs have more CNVs than other cell types. By comparison, we observed significant chromosome 17p gain in two cases of FNB and only chromosome 19p gain in F1 ( Figure 2B ). Similarly, we also focused on mesenchymal cells including fibroblasts and endothelial cells, and the results showed that they have few CNVs, so we considered NEs to be malignant cells.

To study the genomic features of NEs in FNB, we performed a Hotspot analysis ( Figure 2C ) of NEs (37). We obtained a total of 15 modules, and by Jaccard similarity analysis ( Figure 2D ), we performed one-to-one correspondence between gene modules and tumor samples. We found that A8 and A53, A5 and A24 had a similar module expression, but no more consistent module expression was observed between F1 and F2, where F1 mainly expressed module1 (including RPS2, RPL18A, RPS29, RPL39), module2 (including RPL34, RPS27, MTND1P23, RACK1), module6 (including MT-CO2, MTATP6P1, RPL35, MT-CO3), F2 mainly expressed in module5 (including MEG3, JUN, HSPA1A, FOS), and module13 (including CALM2, TMSB10, TUBA1A, TMSB4X). Then, we performed survival analysis of the 15 modules ( Figure 2E ) and we found that the module 13 in highly expressed F2 was associated with better survival. The module13 included genes currently known to be associated with poor prognosis ( Figure 2F ) including MIF (38), DDX5 (39), and also genes associated with good prognosis including UCHL1 (40), and STMN4 (41).

Next, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the module 13, which showed that the genes of module 13 were mainly enriched in the regulation of protein polymerization and amyotrophic lateral sclerosis ( Figure 3A ).

By differential expression gene analysis, we compared the DEGs between FNB and sporadic NB, and we focused on SDHD, FGFR2, and NR4A1 ( Figure 3C ). Among them, SDHD may be associated with the development of family paraganglioma (42), which is also thought to originate from neural crest cells and has the same origin as NB. Therefore, the difference in SDHD protein expression between FNB and sporadic NB was further verified by immunohistochemistry ( Figure 3D ). In FNB, two cases (2/4, 50%) were positive for SDHD protein expression, whereas in sporadic neuroblastoma, all were negative for SDHD protein, and the immunohistochemical results further confirmed the difference in SDHD expression between FNB and sporadic neuroblastoma ( Supplementary Figure 2 ). In addition, as genetic determinants of familial disease, FGFR2 and NR4A1 were associated with familial breast cancer and familial Crohn’s disease, respectively, and we similarly identified significant differences between groups, but FGFR2 was expressed at lower levels in FNB and none in sporadic NB. Finally, we also identified the pseudogene RPL41P5 ( Figure 3B ), but the pseudogene is currently considered non-functional.

To understand the potential origin of NEs in FNB, we used pseudo-time trajectory analysis based on the Monocle 2 algorithm to distinguish whether there was a malignant transformation relationship between ECs, fibroblasts, Schwann cells, and NEs in FNB. The results showed that ECs have the potential to simultaneously evolve into NEs, fibroblasts, and Schwann cells ( Figures 4A, B , Supplementary Figure 3 ). However, the number of Schwann cells was too small for subsequent analysis. Heatmap hierarchical clustering analysis demonstrated the changes of DEGs with the progression of the pseudo-time, and with the increase of the pseudo-time, CHGB, CNTNAP2, ALCAM, BASP1, CALM2, CCND1, CD24, CADM1, CHGA, and other genes increased in expression with increasing pseudo-time ( Figures 4C–E ). We next applied partition-based graph abstraction (PAGA) trajectory analysis to further clarify the differentiation direction of NEs, fibroblasts, and ECs, and the results ( Figure 4F ) were consistent with the previous trajectory analysis, in addition, the differentiation potential of endothelial cells to NEs was stronger than that of endothelial cells to fibroblasts, and finally, we used the branched expression analysis modeling (Beam) algorithm to investigate the gene expression changes during differentiation( Figure 4G ). With the Beam algorithm, we clustered genes with similar expression patterns during differentiation into four clusters, from top to bottom, 4, 3, 1, and 2, respectively, from Pre-branch to cell fate1 representing the direction of differentiation of ECs to NEs, and from Pre-branch to cell fate2 representing the direction of differentiation of ECs to fibroblasts. We then performed GO and KEGG enrichment analyses on each cluster. We noticed that the DEGs expressed in the cluster 2 (ECs to NEs differentiation) were mainly enriched in collagen-containing extracellular matrix as well as PI3K-Akt signaling pathway ( Figures 4H, I ). The PI3K-AKT signaling pathway promotes the proliferation and growth of neuroblastoma (43).

To assess the metabolic heterogeneity of the various cell types in the FNB, we performed single-cell flux estimation analysis (scFEA) (44) on all cell subpopulations of the six samples, which showed that compared to all other cell subpopulations, NEs were highly correlated with the metabolism of spermine, oxaloacetate, and hypoxanthine. The results showed that NEs were highly correlated with the metabolism of spermine, oxaloacetate, and hypoxanthine, while the metabolism of IMP, tyrosine, choline, malate, and dUMP was at a lower level compared to all other cell subpopulations ( Figure 5A ). In addition, we did not observe high levels of lactate metabolism at the level of overall NEs, and the level of lactate metabolism in NEs was lower than that in T cells. In addition, it has been previously described that FNB has a higher proportion of immune cells, which also predicts a more active TME. Next, survival analysis ( Figure 5B ) was performed on selected top100 genes for each subpopulation of non-malignant cells, showing that a variety of immune effector cells, including natural killer (NK) cells, were not associated with prognosis, focusing only on Proliferating-MPs and Proliferating-T Cells, which were associated with poor prognosis, unlike previously observed (45, 46). While further subdivision of immune cells showed that the TME of FNB was mainly composed of M2-like macrophages, which is consistent with the results observed for sporadic neuroblastoma, CellChat analysis then was used to identify the interaction between ligand-receptor pairs and cells, and due to the low number of non-NE cells within the P2 group, P1 was only analyzed. The results showed that in FNB, NE cells communicate most closely with Monocytes via MIF-(CD74+CXCR4) and MIF-(CD74+CD44) ( Figure 5C ), and several ligand receptors including MDK-NCL have important roles in intercellular interactions in TME, which are potential therapeutic targets for NB ( Figures 5D–F ).

Cancer-associated fibroblasts (CAF) are abundant in the tumor microenvironment and are associated with a variety of tumor biological behaviors including tumor invasion, drug resistance, and immune regulation (47, 48). A previous study using markers for CAF in breast cancer (48)defined CAF-S1 as well as CAF-S4 in human NB (49), and to further investigate FNB-CAF, we mechanistically classified fibroblasts in F1 ( Figure 1D ) and obtained a total of four subpopulations, namely Fib-1, Fib-2, Fib-3 and Fib- 4 ( Figure 6A ). In FNB, we did not identify the above two subpopulations (in which αSMA, CD29, and PDGFRβ were not expressed in the matrix, while COL1A1 was expressed in all four subpopulations) ( Figures 6B-D ). The violin plot shows the highly expressed genes in each subpopulation ( Figure 6E ). Next, a survival analysis of the four subpopulations was performed, and the results showed that Fib-4 was associated with a good prognosis in NB ( Figures 6F, G ). Endothelial cells are also currently thought to be a source of CAF (50–52), so we further performed pseudo-time trajectory analysis of Fib-4 with endothelial cells in two FNB cases, and the Monocle 2 trajectory plot showed that in FNB, endothelial cells have the differentiation to Fib-4 potential ( Figures 6H-J ). In conclusion, a kind of prognostic related CAF subtype in FNB was identified by comprehensive bioinformatics, differentiated from endothelial cells, with good prognostic guidance in FNB.

Tumor tissues were obtained from children who were diagnosed with NB in the department of pediatric surgical oncology of children’s hospital of Chongqing medical university. Approval was obtained from the institutional ethics committee at our hospital. The fresh tumor tissues were rinsed with normal saline after surgical resection to remove blood cells. Then, the non-necrotic parts of tumor tissues (0.3 - 0.5 m³ per sample) were moved out, stored in 1 ml GEXSCOPE^® Tissue Preservation Solution (Singleron, China) and transported to the Singleron lab immediately.

The tumor tissues were washed with Hanks balanced salt solution (HBSS) for 3 times and cut into small pieces (1~2 mm), and put into 2 ml GEXSCOPE ^® In tissue dissociation solution (Singleron), stirred gently and continuously at 37 °C for 15 min. After digested, the samples were iltered with 40-micron sterile strainers and centrifuged at 800×g for 5 min. Then, resuspended the samples in 1 ml phosphate buffer (PBS) (HyClone) which added 2 ml GEXSCOPE^® red blood cell lysis buffer (Singleron) at 25 °C for 5-8 min to remove red blood cells. The above mixture was centrifuged at 500 × g for 5 minutes to precipitate cells. Then resuspended cells with PBS. Finally, stained the samples with Trypan Blue to evaluate the cell viability microscopically.

The concentration of single cell suspensions was adjusted to 1 x 10⁵ cells/mL. Then, the suspensions were loaded onto microfluidic chip, and scRNA-seq libraries were constructed according to the instructions of Singleron (GEXSCOPE^® Single-Cell RNA Library Kit, Singleron Biotechnologies). Individual libraries were diluted to 4 nM and pooled to sequence on an Illumina HiSeq X with 150-bp paired-end reads.

Raw reads from scRNA-seq were processed to generate gene expression matrixes using CeleScope (https://github.com/singleron-RD/CeleScope) v1.9.0 pipeline. Briefly, raw reads were first processed with CeleScope to remove low quality reads with Cutadapt v1.17 to trim poly-A tail and adapter sequences. Cell barcode and UMI were extracted. After that, we used STAR v2.6.1a (65) to map reads to the reference genome GRCh38 (ensembl version 92 annotation). UMI counts and gene counts of each cell were acquired with featureCounts v2.0.1 (66) software, and used to generate expression matrix files for subsequent analysis.

Scanpy v1.8.1 was used for quality control, dimensionality reduction and clustering under Python 3.7. For each sample dataset, we filtered expression matrix by the following criteria: 1) cells with gene count less than 200 or with top 2% gene count were excluded; 2) cells with top 2% UMI count were excluded; 3) cells with mitochondrial content 20% were excluded; 4) genes expressed in less than 5 cells were excluded. The raw count matrix was normalized by total counts per cell and logarithmically transformed into normalized data matrix. Top 2000 variable genes were selected by setting flavor = ‘seurat’. Principle Component Analysis (PCA) was performed on the scaled variable gene matrix, and top 20 principle components were used for clustering and dimensional reduction Batch effect between samples was removed by Harmony (67). Finally, UMAP algorithm was applied to visualize cells in a two-dimensional space.

Scanpy v1.8.1 was used for further clustering analysis of Fibroblast. Top 2000 variable genes were selected for PCA analysis. The first 6 principle components and resolution parameter 0.8 were used with louvain algorithm to generate 4 cell clusters.

Differentially expressed genes (DEGs) and pathway enrichment analysis, and cell type annotation

DEGs were determined as genes expressed in more than 10% of the cells in a cluster with an average log (Fold Change) of greater than 0.25, and were selected by scanpy rank_genes_groups function based on the Wilcox likelihood-ratio test with default parameters. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were used to investigate the potential functions of DEGs with the “clusterProfiler” R package version 3.16.1. P-value < 0.05 was considered statistically significant. The cell type identity of each cluster was determined with the expression of canonical markers found in the DEGs using SynEcoSys database. Violin plots displaying the expression of markers used to identify each cell type were generated by Seurat v3.1.2 Vlnplot.

Cell differentiation trajectory was reconstructed with Monocle2 (68). Next, highly-variable genes (HVGs) were used to sort cells in order of spatial‐temporal differentiation. DDRTree was used to perform dimension-reduction. Finally, the trajectory was visualized by plot_cell_trajectory function. To run PAGA (69), a symmetrized kNN-like graph based on PCA data was construct by scanpy.tl.paga, using the approximate nearest neighbor search within UMAP. For each partitioning, a PAGA graph was generated using the connectivity.

The InferCNV package (70) were used to evaluate the CNVs in NE cells Schwann cells, endothelial cells, and fibroblasts. T cells, B cells and myeloid cells were regarded as baselines to estimate the CNAs of malignant cells. Genes expressed in more than 20 cells were classified based on their loci on each chromosome. The relative expression values were centered to 1, using 1.5 standard deviations from the residual-normalized expression values as the ceiling. A slide window size of 101 genes was used to normalize the relative expression on each chromosome in order to remove the effect of gene-specific expression.

Each p- or q-arm level change can be simply converted to an equivalent CNV according to its location by considering genomic cytoband information. Each CNV was annotated as a gain or loss. Then, subclones containing the same arm-level CNVs were folded, and trees were reconstructed to represent the architecture of subclonal CNV.

Hotspot (37) was used to identify functional gene modules which illustrate heterogeneity within NEs subpopulations. Briefly, we used the ‘danb’ model and selected the top 500 genes with highest autocorrelation zscore for module identification. Modules were then identified using the create_modules function, with min_gene_threshold =15 and fdr_threshold = 0.05. Module scores were calculated by using calculate_module_scores function. The Jaccard similarity coefficient was calculated for comparing transcriptional similarity between cell types with hotspot modules genes (15).

CellChat (version 0.0.2) (71) was used to analyze the intercellular communication networks from scRNA-seq data. A CellChat object was created using the R package process. Cell information was added into the meta slot of the object. The ligand-receptor interaction database was set, and the matching receptor inference calculation was performed.

IHC was used to analyze the protien expression of SDHD in tumor tissue. Paraffin sections of 4 familial NB and 4 sporadic NB were collected. The paraffin sections were processed with baking, dewaxing, sodium citrate antigen repair, 3% H₂O₂ incubation and 0.5% BSA closure. Then sections were incubated with anti-SDHD (Abcam, ab189945), and corresponding secondary antibodies (ZSGB-Bio, PV-9001). DBA color development was carried out according to the instructions of the immunohistochemistry kit (ZSGB-Bio, ZLI 9019). Finally, the sections were stained with hematoxylin, sealed with neutral gum, and observed under a light microscope.

scFEA (v1.1.2) (45) is a computational method for inferring cellular metabolic fluxomes and metabolite abundances from scRNA-seq data using flux balance constraints based on novel probabilistic models. All cell types were calcuated the metabolic Fluxomes and Abundances by scFEA in python. After calculating metabolic flux and metabolite abundances, the 70 metabolic modules and all metabolites were selected for visualization by heatmap.