Ammonia solution 0.88 SG, formic acid, Optima^TM LC/MS grade water, and Optima^TM LC/MS grade acetonitrile were purchased from Fisher Scientific (Loughborough, United Kingdom). Ammonium acetate was from Sigma-Aldrich (Gillingham, United Kingdom), and ammonium formate from Thermo Fisher Acros Organics (Geel, Belgium). (Z)-4-Hydroxytamoxifen (4-OHT) and MS222 (ethyl 3-aminobenzoate methanesulfonate salt, also called tricaine) were purchased from Sigma Aldrich (Gillingham, United Kingdom).

Adult zebrafish were maintained in the Bioresearch and Veterinary Services (BVS) Aquatic facility in the Queen’s Medical Research Institute, the University of Edinburgh. Housing conditions were similar to those described in the Zebrafish Book (Westerfield, 2000) with 14/10 h light/dark cycle and water temperature of 28.5°C. The water quality, hardness and conductivity were controlled daily, and the pH was maintained at 7–7.5. All experiments were conducted with local ethical approval from the University of Edinburgh and in accordance with United Kingdom Home Office regulations (Guidance on the Operations of Animals, Scientific Procedures Act, 1986).

The zebrafish lines used in this study were generated in house by the Feng lab using Tol2 mediated transgenesis. Tg(krtt1c19e::KalTA4-ER^T2; cmlc2::eGFP)^ed201 (hereafter K19 Gal4 driver) drives basal keratinocytes (Lee et al., 2014) expression of the KalTA4-ER^T2 transcription factor, together with cmlc2::eGFP for identification of transgene-containing larvae by GFP fluorescence in the myocardium (green heart). Tg(UAS::eGFP-HRAS^G12V; cmlc2::eGFP)^ed203 (hereafter UAS:RAS) carries an UAS driven eGFP-HRAS^G12V expressing cassette and a cmlc2 promoter driven eGFP expression cassette as selection marker. Tg(UAS::eGFP-CAAX; cmlc2::eGFP)^ed204 (hereafter UAS:CAAX) carries an UAS driven eGFP-CAAX expressing cassette and a cmlc2 promoter driven eGFP expression cassette as selection marker.

Adult K19 Gal4 driver fish were set for pair mating with either UAS:RAS or UAS:CAAX fish in breeding tanks that contained a divider. The next morning, the dividers were removed, and the fertilized embryos were collected within the next 2 h to ensure synchronous development. The embryos were maintained in 90 mm petri dishes (maximum of 50 embryos/dish) containing 0.3 × Danieau’s solution (17.40 mM NaCl, 0.21 mM KCl, 0.12 mM MgSO₄7H₂O, 0.18 mM Ca(NO₃)₂, 1.5 mM HEPES), in a 28.5 °C incubator.

The detailed induction protocol is described in Ramezani et al. (2015). In brief, the induction solution consisted of 0.3 × Danieau’s solution with 5 μM 4-OHT and 0.5% DMSO to enhance penetration of 4-OHT. A 10 mM stock of 4-OHT dissolved in 96% ethanol was stored at −20°C protected from light. At 52 h post-fertilization (hpf), larvae were transferred to petri dishes containing 20 mL 4-OHT induction solution, at 50 larvae per dish. The larvae were then maintained at 28.5°C in the dark.

At 22 h post-induction (hpi; corresponding to 74 hpf), the larvae were anaesthetized with MS222 and appropriate phenotypes were screened and sorted. Larvae with 30–60% coverage of GFP positive skin cells were collected as pre-neoplastic cell (PNC) group, and larvae with GFP negative skin were collected as Wild Type (WT) siblings group (Figure 1). During screening, the embryos were placed in fresh 4-OHT induction solution and the sorted larvae were maintained at 28.5°C in the dark for an additional period to yield a total induction time of 24 h.

