HES1 promoter activation dynamics reveal the plasticity, stemness and heterogeneity in neuroblastoma cancer stem cells

Neuroblastoma (NB) is the most common early pediatric solid tumor of the developing peripheral sympathetic nervous system (Gautier et al., 2021; Lumb and Schwarz, 2015; van Groningen et al., 2017). Previously, it has been shown that most NBs include two types of cells – neuroblastic (N-type) NB cells and substrate adhesive non neuronal (S-type) NB cells. Another group of cells, intermediate (I) type cells are widely considered as the stem-like cells in NB, and share morphological and biochemical features of both N and S type cells. I type cells show increased expression of NOTCH1 and other stem cell genes, such as PIGF2 (an alternative isoform of the PGF gene), GPRC5C (also known as RAIG-3), TRKB (also known as NTRK2) and p75NTR (also known as LNGFR, an isoform of NGFR), that maintains them in an under-differentiated state. This subpopulation of NB stem cells can transition to both an N type and S type cell fate (Ross et al., 1995, 2015; Veschi et al., 2019; Walton et al., 2004; Zeineldin et al., 2022). Later studies on the molecular characterization of NB cells revealed the existence of PROM1-negative adrenergic cells (ADR-type), which are comparable to N-type cells, and PROM1-positive mesenchymal cells (MES-type), which are comparable to S-type cells. Committed ADR-type cells and undifferentiated MES-type cells with two distinct gene expression profiles can interconvert upon activation of Notch and PRRX1 (van Groningen et al., 2017; Veschi et al., 2019).

Notch signaling is one of the highly conserved signaling pathways required for the maintenance and proliferation of embryonic and adult stem cells (Kopan and Ilagan, 2009; Misiorek et al., 2021). Any aberration in this pathway is often associated with the development of cancers (Aiello and Stanger, 2016; Karamboulas and Ailles, 2013; Matsui, 2016). Thus, canonical Notch signaling and its direct gene target HES1 play a critical role in regulating and maintaining cancer stem cells (CSCs), similar to embryonic development (Bolós et al., 2009; Zhang et al., 2021). Canonical Notch signaling involves proteolytic cleavages of receptor proteins that initiate the activation of Notch signaling. The binding of Notch ligands to the extracellular domain of Notch receptor (NECD) promotes the cleavage of the Notch intracellular domain (NICD) by γ-secretase followed by the translocation of NICD to the nucleus. There, it binds to the DNA-binding protein RBPjk (also known as RBPJ) or CSL [named for the orthologous mammalian, Drosophila and C. elegans proteins CBF-1 (RBPJ), Su(H) and Lag-1]. The interaction between NICD and RBPjk replaces the repressor complex bound with RBPjk with the activator complex and induces the expression of targeted genes of the HES family, p21 (CDKN1A) and MYC (Kopan and Ilagan, 2009; Misiorek et al., 2021).

Previously, we and others have demonstrated that the Notch downstream target Hes-1 can be activated non-canonically without interaction with the Notch receptor (Dave et al., 2011; Dhanesh et al., 2017; Ingram et al., 2008; Sanalkumar et al., 2010). Therefore, the expression of Hes-1 cannot be entirely blocked by treatment with the γ-secretase inhibitor DAPT (Sanalkumar et al., 2010; Xie et al., 2016). We investigated the requirement of differential Hes-1 activation in the mouse neocortex. We showed that non-canonical Notch-independent Hes-1 (NIHes-1) expression maintained neural stem cells in a slowly dividing state that later transited to the canonical Notch-dependent Hes-1 (NDHes-1)-expressing state. Notch receptor-mediated Hes1 activation was required to maintain neural progenitors and radial glial cells (RGCs) during neocortical development (Dhanesh et al., 2017). As CSCs are maintained through Notch signaling by upregulation of HES1 expression (Abel et al., 2014; Xiao et al., 2017), we hypothesized that there could be a similar scenario of differential HES1 activation that maintained subpopulations within the CSCs.

Here, using our reporter system, we captured the dynamic and heterogenous HES1 promoter activity within NB. We demonstrated that in IMR-32 NB cell line the Notch- and CBF-1-independent (i.e. NIHes-1) cells maintained a unique subpopulation of CSCs along with the NDHes-1-expressing CSCs. Interestingly, both these CSCs had the potential to switch between these two states. During this dynamic transition between NIHes-1 and NDHes-1, each type of CSC maintained its own population and mostly transited to a HES1 promoter inactive state designated as promoter inactive Hes-1 (PIHes-1). The PIHes-1 cells were subclassified based on the source of CSCs and designated as NI_PIHes-1 and ND_PIHes-1 cells. Both of these PIHes-1 cells were capable of reverting to the CSC state, which highlights how neuroblastoma cells possibly could evolve.

Based on our previous report of differential Hes-1 activation and maintenance in neural stem cells or progenitors, we decided to explore whether differential HES1 activation is present in CSCs of neural origin tumors, such as NB (Fig. 1A,B) (Dhanesh et al., 2017; Sanalkumar et al., 2010). To prove differential HES1 activation in NB, we electroporated the IMR-32 NB cell line with the pCBFRE-DsRedExpressDR-mtCBF-1-d2EGFP reporter plasmid (Fig. 1C). Here, the CBFRE-DsRedExpressDR cassette will report the NDHes-1 expression under the influence of 4× repeats of CBF-1-binding site, and the mtCBF-1-d2EGFP cassette with all the three CBF-1-binding sites of the HES1 promoter mutated will report NIHes-1 expression simultaneously (Dhanesh et al., 2017).

Our initial results showed the presence of both NIHes-1 and NDHes-1 expression in electroporated IMR-32 cell line (Fig. 1D,E). Here, the NIHes-1-expressing cells contributed only 1.1% of the total population whereas NDHes-1-expressing cells constituted 20.2% of the total population. In addition, we observed that ∼2.3% of the cells remained in a transition state or intermediate stage where they expressed both d2EGFP and DsRedExpressDR. The remaining 76.4% of cells had either no HES1 expression or were not electroporated with the reporter construct (Fig. 1F). Furthermore, we checked the HES1 expression between NIHes-1- and NDHes-1-expressing cells by quantitative real-time RT-PCR (qRT-PCR) analysis. The result showed that NDHes-1-expressing cells had nearly eight-fold the HES1 transcript expression when compared to NIHes-1-expressing cells (P<0.001, Fig. 1G). This result confirmed the presence of differential HES1 activation in neural origin tumors, similar to what is seen during neocortical development, and formed the basis for further investigating the role of NIHes-1 and NDHes-1 expression in NB. Based on the above findings, we asked the following questions: (1) Do the NIHes-1-expressing cells unidirectionally transit to an NDHes-1-expressing state similar to during neo-cortical development? (2) Do they have any CSC properties? (3) Does the differential HES1 activation have any functional implications? To answer these questions, a stable cell line was generated from IMR-32 NB cell line electroporated with the reporter plasmid. The cell subtypes were enriched based on d2EGFP or DsRedExpressDR expression through multiple FACS sorting and passaging as mentioned in the Materials and Methods, and used for all further experiments.

