Neural crest-related NXPH1/α-NRXN signaling opposes neuroblastoma malignancy by inhibiting metastasis

Neuroblastoma is a pediatric cancer that can present as low- or high-risk tumors (LR-NBs and HR-NBs), the latter group showing poor prognosis due to metastasis and strong resistance to current therapy. NBs are known to originate from alterations to cells in the sympatho-adrenal lineage derived from the neural crest, but whether LR-NBs and HR-NBs differ in the way they exploit the transcriptional program underlying their developmental origin remains unclear. Here, we compared the transcriptional landscapes of primary samples of LR-NBs, HR-NBs and human fetal adrenal gland, and thereby identified the transcriptional signature associated to NB formation that further distinguishes LR-NBs from HR-NBs. The majority of the genes comprising this signature belong to the core sympatho-adrenal developmental program, are associated with favorable patient prognosis and with diminished disease progression. The top candidate gene of this list, Neurexophilin-1 (NXPH1), encodes a ligand of the transmembrane receptors α-Neurexins (α-NRXNs). Our functional in vivo and in vitro assays reveal that NXPH1/α-NRXN signaling has a dual impact on NB behavior: whereas NXPH1 and α-NRXN1 promote NB tumor growth by stimulating cell proliferation, they conversely inhibit the ability of NB cells to form metastases. Our findings uncover a module of the neural crest-derived sympatho-adrenal developmental program that opposes neuroblastoma malignancy by impeding metastasis, and pinpoint NXPH1/α-NRXN signaling as a promising target to treat HR-NBs.


Introduction
Neuroblastoma (NB) represents the most common cancer in infants and the most frequent solid extracranial tumor in childhood, accounting for 10-15% of cancer-related child deaths (1)(2)(3). NBs are nearly all sporadic (98%) and are very heterogeneous (4), half of the tumors being detected in para-spinal ganglia in the abdomen, chest and neck, while the other half are found in the adrenal gland medulla (5). Some 60% of NBs are classified as low-to-intermediate risk (LR-NBs) and correspond mostly to loco-regional tumors that present a good prognosis and can even regress spontaneously (6,7). By contrast, the remaining 40% are classified as high-risk tumors (HR-NBs) and are often already metastatic at the time of diagnosis, displaying a strong resistance to therapy and a high probability of relapse that is reflected in poor patient outcome (8,9). The most predictive genetic events of a poor prognosis for NB patients appear to be amplification of the proto-oncogene MYCN, and mutations in the ALK, ATRX and PHOX2B genes (6,8,10). However, these genetic alterations are only found in up to 16% of all NB cases and thus, they are insufficient to explain the differences in the aetiology and malignancy of HR-NBs relative to LR-NBs (4)(5)(6). Accordingly, the challenge remains to identify the genetic and molecular mechanisms that discriminate LR-NBs from HR-NBs.
Lately, much attention has been oriented towards the identification of the NB cell-oforigin, with the goal to determine whether LR-NBs and HR-NBs differ in the way they exploit the transcriptional program of this cell-of-origin. There is a consensus on the fact that NBs originate from alterations to cells in the sympatho-adrenal (SA) lineage (11)(12)(13)(14)(15)(16). Representing one of the multiple cell fate branches of neural crest cells, the SA lineage is composed of four main cell identities: Schwann cell precursors, bridge cells, chromaffin cells and sympathoblasts (12)(13)(14)(15)(16)(17). The specific identity of the NB cell-oforigin remains however debated (11)(12)(13)(14)(15)(16). Thus, more work is needed to elucidate whether LR-NBs and HR-NBs differ in the way they exploit the transcriptional program underlying their developmental SA origin, and even more importantly, if this SA program could be exploited to identify novel factors capable of impeding the metastatic dissemination and malignancy of HR-NBs.
Here, we addressed these questions by first comparing the genome-wide transcriptomic profiles of primary samples from LR-NB and HR-NB patients and human fetal adrenal gland, and thereby identified a transcriptional signature associated to NB formation that further distinguishes LR-NBs from HR-NBs. We then established that LR-NBs and HR-NBs can be discriminated by a signature corresponding to the core SA lineage, rather than by signatures specific of the distinct SA cell identities.
Remarkably, this core SA signature was composed in vast majority of genes associated with a favorable patient prognosis and with diminished disease progression. From these findings we selected Neurexophilin-1, which encodes the secreted glycoprotein NXPH1, as a candidate for functional studies and demonstrated that its activity, mediated by its receptors α-NRXN1/2, indeed represses NB malignancy by inhibiting the ability of NB cells to form metastases.

