Single-cell landscape analysis reveals distinct regression trajectories and novel prognostic biomarkers in primary neuroblastoma.

Neuroblastoma (NB), which is the most common pediatric extracranial solid tumor, varies widely in its clinical presentation and outcome. NB has a unique ability to spontaneously differentiate and regress, suggesting a potential direction for therapeutic intervention. However, the underlying mechanisms of regression remain largely unknown, and more reliable prognostic biomarkers are needed for predicting trajectories for NB. We performed scRNA-seq analysis on 17 NB clinical samples and three peritumoral adrenal tissues. Primary NB displayed varied cell constitution, even among tumors of the same pathological subtype. Copy number variation patterns suggested that neuroendocrine cells represent the malignant cell type. Based on the differential expression of sets of related marker genes, a subgroup of neuroendocrine cells was identified and projected to differentiate into a subcluster of benign fibroblasts with highly expressed CCL2 and ZFP36, supporting a progressive pathway of spontaneous NB regression. We also identified prognostic markers (STMN2, TUBA1A, PAGE5, and ETV1) by evaluating intra-tumoral heterogeneity. Lastly, we determined that ITGB1 in M2-like macrophages was associated with favorable prognosis and may serve as a potential diagnostic marker and therapeutic target. In conclusion, our findings reveal novel mechanisms underlying regression and potential prognostic markers and therapeutic targets of NB.

Introduction

Neuroblastoma (NB) varies widely in clinical presentation and outcome and has been widely categorized according to histological features as high malignancy neuroblastoma (NBL), intermediate malignancy ganglioneuroblastoma (GNB), or benign ganglioneuroma (GN). NBL is the most common pediatric extracranial solid tumor, accounting for 8%–10% of pediatric malignancies. Despite multimodal therapy, NB may progress to malignancy; or may undergo spontaneous regression, a unique phenomenon that is not well understood. Indeed, tumors from patients with NB at the 4S stage have been observed to spontaneously differentiate into benign GN or normal cells after treatment with low doses of chemotherapy. Potential mechanisms for spontaneous regression in NB include alterations in epigenetic regulation, loss of telomerase activity, neurotrophin deprivation, and immunity regulation. These mechanisms provide a theoretical basis for NB differentiation therapy, whereby malignant tumor cells are induced to differentiate into normal or benign cells. However, additional mechanisms are likely to govern spontaneous regression, and exploring these mechanisms may provide new strategies for differentiation therapy.

Many recent studies have focused on characterizing the impact of intertumoral heterogeneity in NB, and molecular signatures for poor survival have been used to guide its diagnosis and treatment., Genetic molecular analysis of NB is now a relevant component of risk stratification, and biomarkers such as MYCN gene amplification and mutations in PHOX2B are used as prognostic indicators of NB. However, these biomarkers are not universally observed in NB: PHOX2B is predominantly mutated in familial NB and only accounts for 5% of hereditary cases, and the amplification of MYCN occurs in only 22% of sporadic NB. Therefore, there is an urgent and unmet clinical need to identify additional prognostic markers with high specificity and sensitivity to guide clinical stratification diagnosis and prognosis.

The tumor microenvironment for NB is a highly complex ecosystem composed of tumor cells, fibroblasts, endothelial cells, Schwann cells, pericytes, mesenchymal stromal cells, and immune infiltrating cells. The contribution of non-immune cells has also been recently recognized in NB, and cancer-associated fibroblasts have been demonstrated to share characteristics and pro-tumorigenic activity with mesenchymal stromal cells. Furthermore, although pediatric solid tumors are considered to possess low immunogenicity,, numerous epidemiological studies have also demonstrated the importance of tumor-infiltrating lymphocytes in NB,, including potential impacts on metastasis and prognosis. Thus, the potential prognostic value of tumor immune infiltrating cells, such as CD8+ T cells and dendritic cells, has attracted significant attention. Nonetheless, our understanding about the specific molecular mechanism of these immune infiltrating cells in NB is limited.

In recent years, single-cell RNA sequencing (scRNA-seq) analysis has shown great promise in dissecting the transcriptional features and functional roles of different cell populations to elucidate pathways of cancer development and progression., Here, we performed scRNA-seq on different stages of NB and peri-tumoral adrenal gland (AG) to investigate the developmental trajectory for spontaneous regression of NB and explore the potential value of the immune infiltrating cells. Our findings reveal novel mechanisms underlying spontaneous regression of NB, as well as potential prognostic markers and novel therapeutic targets for NB. Further understanding of the transcriptional characteristics associated with clinical heterogeneity in NB should assist in the identification of meaningful biomarkers for the risk stratification and treatment of NB.

Materials and methods

Sample collection and processing

Tissues were collected from the Children's Hospital of Chongqing Medical University. The use of clinical samples was reviewed and approved by the Institutional Review Board of Chongqing Medical University, Chongqing, China. Written informed consent forms were signed by parents or legal guardians, and assent, when appropriate, was obtained from the patients. The study population consisted of 20 samples. Among them, 11 samples were NBL (high malignancy; groups 1–4), including three stage IV without MYCN amplification (group 1), two stage IV with MYCN-amplification (group 2), three stage IVS (group 3), and three stage I/II (group 4); three samples were GN (benign; group 5); and three samples were GNB (intermediate malignancy; group 6). Three AG samples comprising group 7 served as the control group.