We next wanted to understand whether the NIHes-1-expressing cells are capable of a unidirectional transition to a NDHes-1-expressing state similar to that observed in neocortical development. For this, we sorted both NIHes-1- and NDHes-1-expressing cells by FACS, seeded them into a glass bottom dish and carried out time-lapse imaging of both the cells independently. We could observe that NIHes-1 cells that express d2EGFP transitioned to a NDHes-1 state with DsRedExpressDR expression and later divided into two daughter cells with NDHes-1 expression (Fig. 1H; Movie 1). These data indicate that it is the same cell that is capable of transitioning from the NIHes-1 state to the NDHes-1 state and there is not two different types of cell that utilize two differential modes of HES1 activation. We also observed NIHes-1-expressing cells dividing and generating two daughter NIHes-1 cells for maintaining its own population (Fig. 1I; Movie 2). Furthermore, and contradictory to the developmental scenario, we could observe a very few NDHes-1-expressing cells reverting to a NIHes-1-expressing state in vitro (Fig. 1J; Movie 3). We sorted NIHes-1- and NDHes-1-expressing cells by FACS and cultured them separately for 24 h and subjected the cells to FACS again to confirm our findings. NIHes-1 cells sorted by FACS after 24 h in culture transitioned to 0.1–0.2% (mean±s.d. 0.1%±0.1) of NDHes-1 cells and maintained 39.6–45.6% (mean 43.2%±2.73) of its own population and the remaining majority of cells had no fluorescence which was designated as NI_PIHes-1 cells (Fig. S1B). On the other hand, NDHes-1 cells sorted by FACS after 24 h of culture showed that 18.8–26.0% (mean 22.9±3.62%) of cells remained in a NDHes-1-expressing state itself. However, 0.1–0.9% (mean 0.2±0.39%) of cells were found in an intermediate state that could transit back to the NIHes-1-expressing state (Fig. S1C). Similar to what was seen with the NIHes-1-expressing cells, we also found the majority of cells without fluorescence in the NDHes-1 FACS analysis were ND_PIHes-1 cells (Fig. S1C). To rule out the possibility that the loss of fluorescence in PIHes-1 cells could be due to not having a copy of the reporter, we sorted ND_PIHes-1 cells derived from NDHes-1 cells by FACS and carried out PCR analysis against d2EGFP to reconfirm our stable cell line. Our data showed that indeed the reporter construct was present in the PIHes-1 cells (Fig. S1F).

Analyzing the FACS data, it appears that in NB, the minority of NIHes-1- and NDHes-1-expressing cells exist in a trans state where they can interconvert between NDHes-1- and NIHes-1-expressing states. During neocortical development, this is a unidirectional transition where the NIHes-1-expressing neural stem cells transition to the NDHes-1-expressing neural progenitors and RGCs and later differentiate into neurons by completely shutting down Hes-1 expression (Dhanesh et al., 2017). Given that we cannot expect the NB cells to be terminally differentiated, we cannot rule out the possibility of PIHes-1 population reverting to NIHes-1- or NDHes-1-expressing state in a reverse manner. We observed through time-lapse imaging that a very minimal number of NI_PIHes-1 cells retained the potential to retransform into a NIHes-1-expressing state spontaneously and further transitioned to a NDHes-1-expressing state (Fig. 1K; Movie 4). We did not observe any direct transition from the ND_PIHes-1 state to either of the NIHes-1- or NDHes-1 states through time-lapse imaging. To further confirm these findings, we sorted NI_PIHes-1 cells derived from NIHes-1 cells, and ND_PIHes-1 cells derived from NDHes-1 cells by FACS and cultured them for 24 h. These sorted NI_PIHes-1 cells after 24 h in culture transitioned to give 7.8–9.6% (mean 8.9±0.98%) in the NIHes-1-expressing state (Fig. S1D). By contrast, 0.2–0.6% (0.4±0.2%) of cells were observed to be in an intermediate state and 0.1–0.2% (mean 0.1±0.05%) of cells transitioned to the NDHes-1-expressing state. A majority of the NI_PIHes-1 population remained in the same state (Fig. S1D). However, FACS analysis of ND_PIHes-1 cells after 24 h culture showed that 94.6–99.9% (mean 98.4±2.59%) cells remained in ND_PIHes-1 state and a meager 0.1–5.3% (mean 1.5±2.54%) cells transitioned to the NDHes-1 state, which we failed to detect in time-lapse imaging (Fig. S1E). This result indicates that both NI_PIHes-1 and ND_PIHes-1 cells have the potential to revert to a CSC state in a tumor microenvironment.

Self-renewal potential is considered as the crucial hallmark of CSCs. We carried out a limited dilution hanging drop ‘tumorsphere’ formation assay to check whether the NIHes-1- or NDHes-1-expressing cells have CSC properties (Weiswald et al., 2015). For this, NIHes-1- or NDHes-1 cells sorted by FACS were serially diluted and hanging drops were made with 25,000, 2500, 250, 25 and 10 cells, respectively (Fig. 2A; Fig. S2G). Next, we checked for the number of sphere-like structures generated from each hanging drop. It was found that with increasing dilution, both NIHes-1- and NDHes-1-expressing cells started losing the potential to form primary spheroids (Fig. S2G). In order to understand the differential stem cell frequency of both NIHes-1 and NDHes-1 cells we performed extreme limiting dilution analysis (ELDA) (Hu and Smyth, 2009) and interestingly, we could find that stem cell frequency of NIHes-1-expressing cells (1/1783) was higher compared to NDHes-1-expressing cells (1/3570; Fig. 2B; Fig. S2G). Moreover, on the second day of culture, the NIHes-1-expressing cells retained the ability to form budding out secondary spheroids at a dilution of 250 cells (three secondary spheroids budded out from primary spheroids generated from 250 cells). By contrast, the NDHes-1 cells could generate secondary spheroids from primary spheroids until a 2500 cell dilution and did not generate any secondary spheroids from further dilutions (Fig. 2C,D; Fig. S2G).

Next, we analyzed the morphology of both NIHes-1- and NDHes-1-derived spheroids by measuring the sphericity of the primary spheroids. Even though both cell types could generate primary sphere-like structures, they were morphologically distinct. The sphericity was analyzed using AnaSP software, which generates a sphericity index between 0 and 1, where anything above 0.9 is considered highly spherical and values below 0.9 are considered non-spheroidal structures (Fig. 2E) (Zanoni et al., 2016). NIHes-1 cell-generated spheroids were more spherical (Fig. 2F,G), whereas NDHes-1-derived structures were not perfect spheres and showed a more flattening tendency although they were tightly compacted (Fig. 2H,I). The adhesive forces between cells are more significant in a hanging drop, allowing the cells to compact together. Still, we observed a reduced sphere-forming potential for the NDHes-1-expressing cells. The sphericity of NIHes-1-derived primary spheroids derived from all cell dilutions were between 0.9–1, whereas the NDHes-1 derived primary spheroids showed a value between 0.6–0.9 (Fig. 2J). We have also analyzed the self-renewal potential of both NI_PIHes-1 and ND_PIHes-1 cells using a sphere formation assay; these cells could form aggregates of cells but did not grow into spheroidal structures even after 7 days in culture (Fig. S2A–D). A similar observation was made with the unsorted IMR-32 cell line (Fig. S2E,F). These results demonstrate that the NIHes-1-expressing cells have comparatively higher stemness potential than NDHes-1-expressing cells, whereas both PIHes-1 cells lack stemness characteristics.