NB formation is correlated with a neural signature accentuated in LR-NBs.
We set out to define the transcriptional signatures associated to the aetiology of LRand HR-NBs and to the malignant behavior of HR-NBs. As such, we analyzed the transcriptome of primary tumor samples from 18  Table S1) (20). Samples of the human fetal adrenal gland (fAG) were used as a normal reference tissue (Fig. 1A, Supplementary Fig.S1A and Supplementary Table S1).
To identify the transcriptional signatures associated to NB aetiology, we first independently compared the LR-and HR-NB groups to the fAG samples and thereby  Table S2), in agreement with previous findings (11). By contrast, the DEGs found upregulated in these two comparisons were both characterized by a marked enrichment in GO terms related to the nervous system, or to the cell cycle and DNA repair (Fig. 1D, E and Supplementary Table S2). The GO terms related to the nervous system represented 68% of the top 50 GOs for the LR vs fAG comparison but only 18% for the HR vs fAG comparison (Fig. 1D, E and Supplementary Table S2), indicating that the formation of both LR-and HR-NBs is associated to a neural signature, and that this signature is more pronounced in LR-NBs. 6 We next searched for a transcriptional signature relevant to NB formation that furthermore distinguishes LR-NBs from HR-NBs. To this aim, we crossed the two independent comparisons LR vs fAG and HR vs fAG and thereby retrieved 3 Table S3). We then searched among these 3,096 common DEGs for the genes that were further differentially expressed between the LR-NB and HR-NB groups. This comparison identified a list of 503 LR vs HR DEGs that an unbiased hierarchical clustering analysis subdivided into 4 clusters of genes presenting distinct expression profiles ( Fig. 1F and Supplementary Table S3). The clusters c1 and c2 consisted of genes more strongly expressed in both LR-NBs and HR-NBs than in fAG samples. The best represented cluster (c1) contained 338 genes that were more strongly expressed in LR-NBs than in HR-NBs and were associated with a marked enrichment in GO terms related to the nervous system (  Table S3). Conversely, the 33 genes found in cluster c2 (such as CDCA7, C4orf46, THOC4) were more strongly expressed in HR-NBs than in LR-NBs and presented the features of oncogenes, as shown by their GO term enrichment profile and Kaplan-Meyer analyses (Fig. 1F, G, Supplementary Fig. S1E and Supplementary Tables S3 and S4). On the other hand, the clusters c3 and c4 consisted of genes expressed more weakly in both LR-NBs and HR-NBs than in fAG samples. The 32 genes found in cluster c3 presented the features of tumor suppressors, while the cluster c4 contained 100 genes that were more expressed in HR-NBs than in 7 LR-NBs and were mainly correlated with a bad prognosis (Fig. 1F, G, Supplementary   Fig. S1E and Supplementary Tables S3 and S4). Our analysis thus identified a complex transcriptional signature relevant to NB formation that furthermore distinguishes LR-NBs from HR-NBs. Remarkably, 67% of the 503 genes found in this signature, which formed the cluster c1, are both associated with a neural identity and a better patient outcome.
The core sympatho-adrenal signature is enriched in LR-NBs and associates with better patient prognosis.
We next assessed whether the transcriptional signature distinguishing LR-NBs and HR-NBs was related to their sympatho-adrenal (SA) origin. The human SA lineage can be subdivided into 4 main cell identities: Schwann cell precurors (SCPs), bridge cells, chromaffin cells and sympathoblasts ( Fig. 2A) (12)(13)(14)(15)(16). We searched for a possible enrichment of these SA cell identities within our list of LR vs HR DEGs using transcriptional signatures recently identified by single-cell RNA-seq during human adrenal gland development (16). We thereby observed that the 4 SA cell signatures all showed a strong enrichment for genes expressed at higher levels in LR-NBs than in HR-NBs, retrieving as many as 342 of our 503 DEGs (Fig. 2B and Supplementary Table   S5). By contrast, the signatures of the non-adrenal medulla cell types identified in the developing human adrenal gland did not show any enrichment, except the adrenal cortex signature which was enriched in genes expressed at higher levels in HR-NBs (Supplementary Fig. S2A and Supplementary Table S5). We noticed that the full SCP, bridge cell, chromaffin cell and sympathoblast signatures were largely overlapping (Fig.   2C and Supplementary Table S5). We retrieved a strong enrichment for genes expressed at higher levels in LR-NBs when we used a core SA signature consisting of the 4,518 8 genes shared by at least 3 of the 4 SA cell signatures (Fig. 2D and Supplementary Table   S5). By contrast, this enrichment was lost when using the signatures specifically consisting of the genes unique to SCPs, bridge cells or sympathoblasts, and the "unique" chromaffin cell signature showed a converse and weak enrichment towards genes expressed at higher levels in HR-NBs ( Supplementary Fig. S2B, C). Using the core SA signature we retrieved more than half of the 503 LR vs HR DEGs ( Fig. 2D and Supplementary Table S6). Even more strikingly, 92% of this list (242 genes out of 262 genes) corresponded to genes belonging to the cluster c1, thereby representing 72% of this whole cluster c1 (Fig. 2E and Supplementary Table S6). The transcriptional landscapes of LR-NB and HR-NBs can thus be discriminated based on a signature largely corresponding to the core SA lineage, and this developmental program is unexpectedly associated with a favorable prognosis.
Among these 242 genes common to the core SA signature and to the LR vs HR DEGs of the cluster c1 (thereafter named "SA-c1" genes) were found for instance the transcription factors (TFs) SOX6 and TFAP2 involved in the early stages of the NCC lineage (12,32), PRPH and the TFs GATA2 and FOXO3, which are associated with chromaffin cell differentiation (12), and various autism spectrum disorder genes (CNTNAP2, CTNND2, DLG2, NRXN2 and SHANK2) recently associated with NB aetiology (33) (Fig. 2E and Supplementary Table S6). To test the predictive potential of the SA-c1 gene module on NB prognosis, we subdivided the SEQC NB cohort into quartiles based on the combined expression of the 242 SA-c1 genes, whereby the quartiles Q1 and Q4 consisted of the patient samples presenting the highest and lowest numbers of SA-c1 genes expressed above average levels, respectively (see Materials and Methods). These quartiles showed a progressive decrease in event-free survival probability (8-years EFS probability of 0.88, 0.80, 0.49 and 0.27 for Q1, Q2, Q3 and 9 Q4, respectively, Fig. 2F). They also consisted of increasing proportions of INSS4 tumors, which are metastatic and aggressive (31), at the expense of INSS1-3 tumors usually associated with a good prognosis (Fig. 2G). Importantly, the positive correlation between SA-c1 gene expression and better prognosis was also observed when focusing on either INSS1-3 or INSS4 tumors ( Supplementary Fig. S3A, B). Remarkably, the expression levels of the SA-c1 module correlated to patient prognosis more strongly than the INSS classification. For instance, the group of INSS4 samples from the quartile Q1 presented a higher EFS probability than the groups of INSS1-3 samples from the quartiles Q3 and Q4 ( Supplementary Fig. S3A, B). We further sub-divided the SEQC cohort based on the progression of the disease, including within the category "progression" the tumors that did not respond to therapy plus those that were recurrent despite an initial therapeutic response. We thereby observed an inverse correlation between disease progression and the combined expression of the SA-c1 genes (Fig. 2H).
The module of SA-c1 genes related to the core sympatho-adrenal program is thus strikingly predictive of NB malignancy, patient prognosis and disease progression. The SA-c1 genes might thus play a crucial role in restraining NB malignancy.

