Human Melanoma-Cell Metabolic Profiling: Identification of Novel Biomarkers Indicating Metastasis

Melanoma is the most aggressive type of skin cancer, leading to metabolic rewiring and enhancement of metastatic transformation. Efforts to improve its early and accurate diagnosis are largely based on preclinical models and especially cell lines. Hence, we herein present a combinational Nuclear Magnetic Resonance (NMR)- and Ultra High Performance Liquid Chromatography-High-Resolution Tandem Mass Spectrometry (UHPLC-HRMS/MS)-mediated untargeted metabolomic profiling of melanoma cells, to landscape metabolic alterations likely controlling metastasis. The cell lines WM115 and WM2664, which belong to the same patient, were examined, with WM115 being derived from a primary, pre-metastatic, tumor and WM2664 clonally expanded from lymph-node metastases. Metabolite samples were analyzed using NMR and UHPLC-HRMS. Multivariate statistical analysis of high resolution NMR and MS (positive and negative ionization) results was performed by Principal Component Analysis (PCA), Partial Least Squares-Discriminant Analysis (PLS-DA) and Orthogonal Partial Least Squares-Discriminant Analysis (OPLS-DA), while metastasis-related biomarkers were determined on the basis of VIP lists, S-plots and Student’s t-tests. Receiver Operating Characteristic (ROC) curves of NMR and MS data revealed significantly differentiated metabolite profiles for each cell line, with WM115 being mainly characterized by upregulated levels of phosphocholine, choline, guanosine and inosine. Interestingly, WM2664 showed notably increased contents of hypoxanthine, myo-inositol, glutamic acid, organic acids, purines, pyrimidines, AMP, ADP, ATP and UDP(s), thus indicating the critical roles of purine, pyrimidine and amino acid metabolism during human melanoma metastasis.


Introduction
Melanoma is a cancer that originates from melanocytes, the specialized melanin-producing cells that mainly reside in the skin [1]. Incidence of skin cancer, the 3rd most common human of 50 metabolites were identified across all samples, with 39 of them being recognized to significantly differ between the two melanoma cell lines in the statistical tests.

Chemometric Analysis of NMR Data
Multivariate statistical analysis was performed, by employing SIMCA P 14.1. In the first step, Principal Components Analysis (PCA) was applied, in order to visualize the obtained data, with the aid of recognizing possible trends and locating possible outliers. Both Unit Variance (UV) and Pareto scaling exhibited clear separation between the two cell lines, WM115 and WM2664, without the presence of outlying values (UV scaling: R 2 X(cum) = 0.712 and Q 2 (cum) = 0.637; Pareto scaling: R 2 X(cum) = 0.953 and Q 2 (cum) = 0.936) (Figure 2A,D). Next, Partial Least Squares-Discriminant Analysis (PLS-DA) was engaged, using also the UV and Pareto scaling ( Figure 2B,E). Analysis resulted to clear separation of the two cell groups, with generated models showing very good predictive ability (UV scaling: R 2 X(cum) = 0.711, Q 2 (cum) = 0.994, 2 components; Pareto scaling: R 2 X(cum) = 0.953, Q 2 (cum) = 0.998, 2 components). Validation of PLS-DA models was performed, allowing of 100 permutations, in order to assess the prediction power of the classification ( Figure  2G,H). Results proved that obtained generated models were not produced randomly and, therefore, they were considered as valid. Finally, Orthogonal Partial Least Squares-Discriminant Analysis (OPLS-DA) led to an equally good separation of the two melanoma cell lines, in both scaling methods (UV scaling: R 2 X(cum) = 0.711, Q 2 (cum) = 0.994; Pareto scaling: R 2 X(cum) = 0.953, Q 2 (cum) = 0.998) ( Figure 2C,F).
Distant features towards the extremes were selected from the S-plot (OPLS-DA model with Pareto scaling) ( Figure 2I), while those exhibiting Variable Importance in Projection (VIP) > 1 were considered as responsible for the separation between the two cell groups ( Table 1). The identified features, together with their corresponding statistical-test information and boxplots, are shown in

