The study included a total of 172 patients with histologically confirmed NB. A written consent allowing the collection of samples and the use of clinical and nongenetic data for clinicalresearch was signed by the patient’s guardians. The study was approved by the Regional Ethical Committee (ANTECER_Neuroblastoma 16/09/2019).

All NB patients were diagnosed between 2011 and 2018 in 29 institutions belonging to the Associazione Italiana di Emato-Oncologia Pediatrica. NB patients were assigned to the different risk groups according to the International Neuroblastoma Risk Group (INRG) pre-treatment risk schema and addressed to the different treatment protocols (2). Fasting plasma samples were obtained from peripheral venous blood collected in 3 mL EDTA K3-containing tubes, centrifuged at 4000 g for 5 min at 4°C and stored at -80°C until analyzed. Leftover plasma samples of outpatients after routine clinical analyses were used as healthy controls.

Samples were divided in two sets: the first set (n=99) was used for biomarker discovery purposes while the second set (n=122) was used for biomarker validation. The first set of samples was obtained from 27 control subjects (CTR) and 50 NB patients among which 28 patients with localized NB (INRG L1, L2) or MS and 22 with metastatic NB (INRG M) which have been sampled at two timepoints: at diagnosis and after induction chemotherapy. The second set consisted of 122 samples from 69 patients with localized NB and 53 with metastatic NB (INRG M). Table 1 summarizes the clinical demographics of the two sets.

In the first set, sex (male vs female) and age group (<18 months vs >18 months) were tested for confounding between CTR subjects and NB patients using Fischer’s exact test; P value less than 0.05 was considered significant ( Supplementary Table 1 ).

Plasma samples used for untargeted metabolomic analyses were prepared as follows: a 50 µL aliquot of plasma was extracted adding 150 µL of cold (-20°C) methanol containing an internal standard (IS) mixture ( Supplementary Materials, 1.1 Chemicals and consumable) and centrifuged for 10 minutes at 14,000 x g at 4°C. The supernatant was then collected and stored at -80°C until analyzed. A 50 µL aliquot of supernatant was then added to an equal volume of IS kit (MSK-QC-KIT) purchased by Cambridge Isotope Laboratories, Inc (Tewksbury, MA, USA) ( Supplementary Materials, 1.1 Chemicals and consumable), vortex mixed and 5 µl were injected in the UHPLC system. Quality control (QC) samples were prepared by pooling together supernatants obtained from all the samples; external quality control (EQC) samples were prepared by pooling together plasma obtained from healthy adult volunteers extracted as described above.

LC-HRMS analysis was carried out using a Vanquish Horizon UHPLC system coupled to a Q-Exactive Plus Hybrid Quadrupole-Orbitrap Mass Spectrometer (Thermo Fisher Scientific, Milan, Italy) as previously described (11). The liquid chromatographic separation was carried out using two different chromatographic conditions: reversed phase and hydrophilic interaction liquid chromatography (HILIC) ( Supplementary Materials, 1.2 Untargeted metabolomic analysis). Ionization was obtained using a heated electrospray source probe both in positive and negative mode with spray voltages at 3.9 kV and 3.7 kV, respectively. The capillary temperature was set at 300°C. Nitrogen sheath and auxiliary gas were set at a flow rate of 35 and 5 arbitrary units respectively. Extracted samples were splitted into four injection batches, and processed combining the four different conditions mentioned above (two different chromatographic conditions and two polarities). For each batch, analyses were performed using two different acquisition modes with an m/z range 70–1200. The first acquisition mode, used for compound profiling purposes, was a full MS scan mode with the following parameters: resolution of 70000, auto gain control target < 1 × 10⁵ and maximum injection time 100 ms. The other acquisition mode, used for identification of unknown compounds, was a full MS scan acquisition followed by data-dependent fragmentation (MS2) scan (DDA) with a resolution of 17500, an auto gain control target of 1 × 10⁵, a maximum injection time of 65 ms, a loop counts of top 5 peaks and an isolation window of m/z 1.2. All MS2 spectra of the compounds were acquired at 3 collision energies: 20, 40 and 80 eV. The MS1 mass range was divided into 10 mass ranges with a width of 50 m/z each on which the data dependent analysis was carried out. DDA was performed with a priority fragmentation for the m/z of our Accurate Mass Retention Time (AMRT) library (11).