The expression of NXPH1 and its receptors α-NRXN1/2 associates with favorable patient prognosis and identifies NB cells with a neural crest stem cell identity.
To functionally test the impact of SA-c1 genes on NB malignancy, we selected Neurexophilin-1 (NXPH1). NXPH1 was the SA-c1 gene with the second highest fold change enrichment in LR-NBs relative to HR-NBs (fc LR/HR=7.46; Fig. 2E and Supplementary Table S6) and its expression levels strongly correlated with favorable patient prognosis and with diminished disease progression in both the INSS1-3 and INSS4 NB sub-groups (Fig. 3A, B and Supplementary Fig. S3C-E). NXPH1 is mostly expressed in the nervous system and encodes a secreted glycoprotein that specifically binds to and modulates the activity of the α-Neurexin transmembrane receptors (α-NRXN1, 2 and 3; Fig. 3C), which are known to play key roles in synaptogenesis and neurotransmission (34)(35)(36)(37). Interestingly, NRXN2 also came out as a SA-c1 gene ( Fig.   2E and Supplementary Table S6) and both NRXN1 and NRXN2 expression levels associated with a favorable prognosis, whereas NRXN3 levels did not ( Fig. 3D and Supplementary Fig. S3F).
We first characterized the expression of NXPH1 and its α-NRXN1/2/3 receptors in a panel of 10 human NB cell lines with diverse genetic profiles and morphological properties ( Fig. 3E) (19,38). In basal culture conditions NXPH1 and α-NRXN1/2 were expressed weakly, while the expression of α-NRXN3 transcripts was below the threshold of detection. Growing NB cell lines in restrictive culture conditions (see Materials and Methods) for 5 weeks in vitro can lead to the formation of spheroids that are progressively enriched in cells presenting stem cell-like characteristics and a more specific NCC identity (39)(40)(41). Following this protocol, most of the NB cell lines (8 out of 10) showed the capacity to grow as spheroids (Fig. 3E). The formation of spheroids correlated to a stem cell-enrichment, as shown by increased levels of the transcripts encoding the stereotypic stem cell markers CXCR4, NANOG and KIT, the NCC stem cell-specific marker p75NTR and the MDR1 gene associated with chemotherapy resistance (Fig. 3E) (42)(43)(44). Remarkably, NXPH1 and α-NRXN1/2 levels increased in all the NB cell lines harbouring a sphere-forming capacity (Fig. 3E), thereby revealing a strong positive correlation between the expression of NXPH1 and α-NRXN1/2 and the acquisition of a NCC stem cell identity.
Taking advantage of a commercial fluorescence-conjugated antibody recognizing the extracellular region of human α-NRXN1, we identified a small subpopulation of α-11 NRXN1 + cells (less than 1.5%) within the 3 human NB cell lines that showed the highest sphere-forming capacity and in cells dissociated from 3 different NB patientderived xenografts (PDX; Fig. 3F). Purifying α-NRXN1 high cells from the SK-N-SH cell line or from the PDX NB-012 by FACS and growing them in restrictive culture conditions for a week revealed their increased spheroid-forming capacity compared to both α-NRXN1 low and α-NRXN1cells (Fig. 3G). In addition, α-NRXN1 high cells harboured an enhanced proliferative ability when seeded under conditions of extreme dilution, relative to purified α-NRXN1cells or unsorted SK-N-SH cells (Fig. 3H). To assess the importance of α-NRXN1 + cells in vivo, we compared the growth potential of SK-N-SH cells deprived of their α-NRXN1 + cell subpopulation with that of nondeprived cells using a CAM assay, in which cells were seeded on the richly vascularised chorio-allantoid membrane of 10 days-old chicken embryos (Fig. 3I). After 7 days of incubation in ovo, the number of cells quantified per tumor section was decreased by ~50% for the α-NRXN1 + -deprived cells relative to their control (Fig. 3J, K), thus revealing that α-NRXN1 + cells are required to support NB tumor growth in vivo. These findings established a correlation between NXPH1/α-NRXN expression and the tumorigenic potential of NB cells, arguing that NXPH1/α-NRXN signaling could control NB growth and/or aggressiveness.