For metabolite extraction, 10 biological replicates of the PNC group and sibling controls were used. The embryos were anaesthetized, washed twice with 0.3 × Danieau’s solution and transferred into 1.5 mL eppendorf tubes at 25 embryos/tube. Any remaining liquid was carefully removed by pipetting and 150 μL of ice-cold extraction buffer (chloroform:methanol:ddH₂O, 1:3:1; chilled at −20°C prior to use) was added. The samples were placed on ice during handling and incubated at 4°C, 80 rpm for 40 min for extraction with further mixing after 15 and 30 min of incubation. The metabolite extracts had a yellowish color at the end of incubation while larvae remained intact. Samples were vortexed briefly, at high speed and centrifuged (15,000 × g) for 10 min at 4°C. The supernatants were transferred to pre-cooled 1.5 mL eppendorf tubes, which were immediately placed on dry ice. The samples were stored at −80°C until transported for LC-MS-IMS analysis on dry ice.

EdU labeling was carried out using the Click-iT Plus EdU Alexa Fluor 647 Imaging Kit (Life Technologies, C10640, Thermo Fisher Scientific, Loughborough, United Kingdom) as described previously (van den Berg et al., 2019). Briefly, at 62 hpf (10 hpi), 74 hpf (22 hpi), or 86 hpf (34 hpi), 1 nL of 10 mM EdU (5-ethynyl-2′-deoxyuridine, a nucleoside analog of thymidine) was injected into the yolk of the larvae. After 2.5 h incubation at 28.5°C, the larvae were fixed in 4% paraformaldehyde (PFA) for 30 min at room temperature (12 and 24 hpi samples) or o/n at 4°C (36 hpi samples). The larvae were permeabilized in phosphate buffered saline (PBS) containing 0.5% Triton X-100 (PBST) for 5 min three times, and blocked with PBST containing 3% (w/v) Bovine Serum Albumin (Sigma-Aldrich, Gillingham, United Kingdom) for 1 h. This was followed by 30 min incubation with the Click-it Plus reaction cocktail according to manufacturer’s instructions, using 250 μL of reaction cocktail per max 10 larvae. Following EdU labeling, larvae were washed with PBST for 5 min three times and re-blocked with PBST containing 5% (v/v) Goat Serum (Sigma-Aldrich, Gillingham, United Kingdom) for 2 h. eGFP immunostaining was performed with rabbit monoclonal anti-GFP antibody (1:200; cat. 2956, Cell Signaling Technology, London, United Kingdom) and Alexa Fluor 488 Goat anti-Rabbit secondary antibody (1:250; A-11008, Invitrogen, Thermo Fisher Scientific, Loughborough, United Kingdom) as described (van den Berg et al., 2019). Stained larvae were stored at 4°C in a glycerol based antifadent mountant (AF1, CitiFluor, Hatfield, PA, United States) until mounted for imaging.

Confocal imaging was performed using a Leica TCS SP8 AOBS confocal laser scanning microscope attached to a Leica DMi8 inverted microscope, with a 20 × dry lens, and excitation lasers 405, 488, and 633 nm were used. Percentage of EdU positive cells within eGFP positive cells were manually counted using the Imaris 9.0 software. In addition, PNCs were assessed for elongated shape and the presence of filopodia-like structures. Statistical significance was determined by multiple t-test using the Holm-Sidak method (EdU staining) and one-way ANOVA with Tukey’s multiple comparisons test (PNC morphology) using Prism 6 (GraphPad). Z projections of representative image stacks were generated using Fiji software (Schindelin et al., 2012).