The run order used to perform the metabolomic analyses started with 5 procedural blank samples composed by the solvents used during all steps of sample preparation, for contamination monitoring and system conditioning purposes; 18 QCs analyzed in DDA mode, for identification purposes, followed by the runs of the 99 samples that were randomized in the analytical sequence to avoid bias due to instrument drift; finally, 16 QCs and 9 EQCs added every 6 and 12 runs respectively. The experimental design and the analytical workflow are summarized in the Supplemental figures ( Supplementary Figure 1 ).

MS1 full-scan .RAW files have been converted into .ABF files using ABF converter (14) and processed with MS-DIAL ver.4.24 software (15) for deconvolution, peak picking, alignment and compound identification. Compounds, initially, were identified by matching retention times, accurate precursor masses, and MS/MS spectra with a previously published AMRT database (11) to provide a level 1 identification of metabolites (16). For MS2 spectra and precursor mass matching, the freely available library MassBank of North America (17) and mzCloud (18) were also used. The four tables obtained by MS-DIAL, generated combining the two different chromatographic methods and the two different MS polarities were exported as a .txt file.

The dataset, composed of four different sample sets and QCs was imported and statistically analyzed in Perseus (19). The QC variation coefficient was calculated and then used to exclude features variations > 0.4. Subsequently, the expression data was log2 transformed and the QC group removed. Each feature was, therefore, filtered by row with a valid value of at least 70% for each group. The missing data were replaced by the QRILC method, which performs the imputation of left-censored missing data, using random draws from a truncated distribution with parameters estimated by quantile regression. Quantile normalization was instead used to normalize the intensity values which were further filtered with the MS/MS assigned category, considered necessary to continue the analysis. Indeed, for the annotations of in silico compounds, starting from ionic features with associated MS2 events, MS-FINDER ver.3.26 was used (20). The MS-FINDER annotation matrix was merged with the quantitative matrix by employing the matching row by name Perseus option and the Alignment ID was used as the unique classifier. A MS-FINDER Structure Rank score > 5 and a mzCloud identification score > 80 were used to select features that were worthy of statistical investigation.

A differential analysis was then performed using the Limma Test package, the parameters are the same for all chromatographs and use the Voom option for data editing with a span of 0.5, normalization is quantile, with an acceptability threshold for up-regulation and down-regulation of 2 with respect to a Log2 Fold Change and an adjusted p-value (FDR) of less than 0.05.

Compounds with a Human Metabolome Database (HMDB) accession were matched with the HMDB database (“All Metabolites’ ‘ file dated 2021-10-23) (21) to add categorical information to the expression matrix like associated proteins. Statistically significant metabolites (and associated proteins) are used in enrichment analyses such as in ClueGO (22), a Cytoscape App (23), and MetaboAnalyst (24) to visualize their associated pathways.

Targeted quantitative analysis of 3-O-MD was performed on the second set of patients by LC-MS/MS on a TSQ Quantiva mass spectrometer coupled to an Ultimate 3000 UHPLC (Thermo Fisher, Scientific, Milan, Italy). Briefly, 50 µL of plasma were subjected to protein precipitation with 150 µL cold (-20 ∘C) methanol containing 3-o-methyldopa-d3 as internal standard. A 5 µl aliquot of the supernatant was then injected in the LC-MS/MS system using a gradient separation chromatography. Method performance was validated following EMA guidelines for bioanalytical method validation (25) ( Supplementary Materials, 1.3 Targeted analysis of 3-O-methyldopa by LC-MS/MS and 1.4 Mass spectrometry methods validation and quality assurance).