All patients were treatment-naïve and diagnosed through biopsy. The clinical characteristics of the selected patients are shown in Table S1. All samples were collected by surgical resection followed by a brief rinse with saline solution to remove contaminating blood cells. The non-necrotic parts of tumor samples (0.3–0.5 m3 per sample) were dissected out and stored in 1 ml GEXSCOPE® Tissue Preservation Solution (Singleron, China).

Tissue dissociation and preparation

Fresh tumor samples were transported expeditiously to the Singleron lab. The tumor tissues were washed three times with Hanks Balanced Salt Solution (HBSS), cut into 1–2 mm pieces, and digested in 2 ml GEXSCOPE® Tissue Dissociation Solution (Singleron) at 37 °C for 15 min with constant but gentle agitation. The digested samples were filtered using 40-micron sterile strainers and centrifuged at 800×g for 5 min. The pellets were then resuspended in 1 ml phosphate-buffered saline (PBS; HyClone, United States), and 2 ml GEXSCOPE® red blood cell lysis buffer (Singleron) was added at 25 °C for 10 min to remove the red blood cells. The tumor cell-containing solution was centrifuged at 500×g for 5 min, and the cells were resuspended in PBS. The samples were stained with trypan blue (Sigma, United States) to verify cell viability.

Single-cell RNA sequencing

Single-cell suspensions with 1 × 105 cells/mL in PBS were prepared. Then, the suspensions were loaded onto microfluidic devices, and scRNA-seq libraries were constructed according to the Singleron GEXSCOPE™ protocol in the GEXSCOPE™ Single-Cell RNA Library Kit (Singleron Biotechnologies). Individual libraries were diluted to 4 nM and pooled for sequencing. Pools were sequenced on an Illumina HiSeq X with 150 bp paired-end reads.

Primary analysis of raw read data

Raw reads were processed with fastQC and fastp to remove low-quality reads. Unique molecular identifiers (UMI) and cell barcodes were extracted after filtering out reads without poly-A tails. Poly-A tails and adapter sequences were then trimmed by cutadapt, and the reads that met quality control standards were matched to the reference genome GRCh38 (ensembl version 92 annotation) using STAR. Gene counts and UMI counts were grouped by featureCounts software to generate expression matrix files for downstream analyses.

Quality control, dimension-reduction, and clustering

Cells with gene counts between 200 and 5000 and UMI counts below 30,000 were filtered, and cells with more than 30% mitochondrial content were removed. After filtering, 89,882 cells were reserved for subsequent analyses, with on average 1097.5 genes and 2480.9 UMIs per cell. Then Seurat v2.3 was used for dimension-reduction and clustering. NormalizeData and ScaleData were used to normalize and scale all gene expression values. The top 2000 variable genes were selected for principal component analysis (PCA) by FindVariableFeatures. FindClusters was used to separate the cells into nine clusters using the top 20 principal components and a resolution parameter at 1.0. The resolution was then set to 1.2 to isolate subcluster cell types. A uniform manifold approximation and projection (UMAP) algorithm was performed to visualize the cells in two-dimensional space.

Differentially expressed genes (DEGs) analysis

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.2, and were selected by Seurat FindMarkers based on the Wilcox likelihood-ratio test with default parameters.

Cell type annotation

The cell type identification of each cluster was manually annotated according to the expression of canonical markers found among the DEGs combined with knowledge from literature. Violin plots that exhibit the expression of cell-type markers were generated by Seurat Vlnplot.

Pathway enrichment analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed using the “clusterProfiler” R package version to explore the potential functions of DEGs between two clusters. Pathways with a P value of less than 0.05 were considered significantly enriched. GO gene sets included molecular function (MF), biological process (BP), and cellular component (CC) categories.

Trajectory analysis

Pseudotime trajectory analysis was performed by using Monocle2 to determine map conversion trajectories of cell subtypes among endothelial cells, fibroblasts, NEs, and steroidogenic cells. To construct the conversion trajectories, highly variable genes were selected from clusters of endothelial cells, fibroblasts, NEs, and steroidogenic cells using Seurat version 3.1.2 FindVairableFeatures, and dimension-reduction was employed using DDRTree. plot_cell to visualize the conversion trajectory.

RNA velocity

For RNA velocity assessment, BAM files containing gene sets of NE cells, ECs, fibroblasts, and the reference genome were analyzed with velocyto and scVelo in python using default parameters. The results were visualized on UMAP plots from Seurat clustering analysis for consistency.

scRNA-seq based CNA detection

We used the InferCNV package to evaluate the CNAs in NE cells, Schwann cells, steroidogenic cells, endothelial cells, and fibroblasts. T cells, B cells and myeloid cells were considered as baselines to estimate the CNAs of malignant cells. Genes expressed in over 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. To remove the effect of gene-specific expression, a slide window size of 101 genes was used to normalize the relative expression on each chromosome.