Chemometric Analysis of NMR Data
Multivariate statistical analysis was performed, by employing SIMCA P 14.1. In the first step, Principal Components Analysis (PCA) was applied, in order to visualize the obtained data, with the aid of recognizing possible trends and locating possible outliers. Both Unit Variance (UV) and Pareto scaling exhibited clear separation between the two cell lines, WM115 and WM2664, without the presence of outlying values (UV scaling: R 2 X(cum) = 0.712 and Q 2 (cum) = 0.637; Pareto scaling: R 2 X(cum) = 0.953 and Q 2 (cum) = 0.936) (Figure 2A,D). Next, Partial Least Squares-Discriminant Analysis (PLS-DA) was engaged, using also the UV and Pareto scaling ( Figure 2B,E). Analysis resulted to clear separation of the two cell groups, with generated models showing very good predictive ability (UV scaling: R 2 X(cum) = 0.711, Q 2 (cum) = 0.994, 2 components; Pareto scaling: R 2 X(cum) = 0.953, Q 2 (cum) = 0.998, 2 components). Validation of PLS-DA models was performed, allowing of 100 permutations, in order to assess the prediction power of the classification ( Figure 2G,H). Results proved that obtained generated models were not produced randomly and, therefore, they were considered as valid. Finally, Orthogonal Partial Least Squares-Discriminant Analysis (OPLS-DA) led to an equally good separation of the two melanoma cell lines, in both scaling methods (UV scaling: R 2 X(cum) = 0.711, Q 2 (cum) = 0.994; Pareto scaling: R 2 X(cum) = 0.953, Q 2 (cum) = 0.998) ( Figure 2C,F).
Distant features towards the extremes were selected from the S-plot (OPLS-DA model with Pareto scaling) ( Figure 2I), while those exhibiting Variable Importance in Projection (VIP) > 1 were considered as responsible for the separation between the two cell groups ( Table 1). The identified features, together with their corresponding statistical-test information and boxplots, are shown in Table 2 and Figure 3, respectively, with most of them presenting statistically significant differences between the two melanoma cell groups, WM115 and WM2664, herein examined.  Table 2 and Figure 3, respectively, with most of them presenting statistically significant differences between the two melanoma cell groups, WM115 and WM2664, herein examined.     Rapid and reliable identification of human pathology-discrimination metabolites is an issue of major importance in metabolomics, and, thus, heat-maps have the advantage of a direct and systemic overview of the data-matrix analyzed. Hence, NMR data heat-map illustrates an excellent sample clustering and indicates a collection of novel mechanistic biomarkers for human pre-metastatic and Rapid and reliable identification of human pathology-discrimination metabolites is an issue of major importance in metabolomics, and, thus, heat-maps have the advantage of a direct and systemic overview of the data-matrix analyzed. Hence, NMR data heat-map illustrates an excellent sample clustering and indicates a collection of novel mechanistic biomarkers for human pre-metastatic and metastatic melanoma (Figure 4), with significant value in disease prognosis, diagnosis and therapy in the clinic.  An alternative visualization of the statistical results of "omics" data is the Volcano plot, which allows the identification of metabolites that exhibit the strongest combination of Fold Change (FC) (x) and statistical significance (p value). The Volcano diagram is two-dimensional, with the x-axis representing the variation ratio of each metabolite (log2FC) and the y-axis the p value obtained from the corresponding statistical test t (-log10P) [17,18] (Table 2). Metabolites with the most significant contribution to the separation of the two melanoma cell lines, WM115 and WM2664, are highlighted in color, with those in the left segment being elevated in WM115 (green) and those being identified in the right segment presenting an increase in WM2664 (blue) ( Figure 5). An alternative visualization of the statistical results of "omics" data is the Volcano plot, which allows the identification of metabolites that exhibit the strongest combination of Fold Change (FC) (x) and statistical significance (p value). The Volcano diagram is two-dimensional, with the x-axis representing the variation ratio of each metabolite (log2FC) and the y-axis the p value obtained from the corresponding statistical test t (-log10P) [17,18] (Table 2). Metabolites with the most significant contribution to the separation of the two melanoma cell lines, WM115 and WM2664, are highlighted in color, with those in the left segment being elevated in WM115 (green) and those being identified in the right segment presenting an increase in WM2664 (blue) ( Figure 5). The Venn diagram of Figure 6 summarizes the metabolites that were found to be important for separation of the two cell lines, WM115 and WM2664, by multivariate (S-plot, VIPs) and univariate statistical analysis (Volcano plot). Metabolites that emerged from both statistical approaches are presented at the intersection of both sets, and, thus, may serve as mechanistic biomarkers, with strong prognostic, diagnostic and therapeutic potential, for human advanced melanoma undergoing metastasis.  The Venn diagram of Figure 6 summarizes the metabolites that were found to be important for separation of the two cell lines, WM115 and WM2664, by multivariate (S-plot, VIPs) and univariate statistical analysis (Volcano plot). Metabolites that emerged from both statistical approaches are presented at the intersection of both sets, and, thus, may serve as mechanistic biomarkers, with strong prognostic, diagnostic and therapeutic potential, for human advanced melanoma undergoing metastasis. The Venn diagram of Figure 6 summarizes the metabolites that were found to be important for separation of the two cell lines, WM115 and WM2664, by multivariate (S-plot, VIPs) and univariate statistical analysis (Volcano plot). Metabolites that emerged from both statistical approaches are presented at the intersection of both sets, and, thus, may serve as mechanistic biomarkers, with strong prognostic, diagnostic and therapeutic potential, for human advanced melanoma undergoing metastasis.

UHPLC-HRMS/MS
Representative base peak intensity chromatograms of WM2664 cells, in the positive and negative ion mode, are shown in Figure 7.