4-OH tamoxifen induced zebrafish larvae were screened for fluorescently labeled HRAS and CAAX and anesthetised with MS-222. The two groups of larvae, HRAS and CAAX control were dissociated using collagenase IV (25 mg/mL) in medium containing goat serum (2% v/v), HBSS, HEPES (15 mM) and D-glucose (25 mM). The two samples were incubated for 3 × 5 min at 28.5°C until fully dissociated. The cell suspensions were diluted with an equal volume of medium containing goat serum (10% v/v) in HBSS, HEPES, D-glutamine and kept cold on ice. DAPI was added to distinguish dead cells and the samples were sorted on a BD FACS Aria II instrument using a 100 μm nozzle, selecting for fluorescently labeled HRAS or CAAX populations. Batches of 1,000 cells were collected in Eppendorf tubes containing 4 μL 0.2% Triton X and RNAase inhibitor (Thermo Fisher Scientific). The samples were processed for cDNA synthesis using a modified Smart-seq2 protocol (Picelli et al., 2014). Post clean-up, the cDNA quality was assessed using a LabChip GX instrument. Gene expression was analyzed using LightCycler 480 SYBR Green I Master (Roche) on a LightCycler 96 system using the following primers; sat1a.2_for 5′-CC GTTTTATCACTGCGTCGT-3′ and sat1a.2_rev 5′-AAAGC AGCTTCCCAATCCAG-3′, itpa for 5′-AGCCTGAGGGCTAT GACAAA-3′ and itpa_rev 5′-GTTTTGAGCGCTTGGTCT CT-3′, pnp5b_for 5′-TTCTCAGTCTTCATGCGGGA-3′ and pnp5b_rev 5′-CACCTCATGAACAGTGCTCA-3′, rrm2_for 5′-ACTGGGCACTCAACTGGATT-3′ and rrm2_rev 5′-GTC CCCTCTTCTTTAGCCAGA-3′.

Chromatographic separation was performed using a SeQuant ZIC-pHILIC, 5 μm polymeric 4.6 mm × 150 mm column (Merck KGaA 1.50461.0001, Darmstadt, Germany) coupled to a SeQuant ZIC-pHILIC Guard 20 × 2.1 mm column (Merck KGaA 1.50437.0001, Darmstadt, Germany). Two different solvent systems of low and high pH were used to run a 30-min gradient. The solvent system for acquiring data in the positive ionization mode consisted of 10 mM ammonium formate in water with 0.1% formic acid, pH 3 (solvent A) and 10 mM ammonium formate in water/acetonitrile (1:9) with 0.1% formic acid, pH 3 (solvent B). Similarly, the solvent system for acquiring data in negative ionization mode consisted of 10 mM ammonium acetate in water, pH 9 (solvent A) and 10 mM ammonium acetate in water/acetonitrile (1:9), pH 9 (solvent B). The solvent gradient for both ionization modes consisted of 80% solvent B at the start of the run, which was reduced to 20% solvent B in 15 min, and further to 5% solvent B in 1 min, where it was maintained for 4 min. At 21 min the column was returned to initial conditions of 80% solvent B and maintained as such until 30 min. A constant flow rate of 0.3 mL/min was observed during the run and the column was maintained at 30°C. For analysis, 5 μL sample was injected into the column. A quality control sample generated by pooling equal volumes of each sample was injected five times at the beginning of the experiment to condition the column and after every five samples to monitor the instrument state over the course of data acquisition.

The LC-MS-IMS instrumentation consisted of an Agilent 1290 Infinity II series UHPLC system coupled to an Agilent 6560 IM-qTOF (both Agilent Technologies, Santa Clara, CA) with a Dual Agilent Jet Stream Electron Ionization source. In both ionization modes data was acquired in the 50-1700 m/z range, with an MS acquisition rate of 0.8 frames/s. The nebulizer pressure was set to 60 psi, gas temperature to 225°C and drying gas (N₂) flow rate to 13 L/min. Sheath gas was set to 340°C with a flow rate of 12 L/min, and the instrument was operated at a capillary voltage of 3,000 V, nozzle voltage of 200 V, fragmentor voltage of 395 V, and octupole voltage of 750 V. Ionization mode parameters for the IMS are instrument specific and available upon request from the authors. Instrument calibration and tuning was performed separately for each ionization polarity using the ESI-L low concentration tuning mix from Agilent Technologies (Santa Clara, CA). A reference mass solution consisting of 50 μM ammonium trifluoroacetate, 5 μM purine and 1.125 μM HP-0921 was injected continuously into each sample to recalibrate for accurate mass and drift time during data processing. The ES-TOF reference mass solution kit was purchased from Agilent Technologies (Santa Clara, CA).