CNV trees

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. After data conversion, subclones containing the same arm-level CNVs were folded, and trees were reconstructed to represent the architecture of subclonal CNV. For data visualization, we employed the uphyloplot2 algorithm.

Cell–cell interaction analysis (CellPhoneDB)

Cell–cell interaction between NEs and immune cell types was predicted according to known ligand–receptor pairs by Cellphone DB version. The permutation number was defined by calculating the null distribution of average ligand–receptor pair expression in randomized cell identities and was set to 1000. The expression threshold for individual ligands or receptors was a cutoff value based on the average log gene expression distribution of all genes in each cell type. Interaction pairs were predicted to be significant if they had P value <0.05 and average log expression >0.1, and were visualized using the circlize (0.4.10) R package.

Single-cell regulatory network inference and clustering (SCENIC)

A transcription factor network was constructed using the SCENIC R toolkit with the scRNA expression matrix and transcription factors in AnimalTFDB. To predict the co-expression network of regulators and targets, we used the GENIE3 package. RcisTarget was used to search transcription factor binding motifs. Genes involved in the predicted regulatory network were defined as a gene set. We then used the AUCell package to calculate the AUC value of the gene set and evaluate the activity of the regulatory network in cells.

Expression programs analysis

We used the cNMF algorithm to extract transcriptional programs, taking the top50 genes as the meta-signature and calculating the score of each program for each cell. Six meta-programs were calculated, and hierarchical clustering was performed based on the personal correlation between each program.

Jaccard similarity analysis

Signature genes of two cell types were used to compare the transcriptional similarity and calculate the Jaccard similarity coefficient. We assessed transcriptional similarities between six meta-programs and seven markers of cell status by calculating the Jaccard similarity coefficient using the first 50 marker genes.

Survival analysis

The TARGET database neuroblastoma project (n = 154) was utilized to evaluate the prognostic value of each signature. Bulk expression data were downloaded by the R package “TCGAbiolinks”. Kaplan–Meier survival curves, Hazard ratios (HRs) and P values were created with the R survival package to present differences in survival time.

qRT-PCR

Total RNA extraction was performed using TRIzol reagent (Invitrogen). Then, 200 ng of RNA, quantified by Bio-Rad CFX Connect (Bio-Rad, USA), was reverse transcribed using an Evo M-MLV RT Mix Kit with gDNA Clean for qPCR (ACCURATE BIOLOGY, China). The cDNA was amplified with target gene-specific primers using a SYBR Green Premix Pro Taq HS qPCR Kit II (ACCURATE BIOLOGY). The sequences of the primers (designed and synthesized by Shanghai Sangon Biotech company) are listed in Table S5. The expression levels of TUBA1A and STMN2 relative to ACTIN levels were calculated using the 2−ΔΔCT method.

Immunohistochemistry

Immunohistochemistry was used to analyze the expression of CCl2, ZFP36, TUBA1A, STMN2, and ITGB1 (Table S6) in tumor biopsy samples. Paraffin-embedded, serial tissue sections with a thickness of 4.5 μm were deparaffinized and subjected to dewaxing, sodium citrate antigen repair, 3% H2O2 incubation at room temperature to block enzymes, 0.5% BSA closure, and dropwise addition of primary antibody overnight at 4 °C. Then, secondary antibody (Table S6) and immunohistochemistry kit reaction enhancement solution were sequentially added for DBA color development, followed by hematoxylin re-staining, neutral gum sealing and observation under a light microscope.

Statistical analysis

SPSS 25.0 was used for statistical analysis of qRT-PCR results, and statistical significance was determined with a Mann–Whitney test. The P values were calculated, and P < 0.05 was considered statistically significant.

Results

Primary NB displays varied cellular constitution with a potential role for differentiation of cell types

To investigate the developmental trajectory for spontaneous regression of NB, we performed scRNA-seq analysis on primary tumors from 17 NB patients and three peritumoral adrenal gland tissue (AG), which were divided into seven groups (Fig. 1A, Table S1). Among them, 11 samples were from high malignancy NBL tumors, including three stage IV without MYCN-amplification (group 1); two stage IV with MYCN-amplification (group 2); three stage IVS (group 3); and three stage I/II (group 4). Three samples were from benign GN tumors (group 5); and three samples were from moderate malignancy GNB tumors (group 6). The three AG samples comprised group 7.

Figure 1

Cellular constitution of primary NB and AG tissues. (A) Schematic diagram of the study roadmap for 17 neuroblastoma samples and three peri-tumoral adrenal gland tissues. (B) The UMAP plot of the nine main cell clusters in primary NB and AG tissues. (C) The UMAP plot of each cell type in 20 samples. (D) Violin plots of the expression of marker genes in nine cell types. (E) Relative proportion of each cell type in four groups. (F) Relative proportion of each cell type in 20 samples. NB, neuroblastoma; AG, adrenal gland tissue; GN, ganglioneuroma; GNB, ganglioneuroblastoma; NE cells, neuroendocrine cells; ECs, endothelial cells.