Furthermore, we checked for the expression of stemness markers by qRT-PCR analysis among all the cell types (i.e. NIHes-1, NDHes-1, NI_PIHes-1 and ND_PIHes-1 cells). The result showed a significant increase of stemness markers (i.e. NES, NANOG, ABCG2 and KIT) in NDHes-1-expressing cells compared to NIHes-1-expressing cells. However, we did not observe any change in expression of LGR5 between NIHes-1- or NDHes-1-expressing cells (Fig. 2K). NIHes-1 cells had properties of CSCs, such as being able to form spheroids, even though they had surprisingly low expression of common NB CSC markers.

Next, we compared the stem cell marker expression between NIHes-1- and NDHes-1-expressing cells and their respective PIHes-1 cells. For KIT and LGR5, we observed a significant downregulation in both PIHes-1 cells when compared to their respective CSCs. We did not find any significant change in expression of NES, NANOG and ABCG2 between NI_PIHes-1 and NIHes-1 cells, but a significant downregulation was observed in ND_PIHes-1 cells when compared to NDHes-1 cells. Thus, these findings support the reduced stemness observed in both PIHes-1 populations. Based on the above findings, we conclude that a small subpopulation of CSCs maintained by NIHes-1 expression exists within the pool of NDHes-1-maintained CSCs and differential expression of stemness markers is observed within these subclasses.

Next, we were intrigued to look at the expression of Notch components as it plays a crucial role in the regulation and maintenance of CSCs. qRT-PCR analysis showed a significant increase in expression of NOTCH1, NOTCH2 and NOTCH3 receptors, JAG1, DLL1 and DLL3 ligands, and RBPJ in NDHes-1 CSCs which suggested activation of classical Notch signaling in NDHes-1 CSCs (Fig. 2L). We also observed a reduction of Notch component expression in both the PIHes-1 cells compared to its respective NIHes-1- or NDHes-1 CSC state (Fig. 2L). Even though all Notch components were downregulated in NIHes-1 CSCs, we could still find the HES-1 expression in this subclass of CSCs as described in Fig. 2L. We found a reduction of NOTCH1 expression in NIHes-1 CSCs as compared to NDHes-1 CSCs, facilitating the non-canonical activation of HES1 possibly through FGF2 (Sanalkumar et al., 2010). Contrasting with the above findings, we found no difference in cleaved NICD expression between NIHes-1- or NDHes-1 CSCs (Fig. 2M). These results suggest that, in NIHes-1 CSCs, the cleaved NICD might not be recruited to the CBF-1-binding site of the HES-1 promoter as there was a reduction of RBPJ expression. This interpretation was made given the NIHes-1-expressing cells did not show any expression of DsRedExpressDR even though the reporter construct has four repeats of the CBF-1-binding site. Thus, a new subclass of CSCs, that is, NIHes-1 CSCs, was found to be maintained without the canonical Notch dependent HES1 expression.

Next, we determined the tumorigenic potential of both NIHes-1- and NDHes-1-expressing cells in vivo by serial xenotransplantation assay; for this, 10⁵ NIHes-1 or NDHes-1 CSCs sorted by FACS were subcutaneously injected into the right and left abdominal flanks of immunocompromised NOD-SCID mice, respectively (Fig. 3A–C). The primary tumor started forming from both NIHes-1- and NDHes-1-expressing CSCs by ∼30 days post-injection. Furthermore, the NIHes-1- and NDHes-1-derived primary tumors were dissected, dissociated and analyzed by FACS on the 60th day post-injection. Analysis data of NIHes-1-derived tumors showed the retention of 0.1% of NIHes-1-expressing cells, 0.1% cells that has transitioned to the NDHes-1 state, and that the majority of cells were NI_PIHes-1 or bulk cancer cells (Fig. 3E). Similarly, the NDHes-1-derived primary tumor analysis showed 1.3% of cells were in the NDHes-1-expressing state (Fig. 3F). We did not find any detectable NIHes-1-expressing cells from these tumors. However, we cannot rule out that NIHes-1-expressing are generated CSCs as we observed 0.1% of cells expressing both DsRedExpressDR and d2EGFP, indicating that these cells could be transiting to a NIHes-1-expressing state. Here too, we found that the majority of cells were PIHes-1 cells. We injected PIHes-1 cells subcutaneously into abdominal flanks of NOD-SCID mice, but we did not find any palpable tumor.

To further confirm whether the CSCs derived from these tumors can self-renew and retain their stemness and generate secondary tumors, we injected 4×10⁴ FACSorted NIHes-1- or NDHes-1 CSCs derived from primary tumors into the right and left abdominal flanks of NOD-SCID mice (Fig. 3D). Again, we observed tumor generation from both these CSC populations. Although we observed the NIHes-1 CSC-derived tumor to be small compared to the NDHes-1-expressing tumor by 45 days, we could not keep the mice any longer than 58 days since the NDHes-1 CSC-derived tumor was very vigorously growing and was nearing the permissible limit. Therefore, we dissected the NIHes-1- and NDHes-1-derived secondary tumors on the 58th day, dissociated them and subjected them to FACS analysis. We observed an enrichment of NIHes-1-expressing CSCs (17.4%) compared to the primary tumor (Fig. 3G). Moreover, we found 1.4% cells transitioned to the NDHes-1-expressing state and 1.9% cells were in the transition stage having both DsRedExpressDR and d2EGFP expression. Analysis of NDHes-1 CSC-derived secondary tumors showed 35% of cells were NDHes-1-expressing CSCs, whereas there was a negligible amount of NIHes-1-expressing CSCs and 0.3% of cells were in transition stage (Fig. 3H). Again, both the NIHes-1- or NDHes-1 CSC-derived secondary tumors generated the majority of the respective PIHes-1 cells.

Interpreting the data from primary and secondary tumors, we conclude that NIHes-1 CSCs indeed had lesser self-renewing potential as compared to NDHes-1 CSCs in each serial xenotransplantation assay. In addition to that, qRT-PCR analysis showed that NDHes-1 CSCs had higher expression of MYCN (Fig. S2H) and HES1 (Fig. 2L) along with other proliferation markers [i.e. CDK1, CDK2, CDK6, CDC42, CDKN1A (encoding p21) and CCND1 (encoding cyclin D1)] (Fig. S2H). We also checked the differential proliferation marker expression between NI_PIHes-1 and ND_PIHes-1 cells as these comprised the majority of the respective tumors. We found that CDK1, CDK2, CCND1, HES1 and MYCN were significantly upregulated in ND_PIHes-1 as compared to NI_PIHes-1 cells (Fig. S2H). In addition to that, MYCN overexpression is often associated with upregulation of NES, which mediates tumor aggressiveness (Thomas et al., 2004). Here, we observed a significant upregulation of NES expression in both NDHes-1 and ND_PIHes-1 cells compared to NIHes-1 and NI_PIHes-1 cells, respectively, which possibly explains the size variation of the secondary tumors derived from NIHes-1- or NDHes-1 CSCs (Fig. 2K).