Data acquisition and processing were performed using the Agilent MassHunter 10 software suite. Briefly, ion multiplexed data files and calibration files acquired using MassHunter Data Acquisition 10.0 software were demultiplexed using the PNNL PreProcessor v2020.03.23. The default settings for demultiplexing, moving average smoothing, saturation repair and spike removal were applied to the data. The data files were recalibrated for accurate mass and drift time using the AgtTofReprocessUi and the IM-MS Browser 10.0, respectively. Two reference masses, m/z 121.050873 and 922.009798 (positive ionization mode) and m/z 112.985587 and 1033.988109 (negative ionization mode) were used for recalibration.

Molecular features were extracted in Mass Profiler 10.0 with a retention time tolerance of ± 0.3 min, drift time tolerance of ± 1.5% and accurate mass tolerance of ± (5 ppm + 2 mDa). The raw multiplexed data, the reconstructed demultiplexed data and Mass Profiler feature lists (.cef files) were used in the High Resolution Demultiplexer (HRdm) 1.0 beta v41 for further peak deconvolution (May et al., 2020). Molecular features were re-extracted from HRdm files using Mass Profiler 10.0 and annotated using accurate mass and CCS values with the McLean CCS Compendium PCDL (version 20191101; Picache et al., 2019). A CCS value tolerance of ± 1% was applied, and the positive ion species (M+H)⁺, (M+Na)⁺, (M+NH₄)⁺, and the negative ion species (M-H)^– and (M+CH₃COO)^– were searched.

Multivariate statistical analysis and pathway enrichment analysis were performed using the MetaboAnalyst 5.0 web-based platform (Chong and Xia, 2020). The data were log-transformed and auto-scaled. To get an overview of the entire dataset, both annotated and unannotated molecular features were used to generate PLS-DA plots and heatmaps. To perform pathway enrichment analysis, only the annotated features without the relative peak intensities was used. For more detailed analysis of the altered pathways, the list of annotated compounds with their relative intensities was submitted to the pathway analysis tool, log-transformed, auto-scaled and examined against the Danio rerio KEGG pathway library, using global test and relative betweenness centrality methods. The data was visualized as scatter plots. For the pathways identified to be significantly altered between the study groups, the dataset was manually screened to identify the intermediates that were not recognized by MetaboAnalyst 5.0 online platform owing to synonymous annotations. Box and whisker plots for all the identified metabolites were retrieved from the feature details view following PLS-DA analysis on MetaboAnalyst 5.0 online platform, and the normalized values are presented.

The eumelanin biosynthesis pathway intermediates were searched for in the dataset using accurate mass values. The only precursor identified was 5,6-dihydroxyindole-2-carboxylic acid with an accurate mass of m/z 193.03751. To calculate ATP to ADP ratio, the peak intensity values in the annotated dataset were used and relative abundance is presented. The 260 annotated features were classified into compound classes using ClassyFire, the web-based application for automated structural classification of chemical entities (Djoumbou Feunang et al., 2016). The output file was manually curated to unambiguously assign a compound class to the metabolite. The data is presented as a two-dimensional pie chart generated in excel.