UHPLC-HRMS/MS
Representative base peak intensity chromatograms of WM2664 cells, in the positive and negative ion mode, are shown in Figure 7.

MS Spectra Processing
The results obtained as a CSV file, after pre-processing of the mass spectra with MZmine 2.31, were evaluated based on the presence of internal standards in all samples and the repeatability of their integration regions, in order to select the appropriate algorithm parameters. For positive Electrospray Ionization (ESI), precursor ions of all three internal standards were detected, whereas in negative ionization the precursor ions for the two internal standards and the Cl adduct for 2-aminophenol were detected. Subsequently, the data were submitted into SIMCA 14.1 and PCA analysis was performed, for all the 20 distinct samples (10 per cell line) and the 9 replicates of Quality Control (QC) sample ( Figure S1). Scores plot exhibits that repetitions of the QC sample do not follow clustering, therefore suggesting systematic error due to instrumental drift, a serious complication that is often encountered in mass spectrometry. Thereby, the QC-robust spline batch correction (QC-RSC) method of normalization, as implemented within the online platform MetaX, was suitably applied, in order to reduce systematic error between measurements, so that only biological differences between examined samples were highlighted [19,20]. Selection of the normalization algorithm is a critical parameter in data analysis. To evaluate results for normalization algorithms and find the most appropriate one, they were, next, analyzed by PCA, using SIMCA 14.1. The more closely arranged the QC samples are shown, the more accurate the normalization method is considered. The Cross-contribution Compensating Multiple standard Normalization (CCMN) algorithm presented the best results in clustering of control samples and was, therefore, chosen as the most suitable normalization method. Figure S1A,B describes the diagram of PCA in UV scaling, for negative ionization, before and after normalization.

Chemometric Analysis of MS Data
After the initial data pre-processing, multivariate analysis (PCA, PLS-DA and OPLS-DA) was applied, in order to explore the underlying metabolic pathways that are mechanistically related to melanoma metastasis. The PCA generated models are presented as scores plot in Figure 8A,D, for the negative (UV and Pareto scaling, respectively) and in Figure 9A,D, for the positive ionization mode (UV and Pareto scaling, respectively). For PCA analysis, R 2 and Q 2 were 0.502 and 0.207, for UV, and 0.447 and −0.009, for Pareto scaling, respectively, in the negative ion mode, whereas in the positive ionization, R 2 X(cum) and Q 2 (cum) were 0.433 and 0.13, for UV, and 0.585 and 0.016, for Pareto scaling, respectively. The 9th sample of WM2664 cell group appears as an outlier (negative and positive ionization), as it can be observed in PCA diagrams, but it is not rejected, as being within the acceptable 5% of the total samples, at the 95% confidence level. Based on the PCA approach, the overall UHPLC-HRMS analysis, after data normalization, is considered reliable for statistical evaluation, since the control samples follow a close clustering pattern.   By applying the supervised PLS-DA methodology, a clear separation between the WM115 and WM2664 cell groups was observed, as it can be nicely described in the corresponding diagrams for negative and positive ion modes, in Figures 8 and 9. The Q 2 (cum) value was 0.966 for the negative ionization, while, for the positive ionization, was 0.865, choosing the 2 main components in UV scaling. The OPLS-DA supervised model used, for both negative (R 2 X(cum) and Q 2 (cum) were 0.427 and 0.958, respectively) and positive ionization (0.369 and 0.856, respective values), showed excellent fitting and predictive capacity.
Furthermore, to confirm the validity of the obtained model, permutation testing was conducted. Not any random model was produced, exhibiting higher R 2 and Q 2 values than the unpermutated, respective, ones ( Figures 8G,H and 9G,H), thus proving the validity of the obtained PLS-DA model. Next, in order to identify those features that contribute the most to the separation, the S-plot was employed, based on the OPLS-DA model, for both ionization modes, engaging SIMCA 14.1. S-plot features that lie remotely from the main-data cluster were selected, as they are considered to contribute the most, with high reliability, to the group separation ( Figures 8I and 9I). Features in the lower-left quadrant correspond to the WM115 cell group, while features in the upper-right quadrant belong to the WM2664 cell group. To verify, the statistically obtained from S-plot, results, t-test and False Discovery Rate (FDR) values [21,22] were calculated, using Graphpad Prism 7. A Q value of 5% By applying the supervised PLS-DA methodology, a clear separation between the WM115 and WM2664 cell groups was observed, as it can be nicely described in the corresponding diagrams for negative and positive ion modes, in Figures 8 and 9. The Q 2 (cum) value was 0.966 for the negative ionization, while, for the positive ionization, was 0.865, choosing the 2 main components in UV scaling. The OPLS-DA supervised model used, for both negative (R 2 X(cum) and Q 2 (cum) were 0.427 and 0.958, respectively) and positive ionization (0.369 and 0.856, respective values), showed excellent fitting and predictive capacity.
Furthermore, to confirm the validity of the obtained model, permutation testing was conducted. Not any random model was produced, exhibiting higher R 2 and Q 2 values than the unpermutated, respective, ones ( Figure 8G,H and Figure 9G,H), thus proving the validity of the obtained PLS-DA model. Next, in order to identify those features that contribute the most to the separation, the S-plot was employed, based on the OPLS-DA model, for both ionization modes, engaging SIMCA 14.1. S-plot features that lie remotely from the main-data cluster were selected, as they are considered to contribute the most, with high reliability, to the group separation ( Figures 8I and 9I). Features in the lower-left quadrant correspond to the WM115 cell group, while features in the upper-right quadrant belong to the WM2664 cell group. To verify, the statistically obtained from S-plot, results, t-test and False Discovery Rate (FDR) values [21,22] were calculated, using Graphpad Prism 7. A Q value of 5% was used, for both positive and negative ionization modes. As potential biomarkers were accepted the ones whose S-plot features proved to significantly differ between cell groups, according to t-test and FDR calculations. Due to the large number of features that met the aforementioned criteria (86 for negative and 17 for positive ionization), a smaller range of features was selected for identification (12 for negative and 15 for positive ionization) (Table S2).