We analyzed the tumor sections of secondary tumors derived from the NIHes-1- and NDHes-1 primary tumor. As expected, NIHes-1-expressing xenografts showed the presence of both d2EGFP and DsRedExpressDR-expressing cells, indicating that the NIHes-1 CSCs had transitioned to a NDHes-1 CSC state (Fig. 3I,J). However, we observed a limited potential of NDHes-1 CSCs to revert to the NIHes-1 CSC state as shown in Fig. 3H. While analyzing the FACS data from NDHes-1 secondary tumors, we observed 0.3% of cells in a transitioning state, which indicates the possible transition of NDHes-1 cells back to NIHes-1 CSCs in vivo. To confirm this possibility, we analyzed the NDHes-1 CSC-derived tumor sections and surprisingly found a very few patches of d2EGFP-expressing NIHes-1 CSCs in addition to NDHes-1 CSC with DsRedExpressDR expression (Fig. 3K,L). These results indicate that the NDHes-1 CSCs can revert to the NIHes-1 CSC state in vivo.

In order to further understand the transition from the NDHes-1 to the NIHes-1 state, we developed a 3D tumorsphere culture model with 2500 NDHes-1 CSCs sorted by FACS and cultured them for 7 days. Even though the sphericity of the NDHes-1 CSC-derived primary spheroids was below 0.9, we could see many budding secondary spheroids derived from the margins of NDHes-1 primary spheroids (Fig. 3M–Q). Our results showed that a few cells in the margin of these spheroids started expressing d2EGFP, indicating the transition to the NIHes-1 expression state (Fig. 3N,n′) and new secondary spheres with a sphericity index above 0.9 were formed containing both NIHes-1- and NDHes-1-expressing CSCs (Fig. 3O–Q). Therefore, we conclude that NDHes-1 CSCs are capable of reverting to the unique NIHes-1 CSC state in vivo. Our finding also throws light at the plasticity existing among the NIHes-1, NDHes-1 and PIHes-1 populations (Fig. 3R).