Zebrafish larval skin is a bi-layer epithelium, in which the basal layer contains the stem cell compartment that gives rise to all strata of adult skin cells (Lee et al., 2014). A krtt1c19e promoter fragment has been shown to drive specific gene expression in this basal compartment (Lee et al., 2014). To achieve precise temporal control of oncogene expression, we used a tamoxifen inducible Gal4/UAS system (Figure 1A; Ramezani et al., 2015). This system relies on the combinatorial effect of two separate elements: the zebrafish optimized version of the transcription factor Gal4, KalTA4, fused to the modified ligand-binding domain of human estrogen receptor α (ER^T2) expressed under the control of the krtt1c19e promoter, and the Upstream Activating Sequence (UAS) controlling the expression of eGFP fused to either HRAS^G12V or CAAX control (van den Berg et al., 2019). We induced eGFP-HRAS^G12V expression in basal skin cells from 52 hpf (Figure 1A) to drive pre-neoplastic cell (PNC) development. The membrane eGFP labeling of eGFP-HRAS^G12V expressing PNCs allowed us to visualize morphological changes of PNCs compared with control cells (Figures 2A,B,D). At 12 hpi most of PNCs maintain their typical keratinocyte polygonal shape or slightly rounded but around 4.6% of PNCs appear to lose their polygonal shape and became more stretched out, by 24 hpi the morphology of PNCs diversifies further with 24% acquiring a more elongated or irregular shape or developing filopodia-like structures compared with polygonal shaped control eGFP-CAAX expressing cells (Figures 2B,D and Supplementary Figure 4). EdU incorporation experiments indicate that PNCs acquire an enhanced proliferation capacity at 24 hpi (22.75 ± 4.95%) compared with CAAX controls (0.55 ± 0.55%) and 12hpi PNC (12.41 ± 3.46) (Figures 2B,C and Supplementary Data File 1). The cell morphology changes and enhanced proliferation rate of PNCs is maintained at 36 hpi (Figures 2C,D). Since enhanced proliferation and clonal expansion are key features of pre-neoplastic development, and the initial morphological changes that might indicate the beginning of cellular state change started at 24 hpi, we decided to use larvae at 24 h post RAS oncogene induction for our untargeted metabolomics study to capture the early metabolic alterations in tumor initiation (Figure 1).

In order to survey metabolic perturbations in the HRAS^G12V induced skin pre-neoplastic lesions, we developed a whole-mount metabolite extraction protocol (Figures 1B,C) using an ice-cold homogenous mixture of chloroform/methanol/water (1:3:1) as extraction solvent. A group of intact larvae were immersed in this solvent mixture and incubated for 40 min allowing unbiased extraction of both polar and apolar metabolites from the larval skin tissue. Pigment extracted from skin could be seen dissolved in the solvent resulting in a yellowish solution and leaving intact larvae with a pale color (Supplementary Figure 3). The metabolic extract was subjected to chromatographic separation using zwitterionic-phase hydrophilic interaction chromatography (ZIC-pHILIC) coupled to MS-IMS detection in positive and negative ionization modes (Figure 1D). By using retention time, m/z and ion mobility drift time and CCS values, a total of 7,094 molecular features were obtained in positive and negative ionization modes combined. Multivariate statistical analysis using partial least squares-discriminant analysis (PLS-DA) on all the molecular features, separated the pre-neoplastic group from the wild-type sample group (Supplementary Figure 1A). Using the McLean CCS Compendium library consisting of 1,446 metabolites, 1,376 (19.4%) molecular features corresponding to 260 metabolites were annotated in the dataset across the two sample groups, preneoplastic vs. healthy larval skin (Figure 1E). Of these, 170 (65.4%) metabolites were unambiguously classified into 13 compound classes, with an overrepresentation (64 of 260 metabolites, 24.6%) of phospholipids, including the skin lipids phosphatidylcholine (PC), phosphatidylethanolamine (PE), phosphatidylserine (PS), and sphingomyelin (SM) (Figure 3A and Supplementary Data File 2). Additionally, the eumelanin precursor 5,6-dihydroxyindole-2-carboxylic acid was also detected (using accurate mass search) in the extracts of both sample groups at similar level (Supplementary Figure 1B and Supplementary Data File 2), demonstrating efficient and consistent metabolite extraction from the larval skin tissue across samples. To verify that our protocol is restricted to extracting metabolites from skin tissue only, we specifically queried the metabolomics dataset for lactate (using accurate mass search), which is one of the most enriched metabolites in muscle tissue just underneath the skin. Interestingly, we could not detect lactate in our samples, suggesting that our extraction was likely to be limited to skin tissue.