NXPH1/α-NRXN signaling stimulates NB growth
To study whether NXPH1/α-NRXN signaling regulates NB growth, we generated stable clones of SK-N-SH cells that constitutively express sh-RNAs targeting either NXPH1 or α-NRXN1 ( Supplementary Fig. S4A-D). Both NXPH1 and α-NRXN1 knockdowns severely decreased cell viability over a 3 days-period in vitro and caused a complete growth arrest after one week ( Supplementary Fig. S4E, F). To circumvent this effect, we generated clones whose sh-RNA and eGFP production was only induced upon doxycycline treatment ( Supplementary Fig. S4G). Using eGFP + FACS-purified cells pre-treated with doxycyline for 96 hours, we observed that the growth arrest caused by both NXPH1 and NRXN1 knockdown was minimized in the inducible sh-NXPH1 and sh-αNRXN1 clones when grown in basal culture conditions (Supplementary Fig. S4H-O). These sh-NXPH1 and sh-αNRXN1 cells however retained a diminished sphere-forming capacity when grown in restrictive culture conditions ( Fig. 4A-C). We next tested the effect of these inducible sh-NXPH1 and sh-αNRXN1 clones on NB growth in vivo, by xenografting FACS-purified, doxycylinetreated cells into the flanks of immuno-compromised NOD/SCID mice that received doxycycline ad libitum (Fig. 4A). The sh-NXPH1 and sh-αNRXN1 cells produced fewer and smaller tumors than their control and caused a marked delay in the time required for their detection (Fig. 4D, E). Thus, inhibiting NXPH1/α-NRXN1 signaling strongly impairs NB tumor formation in vivo.
We next wondered whether stimulating NXPH1/α-NRXN signaling would be sufficient to increase the growth of NB tumors. To test this idea, we opted for a pharmacological gain-of-function approach and performed a CAM assay in which parental SK-N-SH cells were embedded in Matrigel in the presence or absence of recombinant human NXPH1 (rNXPH1, 10 µg/ml; Fig. 4F). After 7 days, we found that rNXPH1 treatment had increased the number of NB cells per tumor section by 62% Interestingly, the enhanced proliferation rate caused by rNPHX1 was associated with a 2-fold increase in the proportion of NB cells expressing the neural crest stem cell marker p75/NTR ( Supplementary Fig. S5G, H). Together, these data revealed that NXPH1/α-NRXN1 signaling is necessary and sufficient for NB tumor growth in vivo ( Fig. 4K).