After quality control, including doublet removal, we obtained 89,882 cells from primary NB and peritumoral AG, including 10,683 cells from GN samples, 16,982 cells from GNB samples, 47,814 cells from NBL samples, and 14,403 cells from peritumoral AG. The gene numbers (Fig. S1a) and total UMIs (Fig. S1b) of each sample are shown in box plots. Next, we used uniform manifold approximation and projection (UMAP) to visualize cell clusters, which revealed nine main cell types, including neuroendocrine cells (NEs), steroidogenic cells, Schwann cells, T cells, B cells, myeloid cells, fibroblasts, ECs, and pericytes (Fig. 1B,C). A violin plot shows the expression of canonical markers in the 9 clusters (Fig. 1D). Consistent with the literature, the NEs were distinguished by high expression of SNAP25, SYT1, DCX, CHGB, CHGA, NPY, and SCG2,30, 31, 32 while the steroidogenic cells were characterized by high expression of FDX1, CYP11B1, CYP11A1, and MGARP., Moreover, the T cells specifically expressed CD2, TRBC2, CD3D, and CD3E,; the B cells had high expression of MS4A1, CD79B, and IGHM,; and the myeloid cells had high expression of LYZ, C1QC, CD1C, CD68, and MRC1., We also identified fibroblasts by high expression of DCN, LUM, COL1A1, and COL1A2,; ECs by high expression of PECAM1, VMF, and CLDN5, and pericytes by expression of RGS5, ACTA2, and TAGLN., The Schwann cells were characterized by the expression of PLP1, CNP, S100B, and MPZ.

We further compared the proportion of different cell types in each patient and found significant differences in the cellular distribution between groups, especially for NEs, ECs, steroidogenic cells and fibroblasts (Fig. 1E, F and Table S2). The data demonstrate that the proportion of steroidogenic cells in stage IVs was significantly higher than that in other stages, suggesting that steroidogenic cells may play a role in the development of NB. Moreover, in different pathological types, with the conversion of NBL to GNB and to GN, the proportion of NEs was decreased, but the proportions of ECs and fibroblasts were increased. Therefore, compared with malignant NB, the intermediate malignancy and benign NB had more fibroblasts and ECs, but fewer NEs. Because fibroblasts have the potential to differentiate into NEs, we speculated that the differentiation of NEs into ECs and fibroblasts could potentially provide a basis for understanding regression trajectories that occur upon spontaneous regression.

NEs represent malignant tumor cells of NB tumors

To further classify malignant and non-malignant cells in primary NB and peritumoral AG, we used the inferred CNV algorithm to determine the clonal structure of the cells. Our results reveal that there were more CNVs in NEs than in other cell types. As shown in Figure 2A, the NEs of all NBL tumors had chromosome 17q gain, while the MYCN-amplified group had chromosome 1p deletion. These chromosomal aberrations have previously been reported in NB tumors., Thus, these results suggest that NEs may represent malignant cells in primary NB tumors, which is supported by the results of other studies (Fig. S2).

Figure 2

Trajectory of malignant tumor cells. (A) Copy-number variation (CNV) analysis was performed for NEs. (B–D) The Monocle 2 trajectory plot shows the dynamics of ECs, fibroblasts, NEs and steroidogenic cells. (E) Heatmap hierarchical clustering shows genes and pathways regulated during the NE pseudotime trajectory. (F) DEGs were tracked along the pseudotime curve. (G) The CNV trees for fibroblasts, ECs, NEs. (H) UMAP plots are shown for fibroblasts in all samples, and for cells that were classified into four subclusters.

Further examination of the results for steroidogenic cells revealed chromosome 2 aberrations (Fig. S1c), while Schwann cells had chromosome 19 deletions (Fig. S1d). As the chromosomal aberrations of steroidogenic cells were similar in different NB tumor types, these results could not distinguish whether chromosome 2 aberration is associated with malignancy in steroidogenic cells. Additionally, Schwann cells did not show tumor-specific chromosome aberrations. Thus, we believe chromosome 19 variations might represent a characteristic of Schwann cells.

Because ECs (Fig. S1e) and fibroblasts (Fig. S1f) exhibited few CNVs, we considered these two cell types as non-malignant cells. Notably, all results show deletions in chromosome 6 because CNV analysis was performed using myeloid cells as a control, and myeloid cells display obvious gain in chromosome 6 without CNVs in other chromosomes. Thus, taken together, our results indicated that NEs represent malignant tumor cells of NB, and other types of cells represent non-malignant tumor cells.

Malignant tumor cells can differentiate into fibroblasts

To determine whether ECs, fibroblasts, steroidogenic cells, and NEs contribute to malignant transformation, we sought to establish the relationship among these four cell types. First, we utilized pseudotime trajectory analysis based on the Monocle 2 algorithm to estimate the NE maturation course. The results suggest that the ECs displayed the potential to differentiate into fibroblasts and NEs, and fibroblasts displayed the potential to differentiate into NEs (Fig. 2B–D). By contrast, the trajectory analysis showed a low tendency of steroidogenic cells to differentiate to fibroblasts and NEs, suggesting that steroidogenic cells may develop independently and have little correlation with differentiation and development of NEs. Further analysis by heatmap hierarchical clustering identified DEGs that were increased or decreased along the pseudotime curve (Fig. 2E, F). Additionally, a clonality tree for fibroblasts, ECs, and NEs based on inferred CNV results suggests that multiple canonical and non-canonical CNVs in fibroblasts were more similar to those of NEs, rather than those of ECs (Fig. 2G). These results support the idea that the transformation between fibroblasts and NEs may be more likely than the transformation between NEs and ECs.