Identification of Selected MS Features
The R statistical package xMSannotator [23,24], as well as the METLIN [25] and HMDB [16] databases, were suitably employed to identify the selected features. Presence of adducts and MS/MS fragment ions, consistent with METLIN and HMDB data for each particular metabolite, were examined, and the isotope ratio was also calculated, for every identified feature. The following scoring scale was used to evaluate the reliability of each assignment: two points for each parent ion or adduct or pseudo MS/MS fragment, 1.5 points for each MS/MS fragment (as it has been acquired in low-resolution mode, using parallel scans to Linear Ion Trap) and 0.5 points for isotopic ratio showing Relative Standard Deviation < 30% (RSD < 30%) when compared to the theoretically computed one by the Thermo Xcalibur Qual Software Browser [26]. Assignments with a sum of less than or equal to four points were considered to be sufficiently reliable. Table 3 presents the melanoma-derived metabolites identified in negative and positive ionization modes. Identification data, derived from both NMR Spectroscopy and Mass Spectrometry, revealed metabolites critically implicated in major (canonical), or alternative (non-canonical), biochemical pathways. Therefore, targeted screening for the remaining metabolites of these pathways, as well as for other selected metabolites that may essentially contribute to human carcinogenesis, could provide important additional information regarding the human skin oncometabolome. This approach primarily focused on "rare" amino acids, urea cycle, arginine metabolism and cycle of Lynen (ketone bodies). Presence of parent ions, and their respective Na, K and NH 4 adducts, together with corresponding dimers, for all the selected metabolites, were carefully considered in both ionization modes.
Targeted screening was performed with the MZmine 2.31 software and the parameters used are listed in Table S3. Metabolites identified in the negative (Table 4), or positive (Table 5) ionization mode were processed with the same manner as the one chosen for those of untargeted screening.
Boxplots and ROC curves, for the most significant metabolites recognized from MS data, are shown in Figure 10 (negative ionization) and Figure 11 (positive ionization).

Discussion
Metabolic rewiring is a hallmark of malignant transformation in several types of cancer [3,27,28]. Tumor microenvironment, which is often not sufficiently vascularized, compels tumor cells to oxygen lack and nutrient deprivation, orchestrating metabolic changes that promote survival and growth in