Further analysis of the 260 annotated metabolites showed differences in their relative abundance across the sample groups and the expected biological variation across the dataset (Figures 1F, 3B), emphasizing the need for multiple biological replicates to capture the true variation in the zebrafish larval skin model. A pathway enrichment analysis of the 260 annotated metabolites identified the representation of five statistically significant (p < 0.05) pathways corresponding to purine, amino acid and aminoacyl-tRNA metabolism in the dataset (Figure 3C).

To identify the metabolic pathways altered upon HRAS^G12V induction, the 260 annotated metabolites in the two sample groups were carried forward for multivariate statistical analysis using the MetaboAnalyst 5.0 platform. The PLS-DA analysis for WT and PNC samples clearly separated the two groups with principal components 1 and 2 explaining 9.2 and 12.8% of the variation, respectively (Figure 4A). A pathway enrichment analysis identified six metabolic pathways, including pyrimidine, purine and amino acid metabolism to be significantly (p-value < 0.05) altered across the sample groups (Figure 4B). Given the higher number of metabolites identified within the pathways for pyrimidine (6/41), arginine and proline (6/38) and purine (14/66) metabolism (Figure 4B, table), we focused on analyzing the possible role of these metabolic networks in preneoplastic cell development.

Nine metabolites involved in pyrimidine metabolism were annotated in our dataset (Figure 5A). Significant increase in the levels of 5-methylcytidine (p = 0.0081), uridine (p = 0.0285), and its isomer pseudouridine (p = 0.0019) were observed, suggesting a potential role for these metabolites in the development of preneoplastic lesions. Uridine is a key substrate for generating uridine monophosphate (UMP, not annotated in our dataset) by uridine-cytidine kinase in the salvage pathway (Supplementary Figure 2; Connolly and Duley, 1999). UMP is the key entry molecule for generating the full spectrum of pyrimidine nucleotides and nucleic acids (Supplementary Figure 2; Connolly and Duley, 1999; Kanehisa and Goto, 2000). Although UMP is not annotated in our dataset, we detected uridine diphosphate (UDP, p = 0.0826) and uridine triphosphate (UTP, = 0.8281), both of which were present at similar levels in the PNC and WT groups (Figure 5A and Supplementary Figure 2). This suggests the ability of PNC to maintain a pyrimidine nucleotide pool under high proliferative demand through a salvage pathway.

The conversion of uridine to 3-ureidopropionate is an important step for pyrimidine degradation (Kanehisa and Goto, 2000). Our data shows similar levels of 3-ureidopropionate (p = 0.9906) in PNC and WT groups (Figure 5A and Supplementary Figure 2), suggesting similar levels of pyrimidine degradation. Additionally, cytidine (p = 0.1783), 5-methycytosine (p = 0.3454), and thymidine monophosphate (TMP, p = 0.6575) were also annotated in the dataset, and no significant changes in their relative abundance was observed between the sample groups (Figure 5A and Supplementary Figure 2).

Unlike the detectable higher abundance for some pyrimidine metabolites, a reduced relative abundance of guanine (p = 0.0241) and adenosine triphosphate (ATP, p = 0.0018) was detected in the PNC samples compared to the control group (Figure 5B). Fourteen additional purine metabolites were annotated in the dataset with no significant (p > 0.05) differences in their relative abundance between the sample groups, although a slight trend of increase could be visualized for three metabolites (ribose-5-phosphate, adenosine monophosphate and inosine) in the PNC samples (Figure 5B and Supplementary Figure 2).