To further study the trajectory of NEs, we segregated fibroblasts into four cellular subclusters, Fib-C1, Fib-C2, Fib-C3, and Fib-C4, which we visualized by UMAP (Fig. 2H). We analyzed the developmental trajectories of ECs and NEs in the four subclusters by pseudotime trajectory analysis. Our results reveal that Fib-C1, Fib-C2, and Fib-C3 subclusters of fibroblasts and ECs showed a sequential commitment to NEs (Fig. 3A). We next applied RNA velocity with velocyto to further predict the trajectory of NEs, ECs, and the four subclusters of fibroblasts. The Fib-C2 subtype was identified as the terminal state of these cells (Fig. 3B,C). To further understand the potential functions of the Fib-C2 subcluster, we performed GO and KEGG enrichment analyses. The results show that DEGs in the Fib-C2 subcluster were mainly enriched in the IL-17 and TNF signaling pathways (Fig. S3a, b).

Figure 3

Malignant tumor cells can differentiate into fibroblasts. (A) The dynamics of ECs, NEs, and the four subclusters of fibroblasts according to pseudotime trajectory analysis. (B, C) RNA velocity analysis of the trajectory of NEs, ECs, and fibroblasts. (D) UMAP plots are shown for NEs in all samples, and for cells that were classified into six subclusters. (E) Cell network for NEs and fibroblasts form CellPhoneDB analysis. (F) The expression of C2-related molecules at the protein level in tumor cells of GNB as evaluated by immunohistochemistry. The arrow points the positive cells and 40× pictures were enlarged from the 20× picture.

Next, we divided NEs into six cellular subclusters (NE-C1, NE-C2, NE-C3, NE-C4, NE-C5, and NE-C6) that we evaluated by CellPhoneDB analysis (Fig. 3D). The results suggest that NE-C1 has close communication with fibroblasts, suggesting that NE-C1 subtypes may differentiate into the Fib-C2 subcluster (Fig. 3E). Functional enrichment for marker genes suggested that the NE-C1 subgroup is associated with ribosome-related pathways, such as “cytosolic ribosomes” and “structural constituent of ribosomes” (Fig. S3c, d). Finally, we conducted immunohistochemical analysis to verify the expression of Fib-C2-related markers in tumors and found that CCL2 and ZFP36 were highly expressed in malignant tumor cells (Fig. 3F), thus supporting the existence of a regression relationship in which NEs can differentiate toward the Fib-C2 subtype. Collectively, the above results strongly suggest that the NE-C1 subtype of malignant tumor cells has the capacity to differentiate into the Fib-C2 subcluster of fibroblasts, providing a pathway for the spontaneous regression of NB.

Intra-tumoral heterogeneity of malignant NB tumors

Based on the aggregated CNV results, we generated a clonality tree for the NBL groups using the UPhyloplot2 plotting algorithm (Fig. 4A). The results indicate multiple canonical and non-canonical subclones in each NBL group. The canonical CNV subclones mainly contained chromosome gains on 17q, 2p, and 1q, while some non-canonical subclones predominantly had chromosome gains at 7q, 20p, and 4p. Our results also demonstrate that the CNV score of chromosome 1 in the MYCN-amplified group was significantly different compared with those of other groups (Fig. S3e). Furthermore, the CNV score of chromosome 17 was significantly higher in the stage IV group than in other groups (Fig. S3f). Thus, these canonical and non-canonical CNV subclones could contribute to tumor development and progression.

Figure 4

Intra-tumoral heterogeneity of primary NB tumors. (A) CNV trees of single cells from NBL groups. (B) Violin plots and Kaplan–Meier survival curves of PAGE5 and ETV1. (C) Heatmap showing the correlation of all 300 signatures derived from NMF analysis of primary NB tumors; six highly correlated meta-programs are highlighted. (D) Jaccard similarities of six meta-programs with the signature expression of seven cell types. (E) Violin plots present the scores of six meta-program signatures from CNV analysis in malignant cells of 11 NBL samples.

Because NEs are mainly composed of the NE-C1 and NE-C2 subclusters, we further compared the DEGs of NE-C1 and NE-C2 subtypes to screen for biomarkers of high-risk NBL. PAGE5 was significantly overexpressed in the non-amplification MYCN group at stage IV but was extremely under-expressed in the other groups (Fig. 4B). Additionally, survival analysis confirmed that PAGE5 is a potential prognostic biomarker (Fig. 4B). We also determined that ETV1 was highly expressed in the MYCN-amplification group at stage IV and was lowly expressed in the other groups. Unfortunately, there are currently no data suggesting that ETV1 serves as a potential prognostic biomarker. Therefore, our results suggest that PAGE5 is both a potential biomarker for the non-amplification MYCN group at stage IV and a potential prognostic biomarker for NBL, while ETV1 is a potential biomarker for the MYCN-amplification group at stage IV.