Discussion
Metabolic rewiring is a hallmark of malignant transformation in several types of cancer [3,27,28]. Tumor microenvironment, which is often not sufficiently vascularized, compels tumor cells to oxygen lack and nutrient deprivation, orchestrating metabolic changes that promote survival and growth in stressful conditions. Rapid proliferation increases the demand for energy supply that is required for macromolecular biosynthesis, as well as synthesis of amino acids, nucleotides and lipids.
One of the major alterations detected in the present study is the upregulated glycolysis and the increased production of lactic acid, likely derived from pyruvic acid to prevent its oxidation in mitochondria even in the presence of oxygen. This phenomenon was, first, described in 1920 by Otto Warburg and is called the "Warburg effect" [29,30]. Alternatively, enhanced glycolysis may likely result in accumulation of the methylglyoxal (MG) cytotoxic metabolite that is spontaneously formed by the degradation of two glycolytic intermediates [31]. Hence, cancer cells, in order to avoid MG-induced self-destruction, and as a survival strategy, upregulate the glyoxalase-dependent enzymatic machinery to remove MG and convert it into D-lactic acid, thus causing its cellular accumulation [31]. Interestingly, overexpression of glyoxalases has been observed in human melanoma [32]. Beyond this metabolic reprogramming, and since glucose is converted almost quantitatively to lactic acid, the TCA (tricarboxylic acid) cycle functions, regardless of oxygen presence, in the reverse direction, using glutamine as an alternative carbon source [10]. Elevated glycolysis and glutaminolysis could explain the absence of both glucose and glutamine metabolites from the two cell lines, WM115 and WM2664, but most importantly the significant increase in lactic acid (FC NMR = 3.27) and glutamic acid (FC NMR = 4.84; FC MS = 5.81), as well as in the TCA cycle intermediates fumaric acid (FC NMR = 2.01) and succinic acid (FC NMR = 4.82; FC MS = 4.04) [11], in the metastatic melanoma cell line WM2664. A rise in the (metabolite-linker) alanine content [10] from implementation of these two processes is expected, since alanine aminotransferase catalyzes the transfer of an amino group from glutamic to pyruvic acid, with alanine and oxoglutaric acid being the reaction products. Glutamine-derived oxoglutaric acid supplies the TCA cycle with additional carbon, while the excess nitrogen is secreted in the form of alanine. In addition to the approximately 2-fold (x) increase of the detected alanine (FC NMR = 1.80) in the metastatic cells, a similar elevation was also observed for proline (FC NMR = 2.01), whose enhanced generation is a characteristic metabolic change in melanoma [33,34], with rather glutamine than arginine serving as a key metabolic precursor. Synthesis of arginine in melanoma is often suppressed due to lack of argininosuccinate synthetase (ASS) [35] (an enzyme involved in the urea cycle), which could justify the absence of arginine from both cell lines.
The glycolysis flux may also be diverted from the intermediate 3-phosphoglyceric acid to serine synthesis, a process that seems upregulated in melanoma [36]. A significant percentage of serine content (FC NMR = 1.02) produced is converted to glycine (FC NMR = 1.49), the uptake, synthesis and consumption of which have been previously associated with tumorigenesis [37]. Creatine is formed from glycine by two consecutive reactions and, in turn, can be phosphorylated into phosphocreatine, ultimately generating adenosine triphosphate (ATP) [38]. The elevated concentration of creatine (FC NMR = 4.91) and phosphocreatine (FC NMR = 2.53) in WM2664 cells mainly reflects their essential roles in energy storage required under the increased metabolic rates of tumor cells. In addition to synthesis of ATP, these two metabolites are also involved in its trafficking from the internal mitochondrial membrane to the cytoplasm, thus acting as energy distribution regulators [38,39]. However, ATP, the most important molecule in the cell-energy network and homeostasis, is predominantly formed during the glycolysis process, as well as by oxidative phosphorylation in mitochondria. The energy requirements of cancer cells are clearly increased and they rise proportionally to the rate of tumor growth, hence justifying the significantly elevated ADP/ATP contents (FC NMR = 7.04) in the WM2664 metastatic melanoma cells (as compared to the WM115 cell line of reference).
Enhanced glycolysis in melanoma can result in increased glucose flow to the hexosamine biosynthetic pathway (HBP), with UDP-N-acetylglucosamine (UDP-GlcNAc) being the final product, a high-energy donor of the GlcNAc group of protein substrates, which was herein proved as an elevated critical metabolite in the WM2664 cell line (FC NMR = 2.04). It has been previously reported that the levels of β-1,6-branched N-glycans are often increased in malignant tumors and this rise is closely related to tumor progression [40,41]. Interestingly, incubation of B16 melanoma cells with GlcNAc upregulates β-1,6-branched oligosaccharide levels in proportion to the elevation of intracellular UDP-GlcNAc [42].
Another glycolysis product herein shown to be notably elevated in WM2664 cells is myo-inositol (FC NMR = 7.76). This metabolite is intrinsically synthesized from glucose in three steps: (a) glucose is phosphorylated by hexokinase, (b) produced glucose 6-phosphate is, next, converted to myo-inositol 1-D-phosphate (by the MIPS synthase) and (c) myo-inositol 1-D-phosphate is finally dephosphorylated towards free myo-inositol [43]. Myo-inositol plays an essential role in regulating cell volume and osmotic pressure, while in vitro and in vivo NMR studies unveil correlation with tumor-cell density (its concentration is often increased) [9]. Myo-inositol addition to B16 melanoma cells, having been pre-treated with LiCl, partially inhibited the anti-proliferative and morphological effect of LiCl on these cells [44].
Taurine, a metabolite that is biosynthesized via the cysteinesulfinic acid pathway [45], exhibits metabolic activity and role similar to myo-inositol ones in cancer, while it seems increased in WM2664 metastatic melanoma cells (FC NMR = 1.52). Although it is thought to possess anti-oxidant properties, its elevated contents may highlight the importance of oxidative stress in metastatic melanoma and in disrupted balance of metabolic pathways [46].
The metabolite that demonstrated the most striking difference between the two cell lines, WM115 and WM2664, was hypoxanthine, whose concentration was measured at least 14-fold (x) higher in WM2664 (FC NMR = 14.05; FC MS = 94.26), as compared to WM115 (reference) cells. Metabolic studies for gastric cancer and melanoma have described higher plasma and tissue levels, due to higher tumor-cell proliferation rates, conditions under which hypoxanthine acts as a substrate and source of nitrogen [11,47]. Hypoxanthine is formed via the Purine Biosynthesis Pathway (PBP) [48] that is characterized by the conversion of ribose 5-phosphate to inosine monophosphate (FC NMR = 1.45), which is, then, dephosphorylated to inosine, the major hypoxanthine precursor molecule. Thereby, the remarkably increased levels of hypoxanthine are tightly linked to the significantly decreased contents of inosine in the WM2664 metastatic melanoma cells (FC NMR = 0.49; FC MS = 0.48). It seems that primary melanoma favors the metabolic route of inosine monophosphate to guanosine (FC MS = 0.42), a molecule that is notably elevated in the WM115 cells and has been suggested to serve as biomarker [14].
Metabolites that are also presented with a statistically significant increase in the primary melanoma WM115 cell group were the choline (FC MS = 0.37), phosphocholine (FC MS = 0.41) and glycerophosphocholine (FC MS = 0.73), which are key components of cell membranes [9]. Abnormal choline metabolism has proved as a critical metabolic feature being tightly associated with aberrant membrane synthesis and tumorigenesis. In the majority of malignant cases, choline is extensively metabolized to phosphocholine [49]. The observed reduction of phosphocholine levels in WM2664 cells likely derives from augmentation of lipid metabolism in metastatic melanoma [50] and/or from acidic environment(s) of the tumor limiting its formation [51].
In contrast to the rest of the herein identified amino acids, leucine (FC NMR = 0.87), isoleucine (FC NMR = 0.91) and, to a lower extent, valine (FC NMR = 0.95) were presented with decreased contents in the WM2664 cell line. Their degradation leads to production of acetyl-coenzyme A, which is, then, oxidized via the TCA cycle. Branched-chain amino acids (BCAAs) are essential components of proteins, mainly due to their hydrophobicity, and, therefore, in oncogenic conditions of increased protein synthesis and energy requirement, their consumption and oxidation are often intensified [52]. Although results from other metabolic cancer studies converge to the experimental reduction of BCAAs [11,49], differentiation of their levels has not proved consistent [50].
Interestingly, a metabolite that also exhibited a decrease in the WM2664 metastatic melanoma cells, and was identified by non-targeted MS analysis, was the 7-hydroxy-6-methyl-8-ribityl lumazine [53] (FC MS = 0.53), which seems to be implicated in the action of melanosomal tyrosinase [54]. Melanosome constitutes an organelle for melanin synthesis, storage and transfer, but, despite the number of reports, its mechanistic direct correlation with melanoma and metastasis requires further investigation. Since 7-hydroxy-6-methyl-8-ribityl lumazine has been previously proposed as biomarker for prostate cancer [55], it could similarly serve as melanoma biomarker, with strong prognostic, diagnostic and therapeutic value.
The metabolic pathways that are herein unearthed to most likely control metastasis in human melanoma are illustrated in Figure 12 (metabolites with statistically significant increase in the WM2664 metastatic melanoma cells are shown in red, while the ones significantly reduced in blue).