In addition to ATP, the energy charge metabolites adenosine monophosphate (AMP, p = 0.2669) and adenosine diphosphate (ADP, p = 0.6585) were also annotated in our dataset. Since reduced ATP to ADP ratio has been reported to favor enhanced glycolysis in proliferating cells (Maldonado and Lemasters, 2014), we determined this ratio in our study. We calculated the ATP to ADP ratio using the raw peak intensity values, and indeed found a trend of reduced ATP/ADP ratio (p = 0.084) in the PNC group compared to WT group (Supplementary Figure 1C), which likely correlates with enhanced proliferation in preneoplastic lesions.

In arginine and proline metabolism, significantly lower relative abundance of hydroxyproline (p = 0.0165) and proline (p = 0.0221) was observed in PNC samples (Figure 6A). In addition, four other metabolites were detected in our dataset, of which arginine (p = 0.1932) and guanidinobutanoate (p = 0.2886) show a trend of lower abundance (Figure 6A), while creatine and phosphocreatine showed no differences in their relative abundance, compared to the WT group. Given these differences and because other amino acid metabolic pathways were among the pathways identified to be significantly altered (Figure 4B), we analyzed the relative abundance of all amino acids in our dataset. Of the 20 amino acids, 15 were detected in our untargeted analysis (Supplementary Data File 2). Ten of those amino acids were relatively lower in PNC group (Figure 6B). Using PLS-DA analysis and a variable importance in projection (VIP) cut-off score of 1, six amino acids were identified to be changed of which alanine (p = 0.0147), hydroxyproline (0.0165), proline (0.0221), and glutamic acid (0.0442) were significantly lower in the PNC group (Figures 6A,C). Additionally, leucine was found to have a trend of lower abundance (p = 0.0584) and aspartic acid higher (0.0820) in the PNC group (Figure 6C), suggesting alterations in amino acid biosynthetic pathways upon HRAS^G12V induced preneoplastic development.

Many rate-limiting metabolic enzymes are known to be deregulated in cancer (Furuta et al., 2010). The changes detected in our untargeted metabolomics analysis would suggest altered activities of some of these metabolic enzymes. However, it is not clear whether these metabolic enzyme changes occur at the level of gene expression upon oncogene induction. Therefore, we determined the expression level of a handful of candidate genes implicated in cancer and found several of these to be upregulated in PNCs compared to control cells (Figure 7). Relevant to nucleotide metabolism pathways, are the genes rrm2, itpa and pnp5b, all of which were significantly overexpressed in the PNCs (Figures 7A–C). Rrm2 encodes the ribonucleoside-diphosphate reductase subunit M2, which catalyzes the rate-limiting biosynthesis of deoxyribonucleotides from the corresponding ribonucleotides to support DNA synthesis and cell proliferation (Kittler et al., 2004; Furuta et al., 2010). Enhanced expression of rrm2 in PNCs would support their enhanced proliferation. Inosine triphosphatase (itpa) is involved in recycling purines by converting (deoxy)inosine triphosphate [(d)ITP] to (d)IMP or xanthosine 5′-triphosphate (XTP) to XMP, which could then be recycled to contribute to DNA and RNA synthesis (Simone et al., 2013). Increased itpa expression in PNCs would support altered nucleotide metabolism, which we observed in our metabolomics data, although we did not detect changes in substrates or products of itpa, suggesting altered flux through the associated pathways rather than the accumulation of metabolites directly linked to the enzyme. Pnp5b encodes zebrafish purine nucleoside phosphorylase 5b, that is homologous to mammalian purine nucleoside phosphorylase, which catalyzes the reversible phosphorylation of purine nucleosides (Williams et al., 1984). Upregulation of pnp5b could contribute to altered levels of purine metabolites, such as guanine, hypoxanthine, guanosine and inosine, detected in our metabolomics analysis. Interestingly, we also observed increased expression of sat1a.2 (Figure 7D), which encodes spermidine/spermine N(1)-acetyltransferase, a rate –limiting enzyme that catalyzes the N1-acetylation of spermindine and spermine (Allmeroth et al., 2021), suggesting that altered polyamine removal might occur in PNCs. However, these metabolites were not annotated in the metabolomics dataset.