Given that NIHes-1- and NDHes-1-expressing cells displayed stemness potential both in vitro and in vivo, and NDHes-1 CSCs generated vigorous secondary tumor mostly contributed to by the bulk ND_PIHes-1 cells often associated with enhanced proliferation markers, we tried to unravel the transcriptomic signatures of these three cell types. We carried out differential gene expression analysis in two groups, that is, NIHes-1 versus NDHes-1 and ND_PIHes-1 versus NDHes-1. All three populations of cells were sorted by FACS and total RNA-seq was carried out, followed by gene set enrichment analysis (GSEA) using java GSEA 4.1 (Fig. 4A). When comparing NIHes-1 CSCs with NDHes-1 CSCs with the threshold of adjusted P<0.05, we found differentially expressed genes as shown in the volcano plot (Fig. 4B; Table S3). For the GSEA analysis, we ranked the genes based on fold change direction and adjusted P-value. Ontology gene sets from the Molecular Signature Database (MSigDB; https://www.gsea-msigdb.org/gsea/msigdb/) were used and GSEA was performed. It showed that NIHes-1 CSCs were enriched with terms related to the Wnt signaling pathway, apoptosis/autophagy signaling mechanisms, and molecules relevant for cell adhesion (Fig. 4C,D; Fig. S3, Table S4). Interestingly, we did not find the Notch signaling pathway to be enriched in NIHes-1 CSCs. On the other hand, NDHes-1 CSCs were enriched with the Notch signaling pathway as found in the GSEA (Fig. 4C,E; Fig. S3, Table S5). These cells were also found to be enriched with epithelial to mesenchymal (EMT) transition markers along with stem cell differentiation marker genes (Fig. 4C,E; Fig. S3). Moreover, the NDHes-1 CSCs were enriched with neural crest cell (NCC) differentiation markers, which reinforces the NCC origin of the NB cell line used here, that is, IMR-32 cells (Fig. 4C,E). The combined GSEA plots in Fig. 4D,E show the positions for each marker genes for the pathways in the ranked gene list. Wnt signaling was enriched in NIHes-1 CSCs and Notch signaling pathway was enriched in NDHes-1 CSCs (Fig. 4D,E). This validates our transcriptome data and the reporter system that was used. While comparing with NDHes-1 CSCs and ND_PIHes-1 cells, we found 118 genes differentially expressed above the threshold mentioned above (Fig. 4F; Table S6). Although we did not find any ontology gene set enrichment, we found ND_PIHes-1 cells had the signature expression of developing fibroblast-like cells (Fig. 4G; Fig. S3, Table S7). This provided us a clue about the fibroblast nature of the ND_PIHes-1 cells, which had upregulated expression of fibulin family matrix protein FBLN2, collagen molecule COL6A2 and matrix remodeling molecule MXRA8 (Fig. S3). To identify the cell type markers, we carried out GSEA using cell-type signature gene sets curated from publicly available single-cell sequencing studies of human tissues (Fig. 5A). We found that neural stem cell markers, progenitor cells markers and microglial cell markers were enriched in NIHes-1 CSCs, whereas NDHes-1 CSCs showed significant expression of adrenal sympathoblast marker genes (Fig. 5B,C; Tables S8, S9).

We observed significant differential HES1 expression between NIHes-1- and NDHes-1-expressing CSCs, as shown above (Fig. 1G). However, we also observed a sustained low HES1 expression in the ND_PIHes-1 cells as compared to NDHes-1 CSCs (Fig. 2L). Although the HES1 promoter is not active in ND_PIHes-1 cells, the sustained HES1 expression led us to investigate HES1 mRNA stability, possibly regulated by its 3′ untranslated region (UTR). Previous studies have shown that longer 3′ UTR isoforms are more prevalent in quiescent fibroblast cells and are associated with increased expression and mRNA stabilization (Mitra et al., 2018). It has also been shown in mature hippocampal neurons that the neuropil-enriched transcripts tend to have significantly longer 3′ UTR. Surprisingly, those transcripts have a longer half-life than the shorter isoforms (Tushev et al., 2018). Quantitative alternative poly-adenylation (QAPA) analysis with our bulk transcriptome data showed a differential 3′ UTR isoform of HES1 transcript usage among the three cell types (Fig. 5D). It was found that all the three cell types had the highest expression (nearly 80%–90% of total isoform usage) for the shorter 3′ UTR isoform of 387 nucleotide length. At the same time, there was a shift towards the use of the distal poly-adenylation site generating a longer isoform (499 bp) in the ND_PIHes-1 cells, which accounted for nearly 11% of total isoform usage. This longer 3′ UTR isoform was not found in the other two cell types. There was another variant of 3′ UTR having 454 bp length, which was only found in NIHes-1- and NDHes-1-expressing CSCs and accounted for nearly 7% of transcripts (Fig. 5D). The presence of longer 3′ UTR in ND_PIHes-1 cells might be associated with stabilization of the HES1 transcript.

Previous findings have suggested that NBs transcriptionally match with fetal adrenal neuroblasts when comparing the single-cell transcriptomes of NBs and normal developing adrenal glands (Jansky et al., 2021). Here, we wanted to understand the proportion of relevant developmental cell types in three different bulk RNA-seq samples of IMR-32 cell populations (i.e. NIHes-1, NDHes-1 CSCs and ND_PIHes-1 cells) and predict its functional implication and the stemness potential. We used the deep learning-based single-cell assisted deconvolution network (SCADEN) method to achieve that (https://github.com/KevinMenden/scaden). We used the published single-cell transcriptome data (Fig. 5E) of fetal adrenal medulla to train the model and later, it was used to predict the proportion of cell types present in each of the three cell populations that we studied. It was found that cycling neuroblast like signatures (i.e. TOP2A and MKI67), which are associated with high proliferative capacity were enriched in NIHes-1 CSCs whereas neuroblast-like signature genes, that is NEFM, STMN2, GAP43 and ISL1 were highly expressed in NDHes1 and ND_PIHes-1 cells (Fig. 5F; Fig. S3). Markers of increased differentiation, like SYN3, were enriched in NDHes-1 CSCs and ND_PIHes-1 cells. In the GSEA, we found that NDHes-1 CSCs were enriched with stem cell differentiation markers and we did not observe any difference in the levels of these markers in ND_PIHes-1 cells (Fig. 4C). Thus, these results show that NIHes-1-expressing CSCs have reduced expression of differentiation markers. Late-stage developmental cell types, such as late chromaffin, late Schwan cell precursors (SCPs) and late neuroblast signatures, were proportionately low in all three cell types (Fig. 5F). High expression of SOX10 in conjugation with MKI67 and ERBB3 expression in ND_PIHes-1 cells closely resembles the expression seen in cycling SCPs (Fig. 5F; Fig. S3). Thus, these results suggest that the reporter construct distinctly separates three different populations of cells, where NIHes-1 and NDHes-1 CSCs resemble cycling neuroblasts and neuroblasts, respectively. On the other hand, ND_PIHes-1 cell types had mixed features of both cycling SCPs and neuroblasts. Transcriptomic signature matching with cycling SCP cells supports the high proliferative ability of ND_PIHes-1 cells. We cannot expect the NB cells to be terminally differentiated even though they had increased differentiation markers.

We were also curious to know whether the NIHes-1 CSCs were responsible for metastasis. To assess the ability of both NIHes-1 and NDHes-1 CSCs to undergo distant metastasis, we sorted 10⁴ NIHes-1- or NDHes-1 CSCs by FACS and seeded them into a two-chambered migration assay plate and treated them with a non-toxic dose of mitomycin-C to inhibit cell division. We found substantial migration in NDHes-1 CSCs compared with NIHes-1 CSCs from 6 h post wound formation (Fig. 6A,B). Given that our transcriptome data showed that the NDHes-1 CSCs had enriched levels of EMT markers (Fig. 6D), they would be expected to undergo EMT and that could be the reason for the observed migration. Moreover, the NIHes-1 CSCs had a significantly higher expression of cell adhesion molecules (Fig. 6C), which could be preventing migration. As cell migration is a critical factor for metastasis, we also analyzed the metastatic potential of NIHes-1- or NDHes-1 CSCs in vivo. For this, we generated a new stable IMR-32 cell line with constitutively active luciferase (Fig. 6E; Fig. S4). This cell line already had stable expression of the reporter construct depicted in Fig. 1C to allow it to be sorted by FACS. 2×10⁴ NIHes-1- or NDHes-1 CSCs were injected through the tail vein in separate groups of NOD-SCID mice (n=3). The animals were imaged for luciferase activity in a whole-body imager on the 22nd and 28th day after injection. We observed the tumor formation in brain and spinal cord of the group injected with NDHes-1 CSCs, but NIHes-1-injected group did not show any metastasis, which is in line with our in vitro migration assay results (Fig. 6F–I; Fig. S5).

Therefore, our results demonstrate that the CSCs present in NB are heterogeneous regarding their stemness, plasticity and metastatic potential. Given that PIHes-1 cells did not generate any tumor in vivo and did not form any sphere in vitro, we did not conduct an in vivo metastatic assay for this cell type.

From our data, we demonstrate different modes of HES1 expression in cells classified as NIHes-1, NDHes-1 CSCs, and respective NI_PIHes-1 and ND_PIHes-1 cells. Furthermore, we depict a coordinated bidirectional transition between the NIHes-1-expressing state and NDHes-1-expressing state. In addition to that, both NIHes-1- and NDHes-1-expressing cells divide to maintain their own population. Our data support previous findings that demonstrated that CSCs ‘empower’ themselves to undergo reversible transitions in a microenvironment with continued pressures (Easwaran et al., 2014; Knoechel et al., 2014; Liau et al., 2017). Heterogeneity in NB is well explained by the interconversions that were observed between the cellular phenotypes in conjunction with their molecular signatures (Chakrabarti et al., 2012; Veschi et al., 2019). We could find an adaptive plasticity in the use of HES1 promoter among NB CSCs. The slow transition that occurred in primary and secondary tumors, in contrast to the rapid transition that we observed in in vitro culture condition could be due to the change in cell–cell interaction in tumors and 3D cultures as compared to monolayer cultures. Therefore, the rapid dynamics observed in monolayer culture needs to be explored further. We also assessed the differential HES1 promoter activation in mouse NB cell line Neuro2A (N2A) as well as human glioblastoma cell line (U87MG) to understand whether the differential modes of HES1 expression occur in IMR-32 cells alone. The reporter construct was transfected to both the cells and it was found that both NIHes-1- and NDHes-1-expressing cell types existed in each cell line indicating that our findings are not limited to IMR-32 cells (Fig. S6).

We found that both NIHes-1- and NDHes-1-expressing cells give rise to a promoter inactive state of the HES1 gene (i.e. NI_PIHes-1 and ND_PIHes-1 cells; Fig. S1B,C). The ND_PIHes-1 cells were similar to supporting fibroblast cells and highly proliferative in nature, as found by a SCADEN analysis. The most notable finding here was that a few NI_PIHes-1 and ND_PIHes-1 cells preserved an ability to reprogram, where there was a transition from promoter-inactive state of HES1 to the respective CSC state (Fig. 1K) (van Groningen et al., 2019). This plasticity observed between the NIHes-1- and NDHes-1 CSCs and PIHes-1 cells might be the reason for the frequent relapse observed with neural origin cancers, which makes the treatment strategies even more challenging (Basta et al., 2016; Eleveld et al., 2015; Tomolonis et al., 2018). Moreover, we could not find NIHes-1- or NDHes-1-expressing cells that divided into two daughter cells where one had an NDHes-1-expressing state and another with NIHes-1-expressing state. Instead, in our time-lapse imaging, observed a coordinated transition of NIHes-1-expressing cells to NDHes-1-expressing cells and vice versa both in vitro and in vivo.

CSCs are known to have the ability to reconstitute tumors in vivo. Among the CSCs, a subtle line has been drawn between tumor-perpetuating cells (TPCs) and tumor progenitor cells (TProgs), where the TProg cells have the significant proliferative capacity necessary to generate tumors but cannot recapitulate the entire cellular heterogeneity of a tumor (Williams et al., 2013). It is important to note that simple tumorigenic potential does not indicate that cells have a CSC property. Rather loss of heterogeneity in the resulting tumor must be assessed in order to define the CSC property (increased or decreased stemness) of a subpopulation. Moreover, serial transplantation assays must be performed to check whether self-renewing capacity is inherent or lost (Tomolonis et al., 2018). Among the cell subpopulations, NIHes-1 cells showed enhanced stemness compared to NDHes-1-expressing cells in terms of 3D tumorsphere formation ability. The enhanced stemness observed with NIHes-1 cells has to be further explored using a panel of stem cell markers. Furthermore, in the xenotransplantation model, we found that both NIHes-1- and NDHes-1-expressing cells formed primary and secondary tumors. The NDHes-1 CSCs had enhanced self-renewing potential compared to NIHes-1-expressing cells, as revealed by the primary and secondary tumor FACS analysis. Both CSCs gave rise to the majority of non-tumorigenic PIHes-1 cells in both primary and secondary tumors. Thus, both the CSCs had self-renewing ability and potential to generate Notch-related heterogeneity. Although we have not ignored the bidirectional transition during tumor formation, this might have obscured the simplistic comparison between the NIHes-1- or NDHes-1 CSCs-derived tumors.

From the bulk RNA-seq data, we found that mRNA encoding the extracellular matrix protein tenascin c (TNC) was upregulated in NIHes-1-expressing cells, which might support the cell cycle by increasing Wnt signals (Fig. S3). The strongest signal that suggests the stemness of NIHes-1-expressing cells was the expression of LGR5 and downregulation of PMP22 (Fig. S3). LGR5 encodes a G-protein-coupled receptor (GPCR) that acts as a downstream target of Wnt signaling, preventing the differentiation of stem cells (Yang et al., 2020). On the other hand, the downregulation of PMP22, a membrane glycoprotein, restricts stem cell differentiation (Yang et al., 2020). Interestingly, CSCs are generally known for their high intrinsic resistance to undergo cell death through an apoptotic signaling mechanism (Fulda, 2013). Nevertheless, the unique NIHes-1 CSC subpopulation was enriched with transcriptomic signatures of apoptosis and autophagy signaling mechanisms. This distinguishes the unique subpopulation of NIHes-1 CSC from the already reported CSC features. Previous studies have shown that activation of the Notch signaling pathway can promote cell proliferation and metastasis (Hu et al., 2012; Yuan et al., 2015). Aberrant upregulation of BMP4, which we observed in NDHes-1 CSCs, possibly increased Jagged-1 expression, which regulates cell migration and invasion (Yang et al., 2020).

Furthermore, we have found differential HES1 expression patterns between NIHes-1 and NDHes-1 CSCs. However, a sustained low HES1 expression in ND_PIHes-1 cells intrigued us, and we hence examined HES1 transcript stability as regulated by its 3′ UTR region. Currently, we do not know the functional implication of differential 3′ UTR usage among these three cell types. To our knowledge, this is first time that such alternate 3′ UTR usage of HES1 transcript has been reported in CSCs. We need to investigate further the role of alternative 3′ UTR isoforms of HES1 in maintaining CSCs. We also assume that there could be a differential expression of alternative splice variants of many transcripts among these cell types that might define their characteristics and role in the tumor microenvironment.

In summary, we conclude that there are two different subpopulations of cells within the CSC pool examined here. The varying modes of HES1 promoter activation and interplay between these modes of activation possibly reports the Notch-related heterogeneity and plasticity among the CSCs. In addition, the fibroblast-like PIHes-1 cells were also capable of reverting to a CSC state, thereby demonstrating robust plasticity between all the cell types within the NB. This scenario also points towards the importance of eliminating the fibroblast-like cells in addition to the CSCs to ultimately reduce the chance of relapse and metastasis. Most importantly, this model can further be utilized for drug screening purposes, designing therapeutic targets and extrapolation to other cancer scenarios.

Three cell lines were used in the study. The human NB cell line IMR-32 (RCB1895) was procured from Riken and maintained in Dulbecco's modified Eagle's medium (DMEM) (Gibco, 11960-044) supplemented with 1% streptomycin and penicillin (Gibco, 15140-122), 10% FBS (PAN, P30-3302), 1% Glutamax (Gibco, 35050), 1% non-essential amino acids (NEAA) (Gibco, 11140-050). The mouse NB cell line (N2A) was maintained in DMEM supplemented with 1% streptomycin and penicillin, 10% FBS, 1% Glutamax, 1% NEAA and 1% sodium pyruvate (Gibco, 11360-070). The human glioblastoma cell line U87MG was maintained in DMEM supplemented with 1% streptomycin and penicillin, 10% FBS, 1% Glutamax and 1% NEAA. All cells were maintained at 37°C with 5% CO₂ for all further experiments. The cell lines were tested for mycoplasma and other contamination and had not been authenticated recently.

The pGL4.51 LUC/Puro vector was constructed by replacing SV40 neomycin resistance of pGL4.51 (Promega) with puromycin under the PGK promoter from pL552 (Chen et al., 2015; Addgene #68407; Fig. S4). pCBFRE-DsRed Express-DR-mtCBF-1-d2EGFP (∼9 kb) was generated previously by our group as described in Dhanesh et al. (2017). CBFRE-DsRedExpressDR was generated replacing luciferase with DsRedExpressDR in a 4XwtCBF1-Luc construct which was a kind gift from Diane Hayward, Johns Hopkins University, Baltimore, USA (Shawber et al., 1996). The mtCBF1 construct accurately reflected the NIHes-1 expression, which was earlier proven by testing using the Notch signaling blocker DAPT and further it was shown to express independently of NICD expression (Dhanesh et al., 2017).

IMR-32 cells were trypsinized separately with 0.05% Trypsin-EDTA (Gibco, 25300-054) and pCBFRE-DsRedExpressDR-mtCBF-1-d2EGFP plasmid was electroporated into the dissociated cells using a Neon Electroporator (Invitrogen). The pCBFRE-DsRedExpressDR-mtCBF-1-d2EGFP plasmid was previously generated in our lab and can simultaneously track NIHes-1 expression with d2EGFP and NDHes-1 expression with DsRedExpressDR (Dhanesh et al., 2017). Four pulses of 1400 V and 20 ms pulse width were given using the Neon electroporator (Invitrogen) and the cells were cultured in DMEM at 37°C with 5% CO₂. After multiple rounds of passaging, the cells were sorted by FACS, and cells with stable integration of plasmid were expanded and maintained. All further experiments were carried out with this stable cell line. Mouse NB N2A cell line and human glioblastoma U87MG cell line were also transfected with pCBFRE-DsRedExpressDR-mtCBF-1-d2EGFP plasmid using the manufacturer's lipid-based transfection protocol (Invitrogen, Lipofectamine 2000 transfection reagent).

The CBFRE-DsRedExpressDR-mtCBF-1-d2EGFP-LUC/Puro IMR-32 cell line was generated by electroporating pGL4.51 Luc/Puro (Fig. S4) into CBFRE-DsRedExpressDR-mtCBF-1-d2EGFP-IMR-32 cells using the conditions mentioned above. A stably integrated cell line was generated after multiple rounds of passaging followed by puromycin (1.5 µg/ml; Sigma, P7255) selection. The selected cells were expanded for in vivo tail vein injection experiments.

FACS analysis and sorting were carried out using a standard protocol (Hendon-Dunn et al., 2016). NIHes-1-expressing d2EGFP and NDHes-1-expressing DsRedExpressDR CSCs sorted by FACS were used for all experiments. FACS was performed using BD FACSaria III System. 75% confluent CBFRE-DsRedExpressDR-mtCBF-1-d2EGFP-IMR-32 cells were trypsinized using 0.05% trypsin-EDTA and passed through a 40 mm cell strainer to avoid clumping. Then, ∼10,000 cells were analyzed in BD FACS Aria III sorter by exciting the cells with a 488 nm laser. The green fluorescence emission was detected in the 530/30 filter with 502 long pass (LP) and red fluorescence emission with 585/42 filter with 556LP. To set a population gate, non-electroporated IMR-32 cells were first passed through the flow cytometer and analyzed using forward scatter versus side scatter. The voltage for FITC and PE channels were adjusted such that the control non-fluorescent cells will appear within the first order of the logarithmic scale of fluorescence in the fluorescence histogram. Next, CBFRE-DsRedExpressDR-mtCBF1-d2EGFP-IMR-32 cells were analyzed with these settings. ∼10,000 events were analyzed using Diva Software 8.0.3. Fluorescent cells were gated using fluorescence dot plots and sorted. In order to check the efficacy of sorting by FACS, cultures containing both NIHes-1 (d2EGFP) and NDHes-1 IMR-32 cells (DsRedExpressDR) were sorted into d2EGFP-expressing NIHes-1 cells and DsRedExpressDR-expressing NDHes-1 cells. The sorted cells were reanalyzed 0 h post the sorting to confirm the efficacy of sorting that was always maintained at above 99%. NIHes-1, NDHes-1, NI_PIHes-1 and ND_PIHes-1 cells were sorted by FACS and were collected in 50% FBS and expanded/used for transcriptome analysis or other experiments. Non-electroporated IMR-32 cells were used as the control for all experiments (Hendon-Dunn et al., 2016).

NIHes-1 and NDHes-1 cells sorted by FACS were cultured in a 35 mm glass-bottom dish. The cells were maintained at 37°C with 5% CO₂ in a microscope-mounted CO₂ incubator attached to the Olympus OX-11 microscope. Fluorescence images were captured continuously at the same frame in a 30 min interval for 14–24 h and analyzed using CellSense software (Olympus).

To carry out 3D tumorsphere formation assay, NIHes-1, NDHes-1, NI_PIHes-1 and ND_PIHes-1 cells sorted by FACS were serially diluted to obtain a 25,000, 2500, 250, 25 and 10 cell dilution in 30 µl of IMR-32 culture medium (Weiswald et al., 2015). Hanging drops (10–18) were made with 30 µl of IMR-32 medium containing the desired number of cells in a non-adherent Petri dish lid and allowed to grow for 24 h as hanging drops in the inverted dish lid. The desired number of cells accumulates at the accessible liquid–air interface to form tumorspheres. Adhesive forces between cells are more significant in a hanging drop, allowing the cells to compact together. The spheroids were counted using an inverted phase contrast microscope. Data were analyzed using extreme limiting dilution analysis (ELDA) by using a modified script from the ‘statmod’ R package (https://cran.r-project.org/web/packages/statmod/index.html) to compare the stem cell frequencies and their statistical significance (Hu and Smyth, 2009). After 24 h, the spheres were transferred to poly-2-hydroxyethyl methacrylate (HEMA; Sigma, P-3932)-coated Petri plates to reduce attachment of the spheres onto Petri plates. The spheroids were cultured in IMR-32 medium at 37°C with 5% CO₂ for 7 days post-seeding on poly-HEMA-coated Petri plates. Images were captured using an Olympus OX-11 microscope each day for 7 days of the culture period. Inverted Mask files of brightfield images were generated using ImageJ software and data extracted for multiple morphological parameters, including sphericity, using AnaSP tool (version 1.4), an open-source software (Piccinini, 2015). The sphericity index was demarcated based on the score, where a value between 0.9 and 1 is considered a compact regular classical sphere (Zanoni et al., 2016), whereas a value below 0.9 is considered as a non-spheroidal structure.

NIHes-1, NDHes-1, NI_PIHes-1 and ND_PIHes-1 cells sorted by FACS were used for total RNA isolation using the Qiagen RNeasy Mini kit (cat. no. 74104) according to the manufacturer's protocol. 1 µg of total RNA was converted to cDNA using random hexamers (Promega) and Superscript RT-II (Invitrogen). Relative expression of stem cell markers (i.e. NESTIN, NANOG, ABCG2, KIT and LGR5), Notch signaling components (NOTCH1, NOTCH2, NOTCH3, DLL1, DLL3, JAG1, RBPJ and HES1) and proliferation markers (CDK1, CDK2, CDK6, CDC42, CDKN1A, CCND1 and MYCN) were measured using SYBR green premix (Takara) along with mRNA encoding β-actin as an internal control (Table S1). Real-time analyses were done with a 2−ΔΔct method in a QuantStudio 3 (Applied Biosystems) real-time PCR machine (Dhanesh et al., 2017).

NIHes-1 and NDHes-1 CSCs sorted by FACS were lysed with RIPA buffer and total protein concentration was quantified by BCA protein estimation assay. ∼50 µg of protein was loaded onto PAGE gel, and immunoblotting was carried out with cleaved NICD and β-actin (internal control) antibodies (Table S2).

Xenotransplantation experiments were carried with approval from the Institutional Animal Ethics Committee (IAEC) of Rajiv Gandhi Center for Biotechnology (IAEC/626/JAC/2017). The animals were housed in IVC cages under standard temperature, humidity, light cycle and were provided with standard feed and water ad libitium. Tumorigenicity was analyzed by ectopic injection of 10⁵ NIHes-1 and NDHes-1 cells, sorted by FACS, respectively. Both of these cell types were re-suspended in 50 µl of sterile 1× PBS and mixed with Matrigel (20 mg/ml) in a 1:1 ratio (Mollo et al., 2016). The whole 50 µl of cell suspension along with Matrigel (Corning, 356231) was injected into the left (NDHes-1) and right (NIHes-1) subcutaneous flanks, respectively of 6-week-old male NOD-SCID mice (The Jackson Laboratory; RRID:IMSR_JAX:001303) in two separate experimental groups. Tumor formation was assessed regularly by looking for a palpable tumor. The same protocol was followed to inject PIHes-1 cells. The body weight of animals and the size of tumors formed were continuously monitored. The animals were euthanized if the body weight was reduced to more than 15–20% of initial body weight or the tumor size was more than 1.5 cm.

For carrying out migration assay, NIHes-1 and NDHes-1 CSCs sorted by FACS were seeded separately on a 35 mm culture-insert with a 2-well µ-dish (Ibidi 81176). This culture–insert contains a removable silicon gasket with two 70 µl wells on which cells can be seeded as separate chambers. Once the cells reached ∼90% confluence, they were treated with 12.5 µg per ml of mitomycin-C (Sigma, M0503) for 4 hours to inhibit cell proliferation. Following treatment, cells were washed to remove mitomycin-C entirely and then the insert was gently removed, leaving a perfect linear scratch. Images were captured at 0, 6, 24 and 32 h using an Olympus OX-11 microscope and a minimum of 15 arbitrary measurements were taken and analyzed (Bobadilla et al., 2019; Jonkman et al., 2014; Justus et al., 2014). Wound closure percentage was calculated using the formula: wound closure percentage=100−{A(t)/Mean A(0)*100%} where, Â(t) is wound area percentage, A(t) is wound area at the time and A(0) is initial wound area.

The in vivo metastatic potential of NIHes-1 and NDHes-1 CSCs were analyzed with tail vein injection of cells. Tail vein injection was carried out according to the Institutional Animal Ethics Committee guidelines (IAEC/685/JAC/2018). ∼75% confluent CBFRE-DsRedExpressDR-mtCBF-1-d2EGFP-LUC/Puro IMR-32 cell line were trypsinized using 0.05% trypsin-EDTA and passed through a 40 mm cell strainer to avoid clumping and sorted by FACS for d2EGFP and DsRedExpressDR-positive cells. 2×10⁴ NIHes-1 and NDHes-1 CSCs sorted by FACS were reconstituted in ∼100 µl of 1× PBS and ∼100 µl of each cell type was injected into the tail vein of 6-week-old male NOD-SCID mice in two separate experimental groups. For in vivo imaging, the animals were anesthetized with 3% isoflurane and maintained with 2% isoflurane during imaging. Firefly luciferase activity was analyzed after injecting 2 mg luciferin/150 µl of 1× PBS. Bioluminescence was measured using an in vivo whole body luminescence imager (Perkin Elmer, CA, USA; IVIS Spectrum in vivo imaging system) and analyzed using Living image version 4.5.2.18424 software. One group of NOD-SCID mice was imaged without injecting luciferin and was used as a control.

Three different subpopulations of cells were sorted by FACS based on reporter fluorescence, and total RNA was isolated using a Qiagen RNA isolation kit (cat. no. 74104). The total RNA was treated with DNaseI to remove any DNA contamination. RNA was later quantified prior to ribo-zero library preparation. Paired-end RNA-seq was performed on an Illumina HiSeq platform. The reads were aligned to the human reference genome using STAR aligner. The human reference genome was downloaded from UCSC genome server (https://github.com/morrislab/qapa). After the alignment, the count matrix was obtained from the biological replicates using ‘featureCounts’ program of the subread package. Thereafter differential gene expression analysis was performed using DESeq2 package in R. Two differential expression analyses were performed, NDHes-1 versus NIHes-1-expressing CSCs and NDHes-1 versus ND_PIHes-1 cells. Volcano plots were made to illustrate the differentially expressed genes using ‘ggplot2’ package in R. Heatmaps were generated using the ‘pheatmap’ package. Gene Set Enrichment Analysis (GSEA) was performed using java GSEA 4.1. GSEA was performed in pre-ranked mode as described previously (Reimand et al., 2019; Subramanian et al., 2005). Two different gene sets, ontology gene sets and cell-type signature gene sets, were used for the analysis, downloaded from the Molecular Signature Database (MSigDB).

Single-cell assisted deconvolution network (SCADEN) analysis was carried out to perform cell composition analysis from bulk RNA-seq data (Menden et al., 2020). We have used scRNA-seq data of fetal adrenal medulla taken from published scRNA-seq data (European Genome-Phenome Archive accession no. EGAS00001004388; Jansky et al., 2021). The adrenal medulla seurat object was downloaded from https://adrenal.kitz-heidelberg.de/developmental_programs_NB_viz/ and the count matrix was obtained by using ‘seurat’ package. 100 artificial RNA-seq samples were simulated by subsampling the default number of cells from the input scRNA-seq data sets for training data. Later, our bulk RNA-seq data with raw counts were used to predict the cell-type proportion present in each of the bulk RNA-seq data.

Alternative poly-adenylation usage was calculated using a quantitative method called quantitative alternative polyadenylation (QAPA) (Ha et al., 2018). The 3′ UTR library was first downloaded from the Morris lab github repository (https://github.com/morrislab/qapa). 3′ UTR sequences were extracted using the human reference genome in fasta format. Then 3′ UTR isoform usage was quantified using the transcript quantification tool ‘salmon’. Quantified files from each sample were merged using create_merged_table R script and then the relative proportion of each isoform in a gene, measured as poly(A) usage (PAU), was calculated using the compute_pau.R script.

All the statistical tests were performed using GraphPad Prism 9. As mentioned in the figure legends, the n values represent the biological repeats measured independently. We calculated the significance by performing a non-parametric, unpaired, two-tailed Mann–Whitney test. The significance for qRT-PCR data was measured by performing parametric unpaired two-tailed multiple t-tests. Significance was defined at P<0.05. Analysis software other than Prism includes ImageJ, AnaSP, R (v4.1.0) and Microsoft Excel.

We thank Dr SanalKumar RS and Dr Nicolo Riggi, Institute Universitaire de Pathologie, Rue du Bugnon 25, 1011 Lausanne, Switzerland for creative suggestions regarding data mining and transcriptome analysis. We thank Ms Indu Ramachandran, Ms Arya Venukumar Sreelekha, Mr Tilak Prasad, Ms Surabhi Subramoniam Vimala and Ms Tanima Chinnu Tomi for their help in performing FACS analysis. We also thank AgriGenome Pvt. Ltd. (Cochin) for carrying out and providing us with the RNA sequencing data and Mr Biju Surendran Nair for his help with animal experiments. The Bioinformatics facility of Rajiv Gandhi Centre for Biotechnology (RGCB) is gratefully acknowledged for giving access to the high-performance computing facility for carrying out transcriptome data analysis. We thank Dr Rakesh S. Laishram for providing primers for proliferation marker genes and Dr Vazhanthodi Abdul Rasheed for his help in generating schematic representations.

The peer review history is available online at https://journals.biologists.com/jcs/lookup/doi/10.1242/jcs.260157.reviewer-comments.pdf.