Int. J. Mol. Sci. 2020, 21, x FOR PEER REVIEW
The metabolic pathways that are herein unearthed to most likely control metastasis in human melanoma are illustrated in Figure 12 (metabolites with statistically significant increase in the WM2664 metastatic melanoma cells are shown in red, while the ones significantly reduced in blue).

Figure 12.
Flow-diagram showing pathways and cognate metabolites critically involved in the discrimination of primary and metastatic human melanoma cells. Blue: reduced levels (Fold Change < 0.50); Red: increased levels (Fold Change > 2). The metastatic cell group is presented with elevated contents of hypoxanthine, myo-inositol, AXP (X: D, or T), UDPs, amino acids and organic acids, and decreased levels of guanosine, inosine and cholines, thus indicating the perturbed purine, pyrimidine and amino acid metabolism during metastasis in melanoma.

Cell Lines and Culture Conditions
The human melanoma cell lines WM115 (primary tumor) and WM2664 (lymph-node metastases), which derived from the same (female; 58 years old) patient, were purchased from ECACC/Sigma-Aldrich (Munich, Germany), and they were cultured in 1× DMEM medium, being supplemented with fetal bovine serum, L-glutamine, sodium pyruvic, sodium bicarbonate, non-essential amino acids, penicillin and streptomycin, at 37 • C, 5% CO 2 , and > 95% humidity.

Cell Collection and Storage
For each cell line, 10 7 melanoma cells, grown in high-density (high-confluence) cultures, produced samples of the required high quantity and quality for metabolomics studies. Culture dishes were placed on ice, medium was aspirated and cells were washed twice with ice-cold 1× PBS to remove residual medium traces. Ice-cold 1× PBS, up to the minimum necessary volume, to cover the culture plate, was added, for immediate quenching of cell metabolism. Specifically, cells were harvested through a mild scrapping process and the suspensions (1× PBS and (adherent, and suspended) cells) were centrifuged at 550× g for 5 min at 4 • C. Supernatants were aspirated and cell pellets were stored at −80 • C, for further metabolomics processing. Ten samples from each cell line, WM115 and WM2664, were simultaneously analyzed.