NXPH1/α-NRXN signaling inhibits the metastatic potential of NBs
Finally, we set out to define whether NXPH1/α-NRXN signaling regulates the metastatic potential of NB cells in vivo. To this end, we generated clones of SK-N-SH cells carrying a constitutively-expressed luciferase cassette in addition to the doxycycline-inducible sh-RNAs and eGFP. FACS-purified, doxycycline-treated cells from the sh-NXPH1, sh-αNRXN1 and sh-Ctrl clones were injected into the left cardiac ventricle of immunodeficient NOD-SCID gamma (NSG) mice and their dissemination and metastatic growth were monitored non-invasively by bioluminescence (BLI) over 9 weeks (Fig. 5A). A few minutes after injection, cells from the sh-Ctrl, sh-NXPH1 and sh-αNRXN1 clones were seen disseminating within the body of the mice (Fig. 5B).
Their BLI then regressed to background levels during the first 3-4 weeks, indicating that the vast majority of the cells disappeared during that period (Fig. 5B). After 9 weeks, 86% and 57% of the mice injected with sh-NXPH1 and sh-αNRXN1 cells had developed detectable metastatic tumor masses, whereas only 29% of the mice injected with the sh-Ctrl cells had (Fig. 5B, C). Moreover, the metastatic tumors developed by the sh-NXPH1 and sh-αNRXN1 cells were detected earlier and produced a photon count much higher than the sh-Ctrl cells (Fig. 5B, D). Metastases were particularly evident in the liver and bone marrow, in agreement with previous reports on the 14 organotropism of metastatic NBs (5,45). Thus, reducing the expression of NXPH1 or that of its α-NRXN1 receptor stimulated the metastatic potential of NB cells in vivo.
These latter findings therefore revealed that NXPH1/α-NRXN1 signaling represses the dissemination and metastatic potential of NBs (Fig. 5E).