Overall survival (OS) curves were plotted by the Kaplan-Meier method and compared with the log-rank test using GraphPad Prism ver. 8 (26). For improving reliability, survival analysis was performed on alive patients with at least 3 years of follow-up. Patients lost at the follow-up or with missing survival information were excluded from the analysis. Cutoff values distinguishing between low or high 3-O-MD concentrations levels were identified by Elbow method (27), which is an empirical method that, given a numerical variable, selects one or more cut-off values whose slope change is evident. Significance of the expression differences of 3-O-MD across tumor stages was assessed by one-way ANOVA test. Post-hoc analysis was performed by Tukey’s multiple comparisons test.

Raw data were deposited in MetaboLights Data Repository (28) and assigned the identification code MTBLS4294. A total of 99 metabolic profiles were obtained from plasma samples of the first cohort from untargeted metabolomic analyses and metabolite identification as described in previous paragraphs. Filtering of data resulted in a significant reduction of the number of features to be queried by statistical analysis ( Figure 1 ). Features were then selected for CV, retention times and MS2 annotation. Quantitative differences in metabolites expression were detected across different groups of patients and controls by using the Limma test. The statistically significant metabolites were annotated with related proteins and metabolic pathways in order to understand their involvement in biological processes ( Supplementary Figure 2 ). The expression of the metabolism of Tyrosine, dopamine and catecholamines should be highlighted. The main reactions of the metabolic process take place through the activity of, among others, the enzymes Catechol-O-Methyltransferase, Aldehyde Dehydrogenase and Monoamine Oxidase ( Supplementary Figure 3 and Supplementary File 1 ). Another path that is important to mention is the metabolism of Arginine and Proline. Recent studies on cancer metabolism have highlighted the role shown by proline metabolism, in particular its critical role in cancer reprogramming and its clinical relevance (29).

In order to estimate possible confounding effects between CTR subjects and NB patients, we studied the association between sex (M vs F) or age group (<18months vs >18 months) and subject group (CTR vs NB). Fischer’s exact test p value did not show a significant association between sex or age group and the subject group (p>0.05), thereby excluding the possibility that sex or age might be confounding factors in our study.

As shown in Figure 2A , the differential expression analysis between samples from NB patients at the onset and controls revealed the presence of 90 significantly modulated metabolites ( Supplementary Table 2 ) belonging to different metabolic pathways, among which L-cystathionine and 2-Hydroxy-3-methylbutyric acid were evident. The comparison between localized (INRG stage L1 and L2) and metastatic (INRG stage M) patients at diagnosis highlighted 26 significantly modulated metabolites ( Supplementary Table 3 and Figure 2B ). In particular, Homovanillic acid sulfate and 3-O-MD, metabolites involved in the L-DOPA metabolism, were found to be overexpressed in metastatic patients. On the contrary, Metanephrine resulted significantly down expressed in metastatic patients. Noticeably, in metastatic patients, L-cystathionine and 2-Hydroxy-3-methylbutyric acid were overexpressed. As shown in Figure 2C , the differential expression analysis using samples from metastatic patients at diagnosis and after induction chemotherapy highlighted the significant modulation of 30 metabolites ( Supplementary Table 4 ) and pointed out the expression alteration of three metabolites of L-DOPA metabolism pathway: 3-O-MD, Homovanillic acid sulphate and Vanillylmandelic Acid, which resulted upregulated in metastatic condition. Moreover, L-Cystathionine and 3-Hydroxyisobutyric acid were also overexpressed at diagnosis, as well as, a number of metabolites of the methionine metabolism (adenosine, L-cystathionine, methionine sulfoxide and spermidine) and polyamines (Spermine and spermidine).

Since 3-O-methyldopa was found significantly modulated in several comparisons and is a metabolite of the catecholamine pathway mentioned in literature (30), we decided to focus on it in subsequent analyzes.