Extraction Protocol
Metabolites were extracted according to Wu's stepwise protocol [56]. Briefly, 200 µL of ice-cold methanol and 200 µL of chloroform were added to the frozen, metabolically quenched, cells, and the suspension was sonicated for 5 min on ice. Next, 180 µL of ddH 2 O was added and vortexing for 1 min followed. Processed samples remained on ice for 15 min, centrifuged at 16,000 rpm for 20 min at 4 • C and the upper aqueous phase together with the lower organic phase were collected. The whole procedure was repeated once more (except from the sonication step) and the two phases, after being suitably collected, were added to the corresponding ones derived from the first extraction step.

Sample Preparation
Polar extracts were dried in a centrifuge evaporator, and reconstituted in phosphate buffer pH = 7.4 (containing TSP as internal standard and NaN 3 as preservative) and deuterated water (700 µL, 10:90 (v/v), Buffer:D 2 O), for the preparation of NMR samples. 100 µL of each NMR sample were dried under vacuum and washed twice with water:acetonitrile, 95:5 (v/v), in order to exchange the deuterated protons for MS analysis. Dried samples were reconstituted in 100 µL of water:acetonitrile, 95:5 (v/v), containing 1% of internal standards. A pooled sample was prepared by mixing aliquots (16 µL) of each individual sample and was used as Quality Control (QC) sample, for mass analysis performance control.

NMR Analysis
All NMR experiments were performed using the NMR Bruker AVANCE III 600 MHz Spectrometer at 300 K, equipped with a PABBI z-gradient probe. After the optimization of acquisition conditions, a fully automated procedure was accomplished through the ICONNMR program suite (TOPSPIN 2.1 version, Bremen, Germany), controlling a B-ACS 60 system (Bruker BioSpin GmbH, Bremen, Germany). 1 H 1D NMR spectra were acquired with 256 scans of 64 k data points, over a spectral width of 12,335.5 Hz. The acquisition time was 2.65 s, the relaxation delay was 9.0 s and the mixing time was 0.01 s. Solvent suppression of residual water was achieved using the standard NOESY (Nuclear Overhauser Effect Spectroscopy) pre-saturation pulse sequence (noesygppr1D, Bruker library). 2-D J-resolved spectra were acquired with 4 scans of 8 k data points over a spectral width of 10,000 Hz, resulting in an acquisition time of 0.4 s; the relaxation delay was 1 s. For a selected sample, WM2664 cell line, homonuclear (Total Correlation Spectroscopy, TOCSY) and heteronuclear (Heteronuclear Single Quantum Coherence Distortionless Enhancement by Polarization Transfer, HSQC-DEPT 135) 2D spectra were recorded to aid the unambiguous assignment of metabolites. The TOCSY experiment was performed with 128 scans while 2 k data points were collected over a spectral width of 9615.4 Hz. The acquisition time was 0.1 s and the relaxation delay was 1.5 s. HSQC-DEPT 135 experiment was performed for a spectral width of 9615.4 Hz for proton frequencies and 27,164 Hz for 13 C; 320 scans were collected.

MS Analysis
The UHPLC-ESI-LTQ ORBITRAP Discovery System (Thermo Fisher Scientific, Rockford, IL, USA) was used for data acquisition, and the Xcalibur 2.1 software was employed for the subsequent visualization, and the initial MS analysis. For chromatographic separation of the metabolome, a C 18 Hypersil Gold column (100 × 2.1 mm, 1.7 µm) (Thermo Fisher Scientific, Rockford, IL, USA) and a pre-column Waters VanGuard  Table S4A. Analysis was performed using the Fourier Transform Mass Spectrometry mode of the LTQ Orbitrap (FTMS), in the full-scan ion mode, applying a resolution of 30,000, while acquisition of mass spectra was performed in every case using the centroid mode. Data-dependent acquisition (DDA) capability has been also engaged, allowing for MS/MS fragmentation of the 3 most intense ions.

NMR
Pre-processing of the NMR spectra (data reduction and bucketing) was conducted using the AMIX Ver. 3.9.12 suite (Analysis of MIXtures, Bruker BioSpin, GmbH). 1 H NMR spectra from 9.40 to 0.60 ppm were segmented into 0.02 and 0.001 ppm bucket widths. Regions, including residual water signal (δ 5.50-4.631 ppm) and contaminations observed in blank samples (δ 1.260-1.241 ppm and δ 2.250-2.221 ppm), were excluded from the analysis [57][58][59]. When normalization to total intensity is performed, alterations of major metabolites may cause biased results for low concentration metabolites [60]. Hence, in order to avoid biased results, and since there was no need to correct for biological or analytical variability, we decided not to perform any normalization. NMR bins, prior to multivariate analysis, were scaled to Unit Variance (UV) and Pareto and the results were compared.