Discussion
Our study provides fundamental insights into the transcriptional signatures associated with both the aetiology and malignancy of NB. We notably reveal that a signature corresponding to the core SA lineage strongly distinguishes LR-NBs from HR-NBs, this signature being largely composed of genes associated with favorable patient outcome and with reduced disease progression. Encoded by one of these SA-c1 genes, the secreted protein NXPH1 and its receptor α-NRXN1 stimulate the growth of NB cells but markedly restrict their metastatic potential.
The pathogenesis of NB remains puzzling for several reasons. First, some NB tumors can regress spontaneously (46,47). Second, so far there is little evidence showing that LR-NBs ever evolve into HR-NBs once diagnosed (48). This led to the idea that LR-NBs and HR-NBs represent distinct types of tumors and that they likely originate at different stages of development within the neural crest lineage. Various studies recently employed single-cell transcriptomics to tackle this question and to identify the NB cellof-origin, which yielded discrepant conclusions (12)(13)(14)(15)(16)(49)(50)(51)(52). Whereas one study argued that adrenal NBs derive from fetal chromaffin cells and that LR-NBs and HR-NBs can be discriminated based on the degree of chromaffin cell differentation (12), others pinpointed instead fetal sympathoblasts as the NB cell-of-origin (11,(14)(15)(16).
Another study found that the signatures of sympathoblasts and norepinephrinechromaffin cells were highly enriched in samples from NB patients, but suggested that the signature of Schwann cell precursors is the one distinguishing LR-NBs from HR-NBs (13). As we herein focused directly on identifying the transcriptional signature associated to NB aetiology that further distinguishes LR-NBs and HR-NBs, we did not attempt to settle this current debate. We specifically addressed the component of the SA lineage that distinguishes LR-NBs and HR-NBs, using transcriptional signatures for human SCPs, bridge cells, chromaffin cells and sympathoblasts recently identified by single-cell RNA-seq (16). In contrast to previous reports (12)(13)(14)(15)(16), we tested these signatures on a restricted list of LR-NB vs HR-NBs DEGs and observed that these 4 signatures all showed a strong enrichment for genes expressed at higher levels in LR-NBs than in HR-NBs. More importantly, having realised that the SCP, bridge cell, chromaffin cell and sympathoblast signatures that we used were largely overlapping, we then demonstrated that their enrichment for genes expressed at higher levels in LR-NBs is due to the signature shared by these sympatho-adrenal cell identities, rather than to the transcriptional singularities of any of these four cell types. Notably, this core SA signature retrieved more than half of the 503 LR vs HR DEGs, highlighting the importance of this developmental signature in discriminating the transcriptional landscapes of LR-NB and HR-NBs. Remarkably, this core SA signature retrieved 72% of the LR vs HR DEGs forming the cluster c1 (i.e. the genes expressed at higher levels in LR-NBs than in HR-NBs than in fAG samples and associated with a favorable patient outcome). This preferential correlation between the transcriptional identity of the SA lineage and the LR-NB phenotype suggests that this developmental program actually opposes NB malignancy, which was unforeseen. Put together, our findings shed a new light on the complex contribution of the neural crest-derived developmental program to NB pathogenesis. The SA program indeed appears to facilitate NB formation but to concomitantly restrict its malignant potential. The dual impact of NXPH1/α-NRXN signalling on NB cell behaviour presented here is a first experimental illustration of this unexpected model. NXPH1 came out in our genome-wide transcriptomic analysis as the SA-c1 gene with the second highest enrichment in LR-NBs relative to HR-NBs. The intensity of NXPH1 expression stratifies NB patients relative to their tumor stage and to the probability of disease progression, extending findings reported in a genome-wide associated study that specifically compared INSS-3 and INSS-4 NB tumors (53). Interestingly, NXPH1 can also be used as a DNA methylation biomarker associated with a good prognosis for NB patients (54), and represents one of the genes differentially expressed between LR-NBs and HR-NBs with the most restricted expression pattern in postnatal human tissues (16).
Therefore, monitoring NXPH1 levels might be clinically relevant for the therapeutic management of patients suffering from NB.
The secreted glycoprotein NXPH1 specifically binds and modulates the activity of the transmembrane receptors α-NRXNs (34)(35)(36)(37). The fact that NRXN2 also came out as a SA-c1 gene and that the expression levels of both NRXN1 and NRXN2 correlate to a better prognosis is noteworthy. Although the data available did not allow us to determine if these correlations were directly related to α-NRXN isoforms, these findings motivated our choice to focus on NXPH1/α-NRXN signaling for functional studies. The results from our functional assays demonstrated that both NXPH1 and α-NRXN1 regulate various aspects of NB cell behaviour. First, they stimulate NB tumor growth and we provide evidence that NXPH1 fuels the proliferation of NB cells. Second, NXPH1 and α-NRXN1 knockdowns both unleash the ability of NB cells to colonise tissues like the liver and bone marrow, two of the main metastatic organs in HR-NB patients (5,45). The molecular mechanisms causing these two opposite effects remain unclear and their study is made harder by the fact that the signaling cascade triggered downstream of α-NRXNs in response to NXPH1 binding is as yet unknown.
Nevertheless, this unexpected anti-metastatic activity of NXPH1/α-NRXN signaling likely explains why higher NXPH1 and αNRXN1/2 expression levels are correlated with a better patient prognosis. Targeting and modulating NXPH1/α-NRXN signaling might thus represent an appealing strategy to treat metastatic, therapy-resistant HR-NBs.

18
In summary, our work reinforces the importance of the neural crest-derived SA developmental program in NB pathogenesis and especially highlights the necessity to examine how these genes regulate both the growth and the metastatic potential of NBs.
This will facilitate the identification of novel actors, such as NXPH1/α-NRXN1 signaling, that might serve as therapeutic targets to treat children affected by metastatic, therapy-resistant HR-NBs.

Animal Studies
Fertilized white Leghorn chicken eggs were provided by Granja Gibert, rambla Regueral, S/N, 43850 Cambrils, Spain. Eggs were incubated in a humidified atmosphere at 38°C in a Javier Masalles 240N incubator, manipulated at embryonic day 10 (E10) and sacrificed at E17. Sex was not identified. HSJD-NB-011 and HSJD-NB-007 have been previously used for publication (18,19). The HEK-293 cell line was grown in DMEM-Glutamax (Lifetechnologies cat#61965-026) supplemented 10% FBS and 1% penicillin/streptomycin cocktail. Unless indicated, cells were grown in standard culture conditions. The presence of mycoplasma contamination was routinely checked by immunofluorescence. Cell lines have not been re-authenticated for the present paper.

Genome-wide transcriptomic analysis
The HSJD-NB dataset used in this study has been published previously (20) and is available at GEO repository (GSE54720). The full human SCP, bridge cell, chromaffin cell and sympathoblast signatures were obtained from a published study (16), and the complementary signatures were obtained from a Venn diagram analysis performed with a freely available software developed by the group of Dr Y. Van de Peer (VIB, Brussels, Belgium, http://bioinformatics.psb.ugent.be/webtools/Venn/). Genome-wide transcriptomic analyses were performed using the web-based Phantasus software (https://artyomovlab.wustl.edu/phantasus; version 1.5.1) (21). Normalized log expression values were obtained using the Quantile Normalize Adjustment tool. The Maximun Median Probe method was then selected to collapse the different probes of a single gene. Finally, genes were filtered based on raw expression levels and the 12,000 most-expressed genes from each array were selected for subsequent analysis. Sample dispersion was then assessed using a principal component analysis leading to the removal of one outlier (#LR-08) from further analysis. Differentially-expressed genes (DEGs) between groups were identified using the Limma R package. The genes were considered as differentially expressed when the false discovery rate (FDR)-adjusted Pvalue (with Benjamini-Hochberg procedure) was <0.05. Gene set enrichment analyses (GSEA) were performed using the FGSEA tool from Phantasus.

Gene ontology analysis
Gene ontology (GO) term enrichment analyses (biological process) were performed with the PANTHER classification system (http://pantherdb.org) (22,23). The 50 mostenriched GO terms were used to assess GO term distribution and were annotated with the Ancestor Chart tool of QuickGO (https://www.ebi.ac.uk/QuickGO).

R2 genomic analysis and visualization platform
Publically available NB patient microarray SEQC 498 (GSE62564) was obtained from the R2 genome analysis and visualization platform (https://hgserver1.amc.nl/cgibin/r2/main.cgi) (24). The R2-web based application was used to generate Kaplan-Meier event-free survival curves. Data were grouped based on the expression levels of either NXPH1, NRXN1-3 or the full SA-c1 gene module. In the latter case, the SEQC cohort was subdivided into quartiles based on the combined expression of the 242 SA-c1 genes, considering as "high" or "low" the expression of each SA-c1 gene with respect to its own average expression level. The quartiles Q1-Q4 thereby consisted of 123-126 samples showing numbers of SA-c1 genes expressed at high levels as follows: 487≤Q1≤360; 359≤Q2≤250; 249≤Q3≤135 and 134≤Q4≤4. Statistical significance was automatically assessed by the R2-server using a log-rank test. The R2-web based application was also used to compare the distribution of patients among different NB prognosis groups according to expression levels of the SA-c1 gene module or NXPH1.
For the analysis of disease progression, the category "progression" included tumors that did not respond to therapy plus those that were recurrent despite an initial therapeutic response. Statistical significance was automatically assessed by R2-server using the Chi-square + Fisher's test. following manufacturer's instruction and lentiviruses were harvested 48h posttransfection. Transduction of target cells was usually performed readily after lentivirus recovery. 2μg/ml of puromycin (SIGMA #p8833) and 600μg/ml of neomycin/G418 (SIGMA-ALDRICH #A1720) were used for selection of transduced cells.

Fluorescence-activated cell sorting (FACS)
All experiments were conducted using a BD FACSAria Fusion at the Cytometry Core

MTT assay
Approximately 8,000 cells were seeded in quadruplicates. At desired timepoints, 10μl of 5mg/ml MTT solution (SIGMA #M2128) was added to each well and incubated at 37°C for 4 hours. Formazan crystals were dissolved with 100μl DMSO (SIGMA #D2650) and the optical density at 570nm (OD570) was measured. This assay was performed on n=2-4 biological replicates.

S-phase index and BrdU-retention assay
A solution containing 10 μM of 5-bromo-2'-deoxyuridine (BrdU; SIGMA #858811) was added to the growth medium of SK-N-SH cell monolayers and cells were incubated at 37ºC for 2 hours to estimate the S-phase index, or for 24 hours prior to the CAM assay to assess BrdU retention, using n=4-9 biological replicates. Cells were then fixed, permebilized and incubated with 5U/ml DNase I type II (SIGMA #D4527) for 15 min at RT. BrdU incorporation was then assessed by immunofluorescence as described below.
The sphere-forming capacity of the different human NB cell lines was performed on n=2 biological replicates.
To assess mRNA expression in stem cell-enrichment conditions, 3•10 6 cells were grown on uncoated 60mm plates in 2ml SFM. Spheres were allowed to grow for 5 weeks and were photographed using an inverted bright-field microscope (Leica DMIRBE). Viable spheres were recovered by centrifugation at 100g for 5min at 4ºC) and processed for RNA extraction as described below.
To assess the sphere-forming capacity of NB cells, 2,000 dissociated cells from each population of interest were sorted in quadruplicates and seeded directly onto 12-well plates coated with 0.5% agarose (Condisa Laboratory #8014). Cells were allowed to grow for 1 week and the presence of spheroids was monitored every two days. To detect the presence of viable spheroids after 1 week, plates were incubated with 0.5mg/ml MTT (SIGMA #M2128) for 3 hours at 37°C. Viable spheres (black coloured) were manually counted using an inverted bright-field microscope (Leica DMIRBE). The sphere-forming capacity was assessed for n=1-4 biological replicates per experimental condition.

2D Extreme limiting dilution assay (ELDA)
Cells subpopulations were FACS-sorted, seeded in triplicates into 96-well plates and grown until one experimental condition fully covered the well surface. Half of the growth medium was renewed every week. Cells were then fixed with 4% PFA (SIGMA #16005) for 15min at RT, washed with PBS and stained with 0.5% crystal violet (SIGMA #V5365) for 20 min at RT. The dye excess was washed out and plates were air-dried for at least 2 hours at RT. Crystal violet staining was dissolved by adding 200μl methanol to each well, followed by 20min at RT. Finally, the OD570 was measured. This assay was performed on n=3 biological replicates per experimental condition. To test the growth of NB cells depleted from their α-NRXN1 + subpopulation, parental SK-N-SH depleted or not from their α-NRXN1 + subpopulation were FACS-purified and resuspended in grafting media. In both experiments, a concentration of 0.5·10 6 cells/10μl was used and grafted onto the CAM, using n=4-13 biological replicates. At E17 NB tumour grafts were removed from the CAM and tumours were photographed, weighed and measured using a digital calliper (VWR). The final tumour volume (V) was calculated using the following formula: V= 4/3·π·length·depth·width. Tumours were then fixed in 4% PFA for 1h30min at 4°C, washed in PBS1x and sequentially cryo-protected with PBS solutions containing 15% and 30% sucrose. Small cubical tumour pieces (5mm·5mm·5 mm approx.) were then embedded in Tissue-tek (Sakura

Immunofluorescence
Finetek #4583) and frozen on dry ice. Tumour pieces were sectioned on a cryostat (Leica) at 16μm, collected serially on home-made TESPA pre-coated slides and slides stored at -20°C. Immunofluorescence was performed as described above.

Statistics
No statistical method was used to predetermine sample size. The experiments were not randomized. The investigators were not blinded to allocation during experiments.    Availability of data and materials: All data generated or analysed during this study are included in this article and its supplementary information files.

Competing interests:
The authors have no competing financial interests to declare.
Correspondence and requests for materials should be addressed to G.L.D.       significance was assessed using a non-parametric Mann-Whitney test (C) Figure S4: Cellular behaviour of constitutive and inducible sh-NXPH1 and sh-αNRXN1 clones in vitro. (A) Schematic representation of the pLKO.1-Puro construct used to generate stable clones of SK-N-SH cells constitutively expressing sh-RNAs specifically targeting human NXPH1 or α-NRXN1. (B) Experimental design applied in C-F. (C, D) Relative transcripts levels of (C) NXPH1 and (D) α-NRXN1 quantified by RT-qPCR in the sh-Ctrl clone and the two sh-NXPH1 and two sh-αNRXN1 clones generated. The values represent the mean relative mRNA levels ± s.d calculated from 3-4 independent experiments. (E) MTT assay performed over a time-course of 72 hours to assess the cell viability of the sh-Ctrl (black) clone and the two sh-NXPH1 (blue) and two sh-αNRXN1 (red) clones generated. The values represent the mean OD570-750 ± s.d obtained from 4 independent experiments. (F) Bright field images of the sh-Ctrl clone and the two sh-NX-PH1 and two sh-αNRXN1 clones, obtained 96 hours after puromycin selection. (G) Schematic representation of the pSLIK-Neo-TTGmiRc2 construct used to generate stable clones of SK-N-SH cells that produce EGFP and sh-RNAs specifically targeting human NXPH1 or α-NRXN1 in an inducible manner upon doxycyline treatment. (H) Experimental design applied in (I-O). (I) Representative images showing how EGFP production is triggered by doxycycline treatment (+DOX, after 2 days) and can be repressed again if doxycycline is removed (-DOX, 3 days later). (J) MTT assay performed over a 96 hours time-course to assess the cell viability of the inducible sh-Ctrl (black), sh-NXPH1 (blue) and sh-αNRXN1 (red) clones. The values represent the mean OD570-750 ± s.d obtained from 2-4 independent experiments. (K) Representative images of the inducible sh-Ctrl (black), sh-NXPH1 (blue) and sh-αNRXN1 (red) clones, their mean number (L) of EGFP+ cells ± s.d quantified per field and their mean proportions of (M) Ki67+;EGFP+/total EGFP+ cells ± s.d, (N) pH3+;EGFP+/total EGFP+ cells ± s.d and (O) pyknotic nuclei EGFP+/total EGFP+ cells ± s.d, quantified 7 days after seeding. Significance was assessed with the non-parametric Mann-Whitney test (C, D) or a 2-way ANOVA + Dunnett's test (E, J, L-O). *p<0.05, **p<0.01, ***p<0.001, ns: p>0.05. Scale bar: 10µm (K), 100µm (I).