Malignant NBL tumors share common prognostic biomarkers: TUBA1A and STMN2 expression

We next conducted non-negative matrix factorization (NMF) analysis to explore the transcriptional spectrum of malignant NE cells. Six meta-programs identified 300 meta-genes that were preferentially expressed in subpopulations of NEs. Hierarchical clustering analysis verified that the six meta-programs had high concordance (Fig. 4C). Furthermore, several functions were represented among the six meta-programs, including “regulation of neuron development” (meta-program1; e.g., STMN2 and TUBA1A), “pluripotency of stem cells” (meta-program2; e.g., ZFHX3 and RIF1), “axonogenesis” (meta-program3; e.g., NEFL and MAP1B), “cell cycle” (meta-program4; e.g., CCNA2 and CCNB1), “signal transduction” (meta-program5; e.g., JUN), and “ribosome composition” (meta-program6; e.g., RPS27 and RPL23). Signature gene scoring analysis of NEs for the six meta-programs demonstrated that meta-program 6 received the highest score (Fig. 4D), suggesting that signatures in meta-program 6 may be preferentially expressed in tumor cells. We further evaluated the transcriptional similarity between the six meta-programs and the signatures of the seven cell types and confirmed that the malignant meta-programs had the most similarity with fibroblasts (Fig. 4E), supporting the possibility that malignant tumor cells may have a transformational relationship with fibroblasts.

To further analyze the relationship between the malignant signatures of the six meta-programs and clinical prognosis, we employed the TARGET database neuroblastoma project (n = 154). Survival analyses revealed that a high score in meta-program 4 predicted poor survival in the TARGET neuroblastoma database (Fig. 5 panels a & b, S4), indicating that malignant signatures may predict the aggressiveness of malignant cells in NB. We further identified prognostic signatures in each meta-program by screening 11 biomarkers associated with tumor prognosis (Fig. 5C, S5), including STMN2, TUBA1A, NTRK3, RTN4, CCNB1, CCNA2, CCNB2, UBE2C, NUF2, CKAP2, and JUN. We then used violin plots to analyze the distribution of the 11 prognostic signatures in each group (Fig. S6) and identified STMN2 and TUBA1A as potential diagnostic markers (Fig. 5D). Expression of TUBA1A and STMN2 proteins was verified at the protein level by immunohistochemistry of 56 tumor tissues (Fig. 5E). Interestingly, our results suggest that these genes were highly expressed in the NBL group with poor prognosis, which is opposite of the results from the TARGET database. Analysis of RNA expression by qRT-PCR results confirmed the immunohistology results (Fig. 5F), demonstrating that RNA expression of TUBA1A and STMN2 in high malignancy NBL was higher than that in peri-tumoral AG.

Figure 5

Characteristics of survival prognosis for malignant tumors. (A) Hazard Ratio Plot for six meta-programs. (B) Kaplan–Meier survival curves of the meta-program 4. (C) Hazard Ratio Plot for the eleven signatures. (D) Violin plots and Kaplan–Meier survival curves of TUBA1A and STMN2. (E) The expression of TUBA1A and STMN2 at the protein level in stage II and MYCN-amplification samples as evaluated by immunohistochemistry. All pictures are at 20× magnification. (F) The relative expression level of TUBA1A and STMN2 by qRT-CR. ∗∗ represent as P < 0.01. ∗∗∗ represent as P < 0.001. (G) Heatmap showing transcription factors in NEs from non-negative matrix factorization (NMF) analysis.

For a more detailed understanding of the transcriptional pathways that may mediate prognosis in NB, we next carried out SCENIC analysis of potential transcription factors in NEs. The most highly represented transcription factors included SOX4, ETV1, and UQCRB (Fig. 5G). Furthermore, SCENIC analysis predicted that SOX4 regulates the expression of TUBA1A and STMN2. These results indicate that the malignant signatures established by this analysis show common patterns of intra-tumoral transcriptional heterogeneity in NB tumors.

Complex intercellular networks mediate the communication between NEs and immune cells

By comparing the proportion of different immune cells in each group, we found that there were significant differences in cell types among tumor groups, especially for myeloid cells and T cells, and that the myeloid cells mainly constituted the immune microenvironment of NB tumors in terms of number. In the NBL group, there were significant differences in the proportion of T cells and myeloid cells in the stage I/II, stage IV without MYCN amplification and stage IV with MYCN amplification groups. The ratio of T cells and myeloid cells decreased gradually from the stage I/II to stage IV without MYCN amplification and then to the stage IV with MYCN amplification groups, suggesting that T cells and myeloid cells may play roles in inhibiting tumor growth.