MS
Mass Spectrometry (MS)-data pre-processing [19] was performed by employing MZmine 2.31 [61], Microsoft Excel 2010 and the MetaX [62] online platform. Initially, raw data were introduced into MZmine 2.31, via which noise was filtered, while peaks were detected, using the ADAP Chromatogram Builder algorithm [63], and extracted-ion chromatograms (EICs) were obtained for each sample tested. Next, peak de-convolution was carried out with the Wavelets ADAP algorithm, a continuous wavelet transform (CWT)-based peak-detection algorithm, which uses shape information based on calculated coefficients and, consequently, reduces false positives [64]. Subsequently, isotopic distribution search was performed, with the "isotopic-peak grouper" and chromatographic peaks being aligned via employment of the Join Aligner algorithm, according to their elution time. Finally, after filling in the gap entries, the CSV file was exported. The final Table was edited in Microsoft Excel 2010, using the "concatenate", "round" and "transpose" commands. Data were normalized based on QC samples, prior to statistical analysis, with the QC-RSC algorithm of the MetaX platform (Table S4B).

Statistical Analysis
Principal Component Analysis (PCA), Partial Least Square-Discriminant Analysis (PLS-DA) and Orthogonal Partial Least Square-Discriminant Analysis (OPLS-DA) were herein engaged, as implementation via the SIMCA P 14.1 software (Umetrics, Umea, Sweden). Chemometric models were validated, using permutation tests with 100 random permutations. In order to pinpoint the most influential variables for group separation, selection of those with high-magnitude and -reliability of the S-plot derived from the OPLS-DA approach and the corresponding Variable Importance in Projection (VIP) values were taken into account, following a combinational fashion. The GraphPad Prism 7 software was used for the calculation of AUROC (Area Under the Receiver Operating Characteristic) curves and following a series of univariate statistical tests for each identified metabolite ((a) normality test (D' Agostino and Pearson), (b) outlier identification test (Grubbs) and (c) heteroscedasticity test (F-test)) were conducted prior to Student's t-test, or Welch test. Multiple hypothesis testing for false positive results, were minimized with Bonferroni-Dunn and, the less strict, False Discovery Rate (FDR) corrections. Fold Change (x) (FC = (mean WM2664 )/(mean WM115 )) and z-score ((mean WM2664 − mean WM115 )/SD WM115 ) were also calculated, using Microsoft Excel 10.

Conclusions
Melanoma accounts for more than 60% of all deaths associated with skin cancer, thus necessitating its early and accurate detection to improve treatment outcomes. Given the pivotal role of metabolic reprogramming in cancer initiation and progression, identification of specific metabolic patterns that can differentiate pre-metastatic from metastatic tumor cells may serve as a valuable biochemical tool that carries strong prognostic, diagnostic and therapeutic potential.
Using Nuclear Magnetic Resonance (NMR) Spectroscopy and High-Resolution Mass Spectrometry (HRMS) Coupled with Ultra High Performance Liquid Chromatography (UHPLC-HRMS), as part of a non-targeted metabolic approach based on cell activity, metabolomes of the human melanoma cell lines WM115 (primary) and WM2664 (metastatic), both derived from the same patient, were herein profiled. Next to pre-processing of the original data, a number of critical metabolites were identified, via suitable and successful engagement of the ChenomX NMR software, and the HMDB and METLIN libraries for MS and MS/MS spectra, while univariate and multivariate statistical analyses were performed using the SIMCA 14.1 and GraphPad software Prism 7 platforms.
A collection of 50 metabolites were detected through the NMR analysis, 40 of which presented a statistically significant difference between the pre-metastatic WM115 and metastatic WM2664 cells of human melanoma. Moreover, 30 out of 36 metabolites proved to differ significantly between the two respective metastatic states (Tables 3-5) via the UHPLC-HRMS analysis (10 common metabolites with NMR). Differentially synthesized (or consumed/degraded), secreted and/or sequestered metabolites are components of a diverse set of fundamental metabolic pathways, such as the glucose catabolism, glutamine catabolism and choline metabolism, which seem to critically contribute to the survival and growth of advanced melanoma cells. Interestingly, hypoxanthine, ATP/ADP/AMP, myo-inositol, creatine, glutamic acid, succinic acid, NADH, lactic acid, creatine phosphate, UDP-N-acetylglucosamine, proline and fumaric acid were found significantly elevated in WM2664, whereas choline, o-phosphocholine, inosine, guanosine and leucine were presented with notably reduced contents in the same cell line, and, therefore, they could all be considered as powerful tumor biomarkers, carrying strong prognostic, diagnostic and therapeutic values for the metastatic form(s) of human advanced melanoma(s). Mapping of novel metabolic landscapes in melanoma cells that grow in vitro and bear distinct metastatic activities is expected: (a) to illuminate our knowledge and understanding regarding the Molecular Biology of advanced melanoma, and (b) to open new therapeutic windows for the metastatic disease form(s) in the clinic.