SupplementalThe concentrations of 3-O-MD by LC-MS/MS ranged between 23.6 and 6272.8 ng/ml. In order to validate the association between 3-O-MD expression and metastatic NB that was observed in our untargeted analysis, the distribution of 3-O-MD expression was reported across the subsets of patients defined by the INRG stage ( Figure 3 ). One-way ANOVA was used to estimate the significance of the 3-O-MD expression differences across stages. Post-hoc analysis using Tukey’s correction method was used to identify groups of patients with significantly different 3-O-MD expression. ANOVA test showed significant 3-O-MD expression differences across stages (ANOVA p value<0.05). Post-hoc analysis showed a statistically significant difference between stage M and stage L2 tumor subsets (p<0.05), thereby confirming our findings on the first cohort of NB patients. Visual inspection of the boxplot showed a higher 3-O-MD expression in the subset of stage M tumors with respect to other tumor stages ( Figure 3 ).

Since NB patients with stage M tumors are characterized by a disseminated disease and poor prognosis (2), we hypothesized a potential prognostic value of 3-O-MD at the onset for NB.

To estimate whether NB patients might be divided into groups on the basis of 3-O-MD expression, we split patients into two groups using the Elbow method (27), which was able to return the cut-off point from an ordered numerical variable. The application of the Elbow method to 3-O-MD expression identified 454.9 as a candidate cut-off value ( Figure 4A ). This cut-off value divided the cohort into two groups of 91 and 31 patients with high or low expressions of 3-O-MD, respectively. Unpaired t test assessed the significance of 3-O-MD expression difference between the newly defined groups of patients. The analysis highlighted a statistically significant difference of 3-O-MD expression between the two groups of patients (p<0.0001, Figure 4B ) and confirmed that NB patients may be divided into groups based on 3-O-MD expression.

The distribution of the two populations of patients with low or high 3-O-MD expression in the subsets defined by age at diagnosis (<18 months vs. > 18 months), INRG stage, and MYCN status (amplified vs. not amplified) was then evaluated. The number of patients with high 3-O-MD expression was greater than zero in all subsets except for patients with stage L1 tumor, but it was higher in the subsets of patients older than 18 months, INRG stage M, and amplified MYCN tumors ( Figure 5 ). These results indicate that high 3-O-MD expression is associated with unfavorable clinical characteristics. Furthermore, a non-trivial number of patients older than 18 months, INRG stage M, and amplified MYCN tumors had a low 3-O-MD expression suggesting that 3-O-MD may additionally stratify these subsets of patients. In order to test this hypothesis, we assessed the stratification of selected subsets of patients defined by known prognostic markers.

In order to improve analysis reliability, survival analysis was performed on alive patients with at least 3 years of follow-up excluding from the analysis patients lost at the follow-up or with missing survival information. Data from remaining 75 patients were then used for subsequent analyses. High or low 3-O-MD expression was able to stratify the entire cohort of 75 patients into groups with significantly different overall survival (log-rank p < 0.05, Figure 6A ). High expression of 3-O-MD metabolite was associated with worse prognosis indicating that 3-O-MD can be considered as a new unfavorable prognostic factor for NB.

Next, the prognostic value of 3-O-MD was assessed in additional clinically relevant subgroups of patients defined by established prognostic markers. Specifically, 3-O-MD significantly stratified patients with stage M tumors, in which, high expression of 3-O-MD metabolite was associated with worse prognosis (log-rank p < 0.05, Figure 6B ).

We then assessed the clinical utility of 3-O-MD in the subset of high-risk patients older than 18 months with stage M tumor that are difficult to be stratified with actual risk factors (31) and we found that 3-O-MD was able to significantly stratify them (log-rank p ≤ 0.05, Figure 6C ).

Our findings highlighted the ability of 3-O-MD to stratify clinically relevant subsets of patients, thus supporting the clinical utility of 3-O-MD metabolite for NB.