To evaluate the contribution of immune cells in the NB tumor prognosis, we subdivided T cells, B cells, and myeloid cells (Fig. 6A). T cells were classified into five subclusters: naïve T cells, helper T (Th) cells, CD8+ effector T cells, proliferating CD8+ effector T cells, and Tregs. The B cells were categorized into four subclusters: naïve B cells, plasma cells, proliferative B cells, and proliferative plasma cells; and the myeloid cells were classified into seven subclusters: proliferative macrophages, M1-like macrophages, M2-like macrophages, non-classical monocytes, classical monocytes, cDCs, and pDCs. Further subdivision of immune cells revealed that myeloid cells were mainly composed of M2-like macrophages (Table S3), indicating that M2-like macrophages mainly constituted the tumor immune microenvironment in terms of number, suggesting that M2-like macrophages may play an important role in inhibiting tumor growth.

Figure 6

Complex intercellular communication networks in NEs and immune cells. (A) UMAP plot for T cells, B cells, and myeloid cells in all samples. (B) Cell network for major immune cells from CellPhoneDB analysis. (C) Heat map showing correlations between major immune cells. (D) Dot plot showing ligand–receptor pair interactions between NEs and major immune cells; ordinate is ligand–receptor pairs, abscissa is NEs and major immune cells. (E) Hazard Ratio Plot for nine signatures. (F) The expression of ITGB1 at the protein level in stage II and MYCN-amplification samples as evaluated by immunohistochemistry. The arrow indicates positive cells; 40X pictures were enlarged from 20X pictures. (G) SCENIC analysis was used to construct transcriptional factor networks for M2-like macrophages, DCs, and Th cells.

Next, we used CellPhoneDB analysis to identify ligand–receptor pairs and molecular interactions with NEs (Fig. 6B, C). The results suggest that compared with other immune cells, NEs had the closest communications with M2-like macrophages, which was through a large number of receptor–ligand pairs, including CD22–PTPRC and NGF–NGFR (Fig. 6D). In addition, NEs communicated with T cells, B cells, and DCs. NEs and Th cells communicated with each other mainly through CD40LG–CD40, while NEs and effector T cells communicated via GZMA–PARD3, and NEs and naïve B cells communicated via SELL-PODXL. NEs also communicated with DCs via CD48–CD244A. Thus, NEs engage in a repertoire of interactions with immune cells, and key signatures for the interaction mechanism between NEs and immune cells represent potential therapeutic targets.

ITGB1 is a potential therapeutic target for NB

Lastly, to verify the relevance of the immune subtypes in NB outcome, we analyzed possible correlations of the distribution and signaling characteristics of immune cell subtypes with NBL prognosis. Our results showed that the distribution of CD8+ T cells and Th cells was correlated with the malignancy of NB tumors, as more of these cell types were found in the early stages (Table S4), which is consistent with the results of other studies. The pattern of distribution of DCs also was consistent with those observed in previous studies, including a positive correlation with T cell infiltration in NBL at both the transcriptional and protein levels and association with favorable prognosis (Table S3). However, the distribution of B cells and monocytes did not correlate with the prognosis of NB.

We further analyzed the signatures associated with survival prognosis in subtypes of DCs and T cells and found that DCs could predict the survival of NBL patients with the expression of ABHD6 and BOD1L1, and that in Tregs, the expression of UTP23 provided favorable prognostic value in NBL (Fig. 6E, S7). RPA3 and CYC1 in Th cells also were able to predict the clinical outcome of NBL, while the expression of IGHG4 in proliferating CD8+ effector T cells correlated with poor prognosis in NBL. We further analyzed the distribution of these prognostic signatures in each group and found that only BOD1L1 and CYC1 were distributed among all NBL groups, especially in the high-risk group of NBL, suggesting that these two prognostic signatures may serve as potential diagnostic biomarkers for high-risk NB (Fig. S8).

Next, we evaluated macrophage markers that are associated with NEs. Notably, ITGB1, NCAM1 and ERBB3 were associated with NEs and correlated with favorable prognosis (Fig. 6E, S7). Moreover, ITGB1 were distributed in most groups (Fig. S8), suggesting that the signature may serve as potential therapeutic targets. Through immunohistochemical analysis, we found that ITGB1 was highly expressed in patients with II stage, compared with patients with MYCN-amplification (Fig. 6F). As further verification, we constructed a transcription factor network to explore the co-expression of transcription factors and putative target signatures in NEs and immune cells using SCENIC analysis and found that the most represented transcription factors were SPI1, JUNB, FOSB, ATF3, JUN, MEF2C, and CEBPD, among which RAD21 can regulate ITGB1 expression (Fig. 6G). Collectively, our results provide potential diagnostic and/or prognostic signatures, as well as novel therapeutic targets of NB.

Discussion

NB is a life-threatening pediatric disease with heterogeneous clinical presentation. In a subset of patients who do not receive therapy, the primary retroperitoneal tumor is identified by surgical histopathology as intermediate malignant GNB while the mediastinal metastases and cervical lymph node metastases are identified as benign GN. This clinical phenomenon indicates that the NB tumor can undergo spontaneous transformation or regression. Thus, we speculated that there is a trajectory of spontaneous regression from NB to GNB or to GN that could inform the development of novel therapeutic approaches. To elucidate the potential mechanisms that underlie spontaneous regression, we performed scRNA-Seq on 17 NB samples with distinct histological classification and three AG samples. We hypothesized that the distribution of steroidogenic cells, ECs, fibroblasts, and NEs may differ among NB groups. Based on pseudotime trajectory and clonality tree analyses, we revealed a pathway by which malignant tumor cells may have the potential to differentiate into fibroblasts. The similarity in specific transcriptional features supports the idea that malignant tumor cells may have a transformational relationship with fibroblasts. Using RNA velocity and CellPhone DB analysis, we demonstrated that the C1 subgroup of malignant tumor cells may differentiate toward the C2 subtype of fibroblasts, with highly expressed CCL2 and ZFP36. Immunohistochemical results confirmed that these Fib-C2-related markers are also highly expressed in malignant tumor cells. Therefore, our results strongly indicate that the spontaneous regression of NB may be associated with the differentiation of the C1 subgroup of malignant tumor cells to the C2 subtype of fibroblasts. To the best of our knowledge, our findings are the first to provide a trajectory for NB regression, though further investigation is warranted.

By uncovering the transcriptional heterogeneity of malignant tumor cells, we also identified six meta-programs with 300 meta-genes that were preferentially expressed in a subpopulation of NEs, with two independent prognostic factors (STMN2 and TUBA1A) indicated by survival analyses using the TARGET database. Further immunohistochemical results showed that TUBA1A and STMN2 were highly expressed in the NBL group with poor prognosis. Moreover, results of qRT-PCR supported the immunohistochemical results. The expression of these two signatures has been associated with poor prognosis in other tumors., Furthermore, SCENIC analysis predicted that SOX4 regulates the expression of TUBA1A and STMN2, suggesting that these two signatures may provide promising diagnostic and/or prognostic biomarkers, as well as novel therapeutic targets of NB patients. Though these results provide new insight into the molecular mechanisms that underlie disease progression, an expanded sample size is needed to confirm the prognostic value of TUBA1A and STMN2 in NB.

Lastly, we established a complex intercellular communication network to elucidate the cellular tumor microenvironment in NB tumors. Notably, among the immune cells, M2-like macrophages had the closest communication with NEs through a large number of receptor–ligand pairs. Interestingly, we found more M2-like macrophages in stage I/II than in stage IVs, which is in contradiction with other studies suggesting that M2-like macrophages are associated with poor prognosis in NB. We speculate that different subtypes of M2-like macrophages are represented in NB, though further studies are needed to clarify specific phenotypes of M2-like macrophages. We further investigated the prognostic value of signatures in M2-like macrophages. Our results demonstrate that ITGB1 in M2-like macrophages were associated with NEs and correlated with favorable prognosis. This is in contrast with other studies in which the expression of ITGB1 has been associated with poor prognosis in other tumors, suggesting that the role of this gene may vary for different tumor types. Our immunohistochemical analysis confirmed that high expression of ITGB1 protein was associated with favorable prognosis in the NBL groups. Furthermore, SCENIC analysis predicted that RAD21 regulates the expression of ITGB1, suggesting that ITGB1 might serve as a novel diagnostic and therapeutic target for NB. Further studies with extended sample sizes are needed to elucidate the roles of these markers in NB.

In summary, we conducted an extensive analysis of the transcriptomic landscape of NB through scRNA-seq of a panel of different pathological subtypes. Our findings suggest that the C1 subgroup of malignant tumor cells may differentiate into a C2 subcluster of fibroblasts, which reveals a novel mechanism of spontaneous NB regression. Furthermore, we demonstrated that the molecular characterization of malignant and immune cells associated with NB prognosis may provide promising biomarkers for diagnosis and an alternative mechanism for therapy of NB.

Author contributions

S.W. was responsible for the overall conception and design of the study and for revising the manuscript. Q.Q.L., Z.N.W., and Y.J. performed most of the experiments and the data analysis, collected the data, and wrote the manuscript. F.L.S., L.J.C., and X.B.D. collected the samples and the data. Q.L. and Y.B. coordinated the experiments. J.W.Z., L.P., and Z.Z.Z. contributed to the sample collection. M.Z.Z. was responsible for design and revision of immune analyses. T.C.H. revised the manuscript. All authors reviewed, edited and approved the manuscript.

Conflict of interests

We declare no competing interests.

Funding

This work was supported in part by research grants from the Key Project of “Research on Prevention and Control of Major Chronic Non-Communicable Diseases”, the , (No. 2018YFC1313000, 2018YFC1313004), the General Clinical Medical Research Program of Children’s Hospital of Chongqing Medical University (No. YBXM-2019-003), and the (No. cstc2019jscx-msxmX0220).

Ethics declaration

The clinical sample study protocol was reviewed and approved by the Institutional Review Board of the Children's Hospital of Chongqing Medical University. Written informed consent was obtained from the parents or legal guardians, and assent, when appropriate, was obtained from the patients.

Availability of data and material

The data that support the findings of this study are managed by Shan Wang; restrictions apply to the availability of these data, which are not publicly available as they were generated under a license for the current study. The data are, however, available from the author upon reasonable request and with permission of Shan Wang.