Dietary palmitic acid promotes a prometastatic memory via Schwann cells

Fatty acid uptake and altered metabolism constitute hallmarks of metastasis1,2, yet evidence of the underlying biology, as well as whether all dietary fatty acids are prometastatic, is lacking. Here we show that dietary palmitic acid (PA), but not oleic acid or linoleic acid, promotes metastasis in oral carcinomas and melanoma in mice. Tumours from mice that were fed a short-term palm-oil-rich diet (PA), or tumour cells that were briefly exposed to PA in vitro, remained highly metastatic even after being serially transplanted (without further exposure to high levels of PA). This PA-induced prometastatic memory requires the fatty acid transporter CD36 and is associated with the stable deposition of histone H3 lysine 4 trimethylation by the methyltransferase Set1A (as part of the COMPASS complex (Set1A/COMPASS)). Bulk, single-cell and positional RNA-sequencing analyses indicate that genes with this prometastatic memory predominantly relate to a neural signature that stimulates intratumoural Schwann cells and innervation, two parameters that are strongly correlated with metastasis but are aetiologically poorly understood3,4. Mechanistically, tumour-associated Schwann cells secrete a specialized proregenerative extracellular matrix, the ablation of which inhibits metastasis initiation. Both the PA-induced memory of this proneural signature and its long-term boost in metastasis require the transcription factor EGR2 and the glial-cell-stimulating peptide galanin. In summary, we provide evidence that a dietary metabolite induces stable transcriptional and chromatin changes that lead to a long-term stimulation of metastasis, and that this is related to a proregenerative state of tumour-activated Schwann cells. Palmitic acid induces stable transcriptional and chromatin changes that lead to long-term stimulation of metastasis in orthotopic models of cancer through the secretion by tumour-associated Schwann cells of a specialized proregenerative extracellular matrix, the ablation of which inhibits metastasis initiation.


Article
CD36 bright cells, and even CD36 dim cells (that is, with low levels of CD36, previously described as non-metastatic 6  To test whether a PA-induced prometastatic memory could also come from a palm-oil-rich diet, we orthotopically inoculated OSCC cells into mice; once primary tumours appeared, mice were fed short term (10 d) with a diet rich in palm oil or olive oil or a standard diet; the mice were then fed a standard diet until being euthanized at the end point. We next collected primary tumours (all of which contained a similar number of CD36 bright cells), purified OSCC cells and serially transplanted tumour cells into secondary mice that were fed a normal diet ( Fig. 1c and Extended Data Fig. 2a, b). Corroborating our in vitro results for PA and OA, tumour cells from palm-oil-fed primary recipient mice were much more metastatic, and those from olive-oil-fed mice were less metastatic, in secondary recipient mice (Fig. 1d, Extended Data Fig. 2c, d and Supplementary Table 1j). The PA-induced stable boost in metastasis also occurred after serial orthotopic inoculation of melanoma cells under the same dietary conditions (Extended Data Fig. 2h, i). Removing CD36 from the primary tumour (by constitutive CD36 knockdown (KD) or treatment with CD36-neutralizing antibodies) led to a loss of the prometastatic phenotype but did not affect the primary tumour competency (Extended Data Figs. 2e-i and 3a and Supplementary Table 1j-n). CD36 depletion using an inducible short-hairpin RNA (shRNA) in the secondary recipient mice also inhibited the prometastatic memory of PA in the tumour cells (Extended  Data Fig. 3b, c and Supplementary Table 1o, p).
To investigate whether the epigenomes of tumour cells are affected by PA regarding this metastatic memory, we treated OSCC cells in vitro for 4 d with PA or OA, followed by 14 d without PA. Next, using chromatin immunoprecipitation followed by sequencing (ChIP-seq), we studied the genome-wide changes in histone marks that are indicative of promoters (histone H3K4me3), enhancers (histone H3K4me1), active transcription (H3K27ac), polycomb-mediated repression (histone H3K27me3) and non-polycomb repression (H3K9me3) (Fig. 2a-c, Extended Data Figs. 3d-h, 4a and 5a and Supplementary Tables 2a, b, 3a-d and 4a, k-q). We observed stable changes in H3K4me3 in the promoters of specific genes after either PA or OA treatment, but a significantly higher number of genes were affected by PA exposure compared with OA (Extended Data Fig. 4b, e and Supplementary Table 5c, f, g, o-q). A stably altered H3K4me3 profile was also evident in vivo by ChIP-seq analysis of secondary tumour transplants that had retained the prometastatic memory of a palm-oil-rich diet described in Fig. 1 (Extended Data Fig. 4c-d and Supplementary Table 2e). Importantly, the overlap of H3K4me3-marked genes after exposure to PA in vitro or palm oil in vivo was statistically significant, underscoring the physiological relevance of our in vitro studies (Extended Data Fig. 4e, f).
The distribution of H3K4me1, H3K27me3 and H3K27ac was also altered after treatment with PA for 4 d (Extended Data Figs. 4a and 5a  and Supplementary Tables 3a and 4d, o, q), but the majority of these changes had faded 14 d after PA withdrawal (Extended Data Figs. 3d and 5a and Supplementary Tables 3b-e, 4d, p and 5h). Indeed, no stable changes in H3K4me1 and H3K27ac were observed in the enhancers of the genes with long-term differences in H3K4me3 induced by PA, as determined by ChIP-seq analysis, and nascent transcription at enhancers through precision run-on sequencing (PRO-seq) analysis (Extended Data Fig. 3g, h and Supplementary Table 3f-h). Furthermore, promoters with the PA-induced stable changes in H3K4me3 were not bivalent, as they did not show stable changes in H3K27me3 (Extended Data Fig. 3e and Supplementary Table 5d, e). Notably, the H3K9me3 levels were significantly reduced in several genes in a stable manner after PA treatment (even after PA withdrawal); however, there was less than a 2% overlap between these genes with reduced H3K9me3 and those that gained H3K4me3 in their promoters (Extended Data Fig. 3i, j and Supplementary Table 4a-c). Taken together, these results indicate that, although different chromatin marks respond to PA stimulation, these changes are stably maintained mainly by (1) genes with newly marked H3K4me3 promoters and (2) genes with specific regions that lose the repressive H3K9me3 mark.
To distinguish whether PA treatment can induce the CD36 + prometastatic state de novo or whether PA only enhances the potential of CD36 + cells already present within a tumour, we analysed OSCC cells that were stimulated with PA for 4 d in culture using t-SNE-based cluster and pseudotime analyses of single-cell RNA-sequencing (scRNA-seq). Notably, almost all cells responded to PA and acquired the CD36 + signature of metastasis-initiating cells, including genes such as CD36 and IFIT3 (Extended Data Fig. 5b, c and Supplementary Table 4s, t). Furthermore, cells slowed down their cell cycle during PA treatment but were proliferating normally 14 d after PA withdrawal, with no associated apoptotic death at any point, ruling out the selection or out-competition of CD36 + cells over any other cell state (Extended Data Fig. 5d and Supplementary  Table 4e). These results indicate that PA induces, rather than selects, a stable prometastatic memory in OSCC cells.
Gene Ontology (GO) analyses of the genes with the PA-induced stable changes in H3K4me3 revealed many functions related to neurogenesis Re-injection of fat-diet-exposed OSCC cells   (Met.) in mice injected as shown in a, expressed as percentages. n = 15 mice per group. Statistical analysis was performed using two-tailed Fisher's exact tests; for lung metastases, ***P = 0.0007, *P = 0.04. SA, stearic acid. PA/OA denotes 4 d treatment with PA followed by 4 d treatment with OA and then 14 d in culture without treatment. c, Schematic of the in vivo diet experiments in which OSCCs were exposed to a high-fat diet in primary recipient mice and injected into secondary recipients fed a control diet. d, The frequency of developed lymph node (LN) metastasis and lung metastasis of mice injected as shown in c, expressed as percentages. n= 20 mice per group. Statistical analysis was performed using two-tailed Fisher's exact tests; for lymph node metastasis, *P = 0.003, **P = 0.02; for lung metastasis, *P = 0.01, ***P = 0.0004. b and d contain data from two biological replicates. and neural remodelling (Supplementary Tables 2c, d, f and 5a, b, i, j, l). By contrast, the levels of H3K4me3 were downregulated in genes pertaining to the same neural categories after OA stimulation (Supplementary Table 5m, n). We were particularly intrigued by these categories, as perineural invasion and tumour innervation are among the strongest predictors of poor prognosis in oral, prostate, breast and pancreatic carcinomas, among other cancer types 3,4,[12][13][14][15][16] . However, the mechanisms that underlie the contribution of perineural invasion and tumour innervation to tumour aggressiveness remain mostly unknown.
There was a significant correlation between the stable changes in H3K4me3 and differentially expressed genes (DEGs) in OSCC cells stimulated with PA in vitro (Extended Data Fig. 5e and Supplementary Tables 4f, g and 5r, s). Many of these DEGs are related to neural categories, among other biological functions (Supplementary  Tables 4h-j and 5t). Importantly, the same neural categories appeared (in a CD36-dependent manner) in the transcriptome of CD36 bright cells purified from oral tumours and melanomas harbouring the in vivo prometastatic memory of the palm-oil-diet-derived but not in olive-oil-diet-derived secondary tumours (Fig. 3a, Extended Data  Tables 6 and 7). These proneural transcriptional changes were functionally relevant, as chemical sympathectomy through intraperitoneal administration of the noradrenaline-selective neurotoxic drug 6-hydroxydopamine (6-OHDA) strongly reduced the prometastatic memory of PA without affecting primary tumour competency (Extended Data Fig. 7a-c and Supplementary Table 8a-d).
Transcription-factor-binding analysis of the promoters of the in vivo palm-oil-diet-associated, neural-related DEGs of OSCC cells showed a strong enrichment for EGR2 binding sites (Extended Data Fig. 7d and Supplementary Table 9). Furthermore, a search for the in vivo neural signature induced by PA that was common to all of the tumours that we had tested resulted in a list of genes including the glial-inducing peptide galanin (GAL), which was one of the most upregulated genes 17,18 (Extended Data Fig. 7e and Supplementary Table 9). The genomic localization of EGR2 responded to PA stimulation of OSCC cells, and a significant proportion of the genes that gained EGR2 binding, including galanin, were neural-related in a CD36-dependent manner (Fig. 3b, Supplementary Table 10b-h, k, l, n), which was enriched for neural-related genes (Extended Data Fig. 8b and Supplementary Table 5a). Importantly, depletion of either EGR2 or galanin prevented CD36 bright tumour cells from acquiring a stable neural signature (Extended Data Fig. 8c, d and Supplementary Table 9) and strongly reduced the in vivo prometastatic memory of OSCC established by the palm-oil-rich diet (Fig. 3d, Extended Data Fig. 8e and Supplementary  Table 9). This antimetastatic effect was also seen after systemic inhibition of galanin signalling through the intraperitoneal administration of the pan-galanin receptor inhibitor galantide (M15) (Extended Data Fig. 9a-c and Supplementary Table 11a-c). Expression of EGR2 and GAL was significantly correlated with poor prognosis in patients with head and neck squamous cell carcinoma (SCC) as well as other types of solid tumours (Supplementary Table 12). H3K4 trimethylation is deposited by methyltransferases of the MLL/ COMPASS family, including MLL1, MLL2, Set1A and Set1B 19,20 . In OSCC cells 14 d after PA withdrawal, no changes were observed in either the protein expression of these four COMPASS family members (Extended Data Fig. 10a, b) or in the localization of MLL1 and MLL2 on chromatin (as determined by ChIP-seq) (Extended Data Fig. 10c- e and Supplementary Table 13a, d, e). Furthermore, inhibition of MLL1 did not affect the long-term PA-induced metastatic potential of OSCC cells after PA stimulation (Supplementary Table 13b, c). By contrast, OSCC cells depleted of Set1A did not acquire the stable PA-induced H3K4me3 neural signature (Extended Data Fig. 10f, g and Supplementary Table 14a, b, f, g), and their ability to retain a PA-induced prometastatic memory in vitro and in vivo was significantly reduced (Extended Data Fig. 10h). Although the size of the primary tumours was diminished in Set1A-depleted cells, their penetrance was unchanged (Supplementary Table 14c-e).
We next studied whether the predominant neural signature of tumours retaining the prometastatic memory of PA influenced the composition and function of the tumour stroma, specifically relating to its neural component. Notably, the bulk transcriptome of the stroma from these OSCC tumours was distinct from that of tumours derived from all other dietary conditions in a CD36-dependent manner (Extended Data Fig. 11a Table 15). Node pathway interactome predictions further highlighted neuron projection and Schwann cell development as functions significantly enhanced in the stroma of palm-oil-derived tumours (Extended Data Fig. 12a and Supplementary Table 15). We next performed 10x scRNA-seq to specifically characterize the cellular composition of the bulk tumour stroma that we had purified. Uniform manifold approximation and projection (UMAP) analysis identified five clusters annotated to erythrocytes, vascular and lymphatic endothelium, smooth muscle cells and tumour-associated Schwann cells (Fig. 4a, Extended Data Fig. 12b, c and Supplementary Table 16). Within these five cell types, tumour-associated Schwann cells and lymphatic endothelial cells displayed the highest proportional changes in gene expression between both conditions (that is, PA-induced memory compared with control tumours) in a CD36-dependent manner (Fig. 4a). However, the DEGs relating to ECM organization, neurogenesis and gliogenesis/Schwann cells, which we identified as the leading GO terms in the analysis of the bulk transcriptome data, corresponded primarily to the tumour-associated Schwann cell population (Fig. 4b, Extended Data Fig. 12d and Supplementary Table 16).
We were particularly interested in the expression changes in ECM components that are specifically associated with the tumour-associated Schwann cell cluster. Several of these ECM proteins are components of perineuronal nets, a specialized matrix that is secreted by glial cells during nerve injury and is essential for nerve regeneration 21-23  Table 16). Using the TCGA database, we observed a significantly positive correlation between their expression and the presence of the PA-induced CD36 + metastatic signature in human head and neck SCCs and melanomas (Extended Data Fig. 13d). Immunofluorescence analyses further showed their increased expression in the stroma of the tumours with PA-induced prometastatic memory in very close proximity to tumour-associated Schwann cells compared with their control diet counterparts (Extended Data Fig. 13e, f). Spatial transcriptomics enabled us to clearly define the bulk and invasive front of the tumours surrounded by the healthy mouse tissue of the tongue (Extended Data Fig. 14a, b). Lesions that retained the PA prometastatic memory showed a strong increase in the number of Schwann cells infiltrating the tumour compared with the control tumours in a CD36-dependent manner (Extended Data Fig. 14b, c). This was confirmed by immunofluorescence analysis of the Schwann-cell-specific marker S100 (Extended Data Fig. 15a-f). Positional RNA-seq analysis also confirmed high levels of ECM components related to tumour-associated Schwann cells (Extended Data Fig. 14b).
We next expressed a bacterial chondroitinase ABC (chABC; modified to be expressed in eukaryotic cells) in OSCC cells (Extended Data Fig. 16a), which specifically digests the Schwann-cell-derived ECM components, enabling us to directly test their potential role in metastasis 24 . Notably, chABC expression prevented the stable increase in metastatic competency induced by PA in OSCC cells, without affecting the penetrance of primary tumours ( Fig. 4c  of Schwann-cell-specific ECM components, but not of collagen 5A1 (which is not produced by Schwann cells), was confirmed by immunofluorescence staining (Extended Data Fig. 16b, c). Metastatic progression of tumours requires driver mutations, but is also associated with stable (that is, epigenetic) alterations in chromatin. For example, pancreatic adenocarcinoma cells modify large euchromatic regions that are decorated by H3K9me3 downstream of the pentose phosphate pathway to acquire metastatic competency 25 . Moreover, metastasizing melanomas show global increases in H3K27me3 and H3K9me3, and depend on the histone methyltransferases SETDB1 and EZH2 to efficiently metastasize 26 . Our results indicate that dietary PA not only stimulates metastasis, but that it also does so in a long-term stable manner through a transcriptional state that stimulates intratumoural Schwann cells. This prometastatic memory is retained only by tumour cells and not by the tumour stroma (Supplementary Table 18).
Schwann cells are essential for the protection of peripheral nerves as well as for their recovery after damage 27,28 . The presence of intratumoural Schwann cells is associated with metastasis of numerous cancers, although the factors responsible for their stimulation and the mechanism behind how they promote metastasis are ill-defined 3,4,12-16 .
Our data indicate that intratumoural Schwann cells activated by metastatic cells downstream of dietary PA adopt a proregenerative state related to the secretion of a specialized ECM, which occurs at least in part through communication with the neuropeptide galanin that is secreted by CD36 + metastatic cells. Research has shown that the presence of the CD36 + signature is correlated with poor survival in several tumour types 29,30 . We now show that there is a significant positive correlation between the PA-induced CD36 + metastatic signature and the expression of perineuronal net components in many different tumour types (Supplementary Table 19).
Our data suggest that this proneural/prometastatic state is established in part through the transcription factor EGR2 and stabilized by the methyltransferase Set1A/COMPASS downstream of CD36. Notably, CD36 is required for fat tasting and sensing 31 , and Set1A deposits H3K4me3 at neural genes to epigenetically establish the neural lineage during development 32,33 . Thus, one hypothesis is that metastatic cells aberrantly hijack the natural fat-sensing function of CD36 that is normally mediated by the interactions between specialized cells, such as oral keratinocytes and intestinal epithelial cells, and their associated nerves and glial cells, thereby establishing specific epigenetic states involved in fat metabolism and signalling 31,34 . Extensive future work will be necessary to test this possibility.
Different tumour types seem to display specific FA preferences related to metastasis. For example, OA inhibited the metastatic spread of oral carcinoma and melanoma in our orthotopic models, yet it stimulates metastasis of cervical and gastric carcinomas 35 . Moreover, metastasizing melanoma cells that invade through the lymphatic system acquire a high metastatic competency by being exposed to high levels of OA present in the lymph nodes 36 . Different tumours also fuel their growth and progression by converting palmitate into sapienate 37 . However, irrespective of the FA, inhibition of FA metabolism has potential for treating numerous tumours with an unmet clinical need 30 . In summary, our findings not only underscore the long-term health risks associated with a diet rich in PA regarding metastatic progression, but also provide mechanistic insights to identify new epigenetic-and neural/glial-related therapeutic strategies to attenuate and prevent distant dissemination. or LV-chABC-derived cells from primary recipient mice fed a normal diet (CT) or a palm-oil-rich diet. n = 12 (WT VDH-15) and n = 10 (chABC VDH-15) mice per group. Statistical analysis was performed using two-tailed Fisher's exact tests; for lymph node metastasis frequency, *P = 0.03, **P = 0.008; for lung metastasis, **P = 0.002, **P = 0.001. The frequency is expressed as a percentage. c contains data from one biological replicate.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-021-04075-0. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Animal studies
Groups of 8-10-week-old NSG (NOD.Cg-Prkdc scid Il2rg tm1Wjl /SzJ) mice were purchased from Charles River and crossed in-house. All of the mice (of an approximate 50:50 mix of male and female) were housed under a regimen of 12 h-12 h light-dark cycles and specific pathogen free conditions, and all of the procedures were evaluated and approved by the Ethical Committee for Animal Experimentation (CEEA) from the Government of Catalunya. SSC intratongue injection was performed as previously described 38,39 . In brief, mice were anaesthetized by intraperitoneal injection with a mixture of ketamine (50 mg kg −1 ) and medetomidine (0.5 mg kg −1 ), and SCC cells were resuspended in 30 µl PBS and injected into each mouse tongue using a BD ultrafine 6-mm needle. For melanoma orthotopic cell transplant, mice were anaesthetized, and a 0.8-cm incision was made on the dorsal back skin to implant a silicon chamber (ref. CSAT60BMOL1018, Merefsa, S.L.) as previously described 40 . A mixture with 100,000 melanoma cells and 100,000 primary human keratinocytes was introduced into the silicon chamber using a p200 pipet. Mice were housed individually and anaesthetized after 8 d, and the chamber was removed to bring the grafted area into direct contact with air. Mice were monitored for the luciferase bioluminescence signal immediately after injection (T 0 ) and then once weekly with a Xenogen IVIS Imaging System-100 (Caliper Life Sciences). In brief, mice were injected by retro-orbital injection with 50 µl of d-luciferin (Promega) diluted in 1×PBS at 5 mg ml −1 . Continuous administration of isofluorane gas was provided to ensure that the mice were anaesthetized during imaging. Data were quantified using the Living Image software (v.4.4; Caliper Life Sciences). Quantifications were calculated with unsaturated pixels. Colour scale minimum and maximum values are shown in pictures. High-fat diet experiments entailed feeding mice for 10 d with 42% Kcal fat-modified western diet supplemented with either palm oil (TD 150067, Envigo) or olive oil (TD 09820, Envigo). Normal chow diet was used for the control groups.
For histological analysis, mouse tissue was collected, fixed with 4% paraformaldehyde for 2 h at room temperature, dehydrated and then embedded in paraffin.
To induce expression of inducible shCD36, doxycycline hydrochloride (D3447 Sigma-Aldrich) was diluted at dose of 2 mg ml −1 in drinking water with 5% sucrose and continuously administered to mice in secondary recipients.
For blocking experiments with M15 (galantide/galanin receptor antagonist, G1278, Sigma-Aldrich), secondary recipient NSG mice were injected i.p. with either vehicle alone or M15 prepared in dimethyl sulfoxide, injected at dose of 40 nmol kg −1 every 4 d after OSCC orthotopic injection 42 , starting on the day of the OSCC injection and finishing 1 d before the final point of the experiment.
For all experimental functional assays, male and female (50:50 mix) mice were used. No randomization method was used. Female mice were used as tumour recipients for the molecular characterization of the tumour cells. For all of the experiments, mice were euthanized at the approved humane end point and before the tumours reached the maximum tumour size of 1 cm in diameter.

Clinical material
Biological samples and data from patients included in this study were provided by the Vall d'Hebron University Hospital Biobank (PT17/0015/0047) integrated in the Spanish National Biobanks Network, and they were processed following standard operating procedures with the appropriate approval of the Ethical and Scientific Committees of the hospital according to Spanish ethical regulations. The study followed the guidelines of the Declaration of Helsinki, and patient identity and pathological specimens remained anonymous (in the study and otherwise).
Cell culture studies Cells were cultured in a humidified incubator at 37 °C with 5% CO 2 . The human tongue SCC cell line SCC-25 (ATCC, CRL-168) and patient-derived cells (VDH-15) were grown in keratinocyte serum-free medium (KSFM; Gibco) supplemented with 5 µg ml −1 penicillinstreptomycin, 0.025 mg ml −1 bovine pituitary extract and 0.2 µg ml −1 hEGF. The 501mel human cell line (provided by C. Wellbrock, Manchester Cancer Research Centre, The University of Manchester) was grown in DMEM supplemented with 5 µg ml −1 penicillin-streptomycin and 10% FBS (Gibco). Primary human keratinocytes were grown as previously described 44 . PhoenixA and HEK293T cells grown in DMEM/10% FBS were used for retrovirus and lentivirus production, respectively, after transfection using the CaCl 2 method. For selection, 2.5 µg ml −1 puromycin or 0.3 mg ml −1 G418 was added to the medium. No method of cell line authentication was performed. All cell lines tested negative for mycoplasma contamination.
For PA, OA and LA in vitro OSCC treatment, sodium palmitate (P9767, Sigma-Aldrich), sodium oleate (O7501, Sigma-Aldrich), sodium stearate (S3381, Sigma-Aldrich) and linoleic acid sodium salt (L8134, Sigma-Aldrich), respectively, were prepared as a 2.5 mM stock solution by dissolving in 1 ml solution of 0.1 M NaOH and warming at 80 °C until clear. The FA solution was complexed with fatty-acid-free bovine serum albumin (BSA; A7030, Sigma-Aldrich) in a molar ratio FA/BSA of 5:1. In brief, 0.325 g BSA was dissolved in 8 ml of 0.9% NaCl, and the mixture was warmed to 45 °C. The clear solution of palmitate, oleate, stearate and linoleate was added drop-by-drop by pipet with agitation. The final stock solution was filtered at 0.45 µM, aliquoted and stored at −20 °C. SCC cells growing in serum-free medium were treated in vitro with 300 µM or 50 µM palmitate, 50 µM stearate, 50 µM oleate, 50 µM LA in KSFM medium for 96 h (4 d). For intratongue injection of SCC cells and backskin melanoma transplant, cultured cells were flask-detached with 0.25% trypsin-EDTA (Gibco), diluted in PBS/Trypan Blue, counted in a Neubauer chamber and resuspended in 1× PBS. For SCC-25 and VDH-15 cells, 100,000 cells per mouse were injected into the mouse primary recipients. For SCC-25 and VDH-15 injection into secondary recipients, 10,000 SCC DAPI − GFP + purified sorted cells were injected per mouse. Cells were resuspended in cold KSFM/30% Matrigel (BD Matrigel, growth factor reduced, 366231) before being injected into secondary recipients. For melanoma transplant, a mixture of cells containing 100,000 melanoma cells plus primary human keratinocytes (1:1) was introduced into the silicon chamber of primary mice recipients. For secondary recipients, 10,000 melanoma cells plus primary human keratinocytes (1:1) were transplanted.
For doxycycline in vitro treatment, SCC cells were incubated at 37 °C with 5% CO 2 with 2 µg ml −1 doxycycline hydrochloride (D3447 Sigma-Aldrich), diluted in KSFM medium and were collected 24 h to 48 h later for flow cytometry analysis.

Disaggregation of tumours from xenografted mice
Tumours were isolated from mouse tongues or melanoma grafts, and connective tissue was removed to the greatest extent possible. Tissue was chopped in 0.5% trypsin 1-300 (MP Biomedical) in KSFM medium (Gibco) in a McIlwain Tissue Chopper. After complete homogenization, the samples were incubated at 37 °C for 90 min with shaking. Homogenates were filtered sequentially in 100-µm, 70-µm and 40-µm BD strainers and then centrifuged at 1,000 r.p.m. for 10 min at 4 °C. The supernatant was discarded, and each pellet was resuspended in 1× PBS/4% calcium-chelated FBS 44 for antibody staining and subsequent fluorescence-activated cell sorting (FACS) analysis.  44 and centrifuged for the second incubation with APC Streptavidin (405207, BioLegend, 1:50 dilution). Unstained, single-colour and fluorescence minus one controls were used in each case. The samples were analysed using the BD FACS Aria Fusion system. For FACS sorting, cells were selected on the basis of their forward and side scatter, excluding cellular debris. Doublets and dead cells were eliminated by DAPI. GFP-positive cells were gated excluding lineage-BV-positive cells. This population was selected for further CD36 human expression analysis. Tumour-associated stromal cells were gated excluding endothelial (CD31 + ), macrophages (CD45 + ) and human (GFP + ) cells. This population was selected for further mouse stromal gene expression analysis.

FACS analysis
10x scRNA-seq Sample processing. Xenografted tumours were processed and FACS-sorted as previously described (see the 'Disaggregation of tumours from xenografted mice' section) to purify a tumour stroma-associated fraction enriched in glial and neural progenitors. Cells were collected in PBS + 0.5% BSA at 4 °C in LoBind tubes (Eppendorf) and processed immediately.
Sequencing. The cell concentration and viability of the single-cell suspension were verified by counting using the TC20 Automated Cell Counter (BioRad). Cells were partitioned into gel bead-in-emulsions (GEMs) using the Chromium Controller system (10x Genomics) with the aim of a target cell recovery of 5,000 cells. Single-cell gene expression (GEX) libraries were prepared using the Chromium Single Cell 3′ Library & Gel Bead Kit v2 (10x Genomics, 120237) following manufacturer's instructions. In brief, after GEM-reverse transcription clean-up, cDNA was amplified using 11 cycles. cDNA quality control and quantification were performed using the Agilent Bioanalyzer High Sensitivity chip (Agilent Technologies). Libraries were indexed by PCR using the PN-220103 Chromiumi7 Sample Index Plate (10x Genomics). The size distribution and concentration of 3′ gene expression libraries were verified using the Agilent Bioanalyzer High Sensitivity chip (Agilent Technologies). Sequencing was performed using the Illumina Hiseq4000 sequencer to obtain approximately 40,000 reads per cell.
Quality control and normalization. Quality control and normalization were performed independently for each dataset. The distributions of three quality control metrics were assessed both individually and jointly: library size (total unique molecular identifier (UMI)), library complexity (the number of detected genes) and mitochondrial expression. Note that analysing these metrics jointly ensured that cells with high mitochondrial expression were not metabolically active. Cells barcodes with an aberrantly low number of UMI and genes, or with an abnormally high mitochondrial percentage, were classified as poor quality. Doublets were classified as barcodes with an aberrantly large library size and complexity, and by using the Scrublet doublet detection 45 . Thresholds for the aforementioned parameters were set individually for each dataset. Genes that were detected in less than eight barcodes were rules out. Seurat (v.4.0.0) 46 was used to normalize and scale UMI counts, using the function SCTransform 47 .
Data integration. The four datasets were first merged without correcting for any technical artifacts, which revealed marked dataset-dependent clustering. These datasets were then integrated using Harmony 48

Downstream analysis.
To annotate the cell identity, cells were clustered using the FindClusters function. The resolution parameter was set to 0.05 for a high-level annotation. Cluster 5 was then further split into two subclusters on the basis of the cluster labels obtained at a resolution of 0.5. Marker genes were computed for each cell identity, and a Wilcoxon signed-rank test was performed to test for differential expression of each gene; note that this test is among the best performing for scRNA-seq differential-expression analysis 49 . Finally, a combination of differentially expressed and well-known marker genes was used to annotate each cluster to its specific cell type. Finally, the differential abundance analysis between the control and palm-oil diet in the different conditions was performed using a binomial logistic regression. The control diet group was set as the reference and obtained the odds ratio for each type.
Sequencing. Frozen tissue samples were tested for RNA quality with an RNA integrity number > 7.0 (RNA pico, Agilent), and a tissue optimization experiment (10x Genomics, Visium Spatial Tissue Optimization, Rev A) was performed with imaging of the fluorescence footprint on upright epifluorescence (Nikon Eclipse E800), which identified 10 min as an optimum permeabilization time.
The samples were then processed for full spatial transcriptomic (ST) experiment according to the manufacturer's instructions (10x Genomics, Visium Spatial, Rev D), cut in a precooled cryostat at 10 µM thickness onto four 6.5 mm × 6.5 mm capture areas with 5,000 oligo-barcoded spots. The slides then were fixed and stained by H&E with immediate imaging on NanoZoomer S60 (Hamamatsu). Tissues were permeabilized with a proprietary enzyme (10 min), and reverse-transcription and second-strand synthesis were performed on the slide, with cDNA quantification using quantitative PCR with reverse transcription (RT-qPCR) using the KAPA SYBR FAST-qPCR kit (KAPA Biosystems), and analysed on the LightCycler480 Real-Time PCR System (Roche). RT-qPCR results (C q value at 25% of peak fluorescence) revealed cDNA amplification. The Visium slide serial number was V19S23-040.
Libraries were constructed according to the manufacturer's instructions. ST libraries were then quantified using the High-Sensitivity DNA Assay 2100 (Agilent). Libraries were sequenced on the NovaSeq6000 platform (Illumina) using a 100 bp paired-end dual-indexed set-up. Around 200 million paired-end reads were obtained for each tissue section. Read 1 was sequenced with 28 cycles; i7 index, 10 cycles; i5 index, 10 cycles; and read 2, 82 cycles.
Primary processing. Reads were mapped against a mixed human (GRCh38)-mouse reference dataset downloaded from 10x Genomics (https://support.10xgenomics.com/single-cell-gene-expression/ software/downloads/latest) using 10x Genomics Space Ranger v.1.1.0 with the default parameters. Images were processed using Inkscape to divide them into individual capture areas, and then rotated to align the printed fiducial spot pattern.
Quality control and normalization. Quality control was performed by examining the number of UMIs and genes, and the mitochondrial and ribosomal percentage. For shCD36/palm-oil diet samples, spots that did not overlap with the tissue sections were removed. Expression matrix counts were normalized using Seurat's NormalizeData, which uses log1p to natural-log-transform the data and then multiplies the values by a scale factor of 10,000. A subset of the 2,000 most highly variable genes was identified using the FindVariableFeatures function using the default parameters; this subset was further scaled to mean 0 and variance 1 using the function ScaleData.
Downstream analysis. PCA was used to reduce the dimensionality using RunPCA, clustering and nonlinear dimensionality reduction were run using the functions FindNeighbors, FindClusters and RunUMAP, respectively, on the top 10 principal components. The number of principal components was determined by looking at the elbow plot of the variance percentage explained by each principal component.
The percentage of mouse versus human for each spot was set by determining the fraction of UMIs mapped to mouse genes over human genes at each capture location. Tissue region was stratified by clusters overlapping pure mouse, pure human and intermixing regions as being healthy, tumour and tumour front, respectively.
To map the single-cell reference dataset, some granularity was added to the basal cells by dividing it into five subclusters and removing erythrocytes. SPOTlight 50 was used to map the cell types found in the reference dataset to the spatial transcriptomics data. To train the model, up to 100 cells per cell type were used, with the gene set being the union between the marker genes of the cell types and the 3,000 most-highly variable genes. Marker genes for each cell type were obtained using the Seurat function FindAllMarkers, whereby only markers with positive log 2 -transformed fold change were considered. Cell types contributing <4% to the spot's predicted composition were considered to be fitting noise and were set to 0. To quantify the changes in tumour-associated Schwann cells across the healthy, tumour front and tumour, a pairwise Wilcoxon test was used between the regions; the Kruskal-Wallis test was used to compare the three distributions simultaneously.

RNA-seq analysis of in vitro-cultured OSCC cells
Sample processing. In vitro SCC-25 untreated/PA-treated cells at 4 d or 14 d after PA withdrawal were collected for FACS-sorting. We collected 500,000 CD44 + cells per condition and 10 µg of total RNA was purified using the RNeasy Mini Kit (74104, Qiagen) and further processed for mRNA-seq, applying the Illumina sequencing technology (https:// www.illumina.com/science/technology/next-generation-sequencing/ sequencing-technology.html).

ChIP-seq assays
For the 4 d chromatin time point, OSCC cells were either stimulated with PA (300 µM) or OA (50 µM) for 4 d in vitro, and chromatin was collected on the fourth day. For the 14 d time point, the corresponding FA was withdrawn from the medium after the 4 d stimulation, and cells were passaged normally for 2 weeks until the collection time for fixation. All of the ChIP-seq experiments targeting histone marks were performed using the ChIP-IT High Sensitivity kit from Active Motif (53040) according to the manufacturer's instructions. Sonication was performed for 30 min in a Bioruptor Pico (Diagenode). A list of the optimal conditions for each of the histone marks and antibodies used is provided in Supplementary Table 20c. Optimal chromatin fixation for MLL1 and MLL2 ChIPs consisted of double cross-linking using ChIP-IT High Sensitivity kit fixation reagents and the ChIP Cross-link Gold reagent from Diagenode (C01019027). Samples were sonicated for 15 min in a Bioruptor Pico (Diagenode) system. About 60 µg chromatin was used for both MLL1 and MLL2 IPs according to the instructions of the ChIP-IT High Sensitivity kit; anti-MLL1 and anti-MLL2 antibodies were produced in-house in the A. Shilatifard laboratory (Northwestern Medicine).
EGR2 ChIPs were performed according to the recommendations of the ChIP-IT High Sensitivity kit. Chromatin was sonicated for 15 min in a Bioruptor Pico (Diagenode), and 45 µg of chromatin and 8 µg of anti-EGR2 antibodies (Ab43020) were used in the immunoprecipitation step.
H3K4me3 ChIP-seq from sorted human tumour cells was performed using the True MicroChIP kit from Diagenode (C01010130). In brief, 30,000 cells per condition were FACS-sorted, fixed and sonicated together. Optimal sonication was achieved after 20 min in a Bioruptor Pico (Diagenode). Samples were then split in three, and each IP was performed from 10,000 cells using the ChIP-seq-grade H3K4me3 positive control antibodies provided. All libraries for sequencing were prepared using the NEBNext Ultra DNA Library Prep Kit from Illumina (E7370L) according to the manufacturer's instructions.
Peak calling to determine the enrichment over the background was performed using MACS (v.2.1.2) callpeak 54 . The broad flag method was used for H3K9me3, H3K4me1, H3K4me3 and Set1A with a q-value narrow and broad cut-off of 0.05, and for H3K27me3, MLL1 and MLL2, with a q-value narrow cut-off of 0.05 and a q-value broad cut-off of 0.0001. Signal tracks on ChIP-seq peak enrichment levels over the whole genome were generated using MACS (v.2.1.2) bdgcmp on the bed-Graph files generated by MACS (v.2.1.2) callpeak. For H3K27ac, peaks were called using GeneRich (v.0.6) and FDR ≤ 0.05. GeneRich was also used to remove duplicate reads. Differential peaks were determined using the DiffBind package (v.2.4.8) in R (v. 3.54.24) and FDR ≤ 0.05.
The distribution of the differential regions within genomic features was determined using the annotatePeaks.pl script of the HOMER suit of tools (v.4.10) 55 . R (v.3.4.4) 56 and the Bioconductor R package DiffBind 57 were used for differential binding analysis of H3K4me1, H3K4me3, H3K27me3, MLL1 and MLL2. The binding matrix was calculated with affinity scores based on trimmed mean of M values (TMM) normalization (EdgeR), using read counts minus control read counts and effective library size. PCA plots were generated using affinity data for all sites. Differential analysis using DESeq2 was used for H3K27ac, H3K4me3 (SCC-25 EGR2-KD cells) and EGR2 ChIP-seq datasets. The FDR threshold was fixed as either 0.1 or 0.05 (as specified in each figure legend).
Heat maps and average profiles of ChIP binding to TSS regions were plotted using the ChIPseeker R Bioconductor package 58 .

ChIP-seq data and RNA-seq data combined analysis
To assess transcription in the genomic loci bearing PA-driven H3K4me3 changes at the 4 d/14 d time points, the already-defined set of H3K4me3 peaks for each condition was used (as described in the ChIP-seq data processing and differential binding analysis section). The H3K4me3 peaks of all replicates were merged for each condition using the BEDTools Intersect (v.2.28.0) function. The BEDTools Window (v.2.28.0) function was then applied to overlap the TSS of all expressed genes in each experimental condition, as determined by the RNA-seq data (see the RNA-seq analysis of in vitro-cultured OSCC cells section), with all H3K4me3 peaks in that condition in a window of ±300 bp. Regions containing DEGs, as defined by RNA-seq using a fold-change cut-off > 2 and P < 0.05 were selected. The RNA-seq TMM of all of the DEGs in each time point (4 d or 4 d + 14 d) was used as the input for the corresponding box plots.

PRO-seq studies
The PRO-seq assays were performed at A. Shilatifard's laboratory according to a previously detailed protocol 61 . In brief, nuclei from OSCC cells were isolated using a dounce homogenizer with a loose pestle, and the nuclear run-on reaction was performed at 30 °C for 3 min using 25 µM biotin-11-ATP/GTP/CTP/UTP (Perkin Elmer). Total RNAs were extracted and hydrolysed using a 0.2 M NaOH solution for 10 min on ice. Biotinylated nascent RNAs were purified using Dynabeads M-280 Streptavidin (Invitrogen), followed by 3′-VRA3 and 5′-VRA5 adapter ligation and two additional rounds of purification with Dynabeads. cDNA was then generated by reverse transcription using SuperScript III (Invitrogen), and a ten-cycle amplification step of the cDNA libraries was performed. The indexed DNA libraries were size-selected (140-350 bp) and sequenced on the NextSeq 500 system (Illumina) with single-read runs. Two biological replicates per experimental condition were sequenced for the SCC-25 cell line.

PRO-seq data processing
For eRNA mapping, the BEDTools Window (v.2.28.0) function was used to overlap H3K27ac and H3K4me1 peaks for each time point (4 d or 4 d + 14 d) and each condition (PA or untreated) in a window of ±500 bp, after the merging of all of the common peaks using the BEDTools merge (v.2.28.0) function. Homer annotatePeaks (v.4.9.1-5) was used for peak annotation and only intergenic peaks were selected, as described previously 55 . The H3K27ac and H3K4me1 intergenic peaks were overlapped using the BEDTools Window (v.2.28.0) function in a window of ±250 bp, and a set of intergenic enhancers was defined. After identifying intergenic regions containing H3K27ac/H3K4me1 peaks, the BEDtools Intersect function was used with the option --v to discard peaks overlapping with any gene from the hg19 reference genome (ENSENBL v87), including an additional threshold of 2 kb upstream of every TSS. To determine which intergenic enhancers presented transcriptional activity, we separated the PRO-seq mapped reads that were separated by strands and individually called peaks using MACS2.1.2 with the parameters --nomodel --extsize 200 --broad-cutoff 0.1 and --broad to identify broad peaks. The common peaks were merged from each sample (that is, 4 d or 14 d) for each stranded peak; all peaks from the common plus strand were then merged with the common minus-strand peaks. BEDtools intersects between the intergenic enhancer peaks and the common peaks from PRO-seq were used to define enhancers that were actively transcribing enhancer RNAs (eRNAs). All transcripts shorter than 800 nucleotides were excluded to minimize any bias in the Pol II travelling ratio analysis (as these eRNAs are short) to precisely separate promoter and gene body regions. Finally, the TSS coordinates were defined for each eRNA based on the position of the highest signal occurrence in the PRO-seq peak window, and those with a counts per million read value of >0.5 were discarded. For the eRNA travelling ratio calculations, the promoter region was defined as −100 nucleotides from the PRO-seq summit to +150 nucleotides, and the gene body as +151 nucleotides from the PRO-seq summit to transcription ending site plus 1 kb, for all transcripts used in the analysis (n = 2,182). The BEDTools closest (v.2.28.0) function was used to identify the closest expressed upstream and downstream genes.

scRNA-seq assays using SMART-seq
Full-length single-cell RNA-seq libraries were prepared using the Smart-seq2 protocol 62 with minor modifications. In brief, freshly collected single cells were sorted into 96-well plates containing the lysis buffer (0.2% Triton-100, 1 U µl −1 RNase inhibitor). Reverse transcription was performed using SuperScript II (Thermo Fisher Scientific) in the presence of 1 µM oligo-dT30VN (IDT), 1 µM template-switching oligonucleotides (Qiagen), and 1 M betaine. cDNA was amplified using the KAPA Hifi Hotstart ReadyMix (Kapa Biosystems) and IS PCR primer, with 23 cycles of amplification. After purification with Agencourt Ampure XP beads (Beckmann Coulter), product size distribution and quantity were assessed on a Bioanalyzer using the High Sensitivity DNA Kit (Agilent Technologies). A total of 120 pg of the amplified cDNA was fragmented using Nextera XT (Illumina) and amplified with double-indexed Nextera PCR primers (IDT). The products of each well of the 96-well plate were pooled and purified twice with Agencourt Ampure XP beads. The final libraries were quantified and checked for fragment size distribution using the Bioanalyzer High Sensitivity DNA Kit. Pooled sequencing of Nextera libraries was performed using the HiSeq2500 (Illumina) system to an average sequencing depth of 50,000 reads per cell. Sequencing was carried out as paired-end (PE75) reads with library indexes corresponding to cell barcodes.
scRNA-seq analysis of SMART-seq data scRNA-seq data analysis was performed using the R package Seurat (v.2.3.0) 63 . Gene expression count matrices were filtered both to exclude cells with a mitochondrial percentage higher than 15, percentage of mapped reads lower than 85 and small library size (corresponding to the first quartile from total sum of counts across all features), and to exclude features that were not expressed in at least 5 cells and had a minimum of 90 counts. Gene expression levels for each cell were normalized to total expression, multiplied by a scale factor of 10,000 and log-transformed. Before clustering, highly variable genes, which represent 10% of total gene features, were identified on the basis of their average expression and dispersion and data were scaled regressing out batch effects. The number of significant principal components was assigned by the randomization approach 'jack straw', and clustering was performed on the top 10 principal components with a resolution of 0.5 by a shared nearest neighbour modularity optimization-based clustering algorithm. t-distributed stochastic neighbour embedding (t-SNE) was applied to display the data to cell loadings of ten selected principal components. Cells were scored by PA responder genes and projected onto the t-SNE. Gene expression markers were found for all of the clusters using Wilcoxon rank-sum tests with a return threshold of 0.01. To specifically find genes that were different between each time point and condition, Wilcoxon rank-sum tests were applied between untreated cells at 4 d and PA-treated cells at 4 d; only genes that were detected in a minimum fraction of 0.1 in either of the two populations were tested.
Trajectory analysis was performed using Monocle (v.2.10.1) 64 with the highly variable genes obtained in previous steps for each time point separately. The DDTree method was used for dimensionality reduction and batch effects were subtracted from the data to avoid their contribution to the trajectory. Finally, clusters and different treatment conditions were projected on the inferred trajectory.

Western blotting
After a wash with cold PBS, cells were lysed for 30 min at 4 °C using buffer containing 1× EDTA, 150 mM NaCl, 20 mM HEPES, supplemented with protease and phosphatase inhibitors.
All of the samples were pipetted several times to ensure proper lysis. Cell lysates were centrifuged at 4 °C for 15 min at 12,000 r.p.m., and the supernatant was retained as the protein extract. All of the samples were diluted in Laemmli buffer and boiled at 95 °C for 5 min. Equal amounts of protein for all of the samples were loaded into 6% polyacrylamide gels for protein separation, followed by transferring to LF-PVDF membranes for 90 min at 4 °C. Membrane blocking incubation was carried out with 5% milk-supplemented TBS for 30 min at room temperature. In-house-generated anti-MLL1, anti-MLL2, anti-Set1A or anti-Set1B antibodies, or loading control anti-Hsp90 antibodies (Ab13495, Abcam) were incubated at 4 °C overnight on a rocking platform. Next, membranes were incubated with anti-rabbit HRP-conjugated secondary antibodies at 1:5,000 dilution for 1 h at room temperature and washed. Blots were visualized using SuperSignal West Pico Chemiluminescent Substrate (Thermo Fisher Scientific) on the ChemiDoc Gel Digital Imaging System (BioRad).
For chABC western blot analysis, 5 µg total protein lysate from in vitro VDH-15 cells infected with LV-chABC was loaded into 12% polyacrylamide gel for protein separation and then transferred to LF-PVDF membrane for 90 min at 4 °C. Membrane blocking incubation was performed with 5% milk-supplemented TBS for 30 min at room temperature followed by incubation with anti-chABC antibodies (1E10) (NBP 96141AF594, Novus Biologicals). Membranes were then incubated with anti-rabbit HRP-conjugated secondary antibodies at 1:5,000 dilution for 1 h at room temperature and washed. Blots were visualized using SuperSignal West Pico Chemiluminescent Substrate (Thermo Fisher Scientific) in a ChemiDoc Gel Digital Imaging System (BioRad). Western blot analysis was performed according to standard protocols.

Gene overlap analysis
For all gene overlap analyses (performed on expression and/or ChIPseq data), a hypergeometric test was used, considering a background of 25,000 genes in total. Hypergeometric tests yielded a representation factor and P value for each comparison, indicative of the gene overlap enrichment and randomness, respectively.

Gene expression analysis
RNA isolation and cDNA amplification were performed as previously described 65 . In brief, RNA was isolated using magnetic beads (RNAClean XP beads, Agencourt). RNA was reverse-transcribed and amplified using Whole Transcriptome Amplification chemistry (WTA2, Sigma-Aldrich). For monitoring amplification, SYBR Green was added to the reaction, and the reactions were stopped when the SYBR Green signal reached a plateau. cDNA was purified using the PureLink Quick PCR Purification Kit (Invitrogen) and quantified using the Nanodrop ND-1000 spectrophotometer (Thermo Fisher Scientific).
Subsequently, 8 µg of cDNA were fragmented and labelled using the GeneChip Mapping 250K Nsp Assay Kit (Affymetrix, 900766) according to the manufacturer's instructions. Mouse MG-430 PM arrays were hybridized, washed, stained and scanned according to the protocol described in the GeneAtlas Hybridization, Wash, and Stain Kit for 3′ IVT Arrays user manual.
Samples that were ready for hybridization were denatured at 96 °C for 10 min before incubation on the Mouse MG-430 PM Array Strip. Hybridization was performed for 16 h at 45 °C in a GeneAtlas Hybridization Oven (Affymetrix). Washing and stain processing were performed at GeneAtlas Fluidics Station (Affymetrix), according to the specific script for Mouse MG-430 PM Arrays. Subsequently, arrays were scanned with GeneAtlas Scanner (Affymetrix). Normalized expression signals were calculated from Affymetrix CEL files using the RMA algorithm 66 . qPCR qPCR using TaqMan gene expression probes (Applied Biosystems) (Supplementary Table 20d) was performed and analysed using a 7900-HT Fast Real-Time PCR Instrument (Applied Biosystems). Relative expression levels were determined by normalization to beta-2-microglobulin (B2M) and/or beta-actin (ACTB) using the ΔΔC t method.
Gene expression analysis in the microarray data was performed using R 75 and Bioconductor 76 . Microarrays were processed using RMA normalization as implemented in the Bioconductor R package affy 77 . Data were annotated using probeset information provided by Affymetrix on its product support webpage (http://www.affymetrix.com/ support/ downloaded 10/05/2017). Standard quality controls were performed to identify abnormal samples 78 regarding (1) spatial artefacts in the hybridization process (scan images and pseudo-images from probe level models); (2) intensity dependences of differences between chips (MvA plots); (3) RNA quality (RNA digest plot); (4) global intensity levels (box plot of perfect match log-intensity distributions before and after normalization and RLE plots); (5) anomalous intensity profile compared with the rest of the samples (NUSE plots, principal component analyses); and (6) effect of quality metrics 79 on expression measures. For the CD36 bright /CD36 dim and the neural stromal fraction datasets, expression values from samples sharing the same hybridization batch were a priori corrected genewise by scanning the batch using a linear model, in which sample group was included as a covariate. Furthermore, in these datasets, hybridization batch was also a priori corrected using the same approach. Expression for each specimen was summarized genewise using the mean across its technical replicates.
Differential expression analysis was performed at the probeset level using moderated t-statistics by empirical Bayes shrinkage as implemented in the limma R package 80 , including hybridization batch (CD36 bright /CD36 dim and neural stromal fraction datasets) or scanning batch (shEGR2/shGAL datasets) as a covariate in the model. In the shEGR2 and shGAL human dataset, an additional control by RMA. IQR metric 79 was required and was therefore included in the model. The Benjamini and Hochberg correction method 81 was applied for multiple-contrast adjustment.
Functional enrichment analyses in the microarray data were performed using a modification of ROAST 82 -a rotation-based approach implemented in the R package limma 80 that is especially suitable for small-size experiments. Such modifications were implemented to accommodate the restandardized maxmean statistic proposed by ref. 83 , in the ROAST algorithm, to enable it for competitive testing 84 . When no biological replicates were available, functional enrichment was performed using the preranked version of the GSEA 69 , for which genes were ranked according to the base 2 logarithm of the fold-changes obtained in the differential expression analysis. In both cases, data were collapsed to the gene level by selecting the probeset showing the highest s.d. within each gene. Genesets for these analyses were derived from GO 85 , as collected in the R packages org.Hs.eg.db 86 and org.Mm.eg.db 87 .

Statistical analysis
For all of the experiments, an adequate sample size was determined on the basis of the results of pilot studies. No statistical methods were used to determine sample size. All mice that fulfilled proper experimental conditions during the experimental procedures were included in the analysis. On the basis of the results of pilot studies, homogeneous groups of males and females between 8 and 10 weeks, and their control littermates, were used for the experimental studies. No statistically significant differences were found with respect to the sex of mice, and no previously described randomization method was used. Data are generally shown as the mean ± s.e.m. Statistical significance was analysed using Prism 6 (GraphPad) using two-tailed t-tests and Fisher's exact tests. P ≤ 0.05 was considered to be statistically significant.

Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this paper.

Data availability
All raw data from gene expression RNA microarrays, ChIP-seq, RNA-seq, scRNA-seq, PRO-seq and spatial transcriptomics are available at the Gene Expression Omnibus (GEO) repository under accession code GSE148321. Source data are provided with this paper.

Code availability
Extended Data Fig. 1 | Palmitic acid, but not linoleic, oleic or  The BLI signal is expressed as the relative normalized photon flux. Data are given as the mean and s.e.m. (untreated, n = 14; 14D after PA wdl, n = 14, *P = 0.05, two-tailed t-test). f, Frequency of developed LN metastases from animals in e (*P = 0.02, two-tailed Fisher's exact test). g, ex vivo BLI lung metastasis quantification of mice injected with VDH-15 cells at 14D fatty acid withdrawal (PA, palmitic acid; SA, stearic acid; OA, oleic acid; LA, linoleic acid). BLI signals are expressed as the relative normalized photon flux. Images are representative of n = 15 mice per group. h, Frequency of developed lung metastasis of mice injected with SCC-25 after FA withdrawal expressed as percentages (n = 10 mice per group, lung metastases: *P = 0.05; two-tailed Fisher's exact test). PA/OA denotes a 4-day treatment with palmitic acid followed by a 4-days treatment with oleic acid, followed by 14-day with no fatty acid. i, BLI quantification of FACS-sorted and serially-transplanted PT populations (as indicated) from SCC-25 primary recipients. BLI signals are shown as the normalized photon flux (n = 10 mice per group. *P = 0.03, **P = 0.05, ***P = 0.01; two-tailed t-test). j, Frequency of developed LN metastases of mice in i (*P = 0.05, two-tailed Fisher's exact test). Fig. 2 | Tumour cells from palm-oil-fed primary recipient mice display a prometastatic memory in secondary recipient mice in a CD36 dependent manner. a, Schematic diagram representing the in vivo diet experiments in which OSCC that were exposed to a high-fat diet in primary NSG mice recipients were injected into secondary recipients that were maintained on a control diet. At the final point primary tumour (PT) cells were purified by FACS-SORT for molecular characterization (created with BioRender.com). b, Flow cytometry analysis of PTs from VDH-15-and SCC-25-injected primary recipients fed a high-fat or control diet, as schematized in a. Numbers indicate [CD44 bright CD36 bright ] and [CD44 bright CD36 dim ] populations in the represented gate, expressed as percentages of total DAPI-GFP+Lin-cells. Samples are representative of n = 3 independent experiments. c, d, BLI tumour monitoring of secondary recipients injected with VDH-15-derived cells from primary recipient mice fed with palm oil-rich, olive oil-rich diet and normal diet (as control). BLI signals are shown as the normalized photon flux (n = 20 mice per group, two independent experiments. LN metastasis, P = 0.04; lung metastasis, *P = 0.04, **P = 0.007; two-tailed t-test). e, Frequency of developed LN metastases from animals in c expressed as percentage (for LN metastasis, *P = 0.003; **P = 0.02; for lung metastasis: *P = 0.01, ***P = 0.0004, two-tailed Fisher's exact test). e, f, BLI metastasis quantification of secondary recipient mice injected with control pLKO.1 or shCD36 VDH-15-(f) and SCC-25-(g) derived cells from primary recipient mice fed with a high-fat palm oil or a control diet. BLI signals are expressed as the normalized photon flux (for VDH-15, n = 10 mice per group; for lymph node and lung met, *P = 0.04; n.s., not significant; two-tailed t-test, data are mean ± s.e.m.; for SCC-25, n = 20 mice per group, data from two independent experiments; *P = 0.03, ****P < 0.0001; two-tailed t-test, data are mean ± s.e.m.). g, Frequency of developed metastases from animals in g expressed as percentage (for lymph node met, ***P = 0.003, **P = 0.0004; for lung met, ****P < 0.0001, two-tailed Fisher's exact test). h, Ex vivo BLI lung metastasis monitoring from secondary recipient mice injected with control pLKO.1 or shCD36 melanoma-derived cells from primary recipient mice fed with a high-fat palm oil or a control diet. Pictures are representative of n = 5 mice per group. i, BLI metastasis quantification of animals in i. BLI signals are expressed as the normalized photon flux. (***P = 0.0009, ****P <0.0001, twotailed t-test, data are mean ± s.e.m.). Fig. 3 | CD36 blocking metastatic assays and assessment of the impact of PA on the OSCC chromatin landscape. a, Frequency of developed lung metastasis from secondary mice injected with SCC-25-derived from primary recipient mice fed a normal or palm oil-rich diet; Primary recipients were treated with an anti-CD36 neutralizing antibody (JC63.1) or a control isotype (IgA). Data is expressed as the percentages. (n = 13 mice per group; *P = 0.01, *P = 0.03, two-tailed Fisher's exact test). b, c, ex vivo BLI imaging and quantification of lung metastasis from VDH-15 secondary recipients injected with a doxycycline-inducible shCD36 (shCD36 #23 or #76) VDH-15-derived tumour cells from primary recipient mice fed a normal or palm oil-rich diet. The secondary mice were untreated or continuously doxycyclinetreated. (Data are mean ± s.e.m. Images are representative of n = 10 mice per group; *P = 0.01; **P = 0.001; ****P < 0.0001; two-tailed t-test).  Fig. 7 | Impact of chemical sympathectomy on metastasis and identification of EGR2 and galanin signalling in tumour cells with a Palm Diet-induced metastatic memory. a, Diagram of the experimental setup. Secondary (2ary) recipient mice were injected on day 0 with OSCC from primary tumours (PTs) ± palm oil memory. 2ary recipients were treated with the neurotoxin 6-hydroxydopamine (6-OHDA; to induce the apoptosis of dopaminergic neurons) or vehicle on days -6, -3 and +3 relative to the day of injection. b, Frequency of developed lung metastasis from animals treated as shown in a, expressed as the percentages (*P = 0.03, **P = 0.001; two-tailed Fisher's exact test). c, Tyrosine hydroxylase (TH) immunofluorescence analysis of VDH-15 primary tumours (PTs) from secondary recipients ± palm oil memory after vehicle or 6-OHDA treatment. Nuclei are stained with DAPI. KT-14, cytokeratin 14. Note that 6-OHDA-lesioned-tumours display a marked reduction in the expression level of TH. Images are representative of n = 3 biological replicates per group. d, Top common predicted binding site motifs in promoter regions of the co-regulated neural-related genes in SCC-25 or VDH-15 primary tumours with a palm memory. Z-scores and P values are shown for each cell line. e, Integrative gene set enrichment analysis (GSEA) from PTs of secondary recipients SCC-25-, VDH-15-or melanoma tumour-derived cells ± palm oil memory. The graph shows the biological processes enrichment in palm oil memory compared to control diet. On the right side, detail of the neuropeptide signalling pathway-associated gene-set assayed by integrative GSEA showing galanin (GAL) as the top significantly represented gene in palm oil memory tumours. f, g, PCA plots of differentially bound signal (DBS) regions detected for 4D Untreated/PA-treated EGR2 ChIP-seq samples in f) SCC-25 pLKO.1 and g) SCC-25 CD36-KD cells. h, Heat maps showing the EGR2 DBS regions (differential peaks) detected for 4D SCC-25 pLKO.1 (left) and SCC-25 CD36 KD (right) shown in f) and g) respectively (FDR<0.05). Fig. 11 | Tumour cells with a Palm Diet-induced metastatic memory alter the tumour-associated stroma in a CD36-dependent manner. a, Flow cytometry analysis of SCC-25-primary tumours (PTs) from secondary (2ary) recipients ± palm oil memory. [GFP-/CD31-/CD45-] tumour stroma-selected cells were purified and processed to perform RNA microarrays. Data are representative of n = 3 independent experiments. Numbers indicate the population in the represented gate, expressed as the percentages. b, Principal component analysis (PCA) of the microarray data of VDH-15-associated bulk stroma purified from secondary (2ary) recipients. Diet samples are indicated as (B) primary recipient mice fed a normal (red), olive oil-rich (green) or palm oil-rich (blue) diet, and 2ary mice fed with a normal diet and (A) primary mice fed with a normal diet and 2ary mice fed with a palm oil-rich diet (fuchsia). The axis shows the percent variability covered by each of the represented components. c, Principal component analysis (PCA) from microarray data of 2ary recipient tumour-associated stroma derived from injection in primary recipients of control (pLKO.1)-or shCD36-VDH-15 or SCC-25 cells, as indicated by the colour-code. The axis shows the percent variability covered by each of the represented components. d, Volcano plots for the differential gene expression microarray analysis and GO analysis from a purified tumour stroma from 2ary recipient mice injected with OSCC from PTs ± palm oil memory. Volcano plots show the up-or downregulated genes in control-palm oil memory tumours compared to control-normal diet tumours (top) and the expression of these genes in the shCD36-palm oil (lower plot). Data from n = 2 biological replicates; fold-change, FC > 2, P < 0.01. The GO analysis shows the top biological processes categories upregulated in the palm oil memory tumour stroma, FC > 2, P < 0.05. (Integrated analysis, derived from VDH-15 or SCC-25 cells). e, Neural-mouse stromal signature induced in 2ary recipients after orthotopic injection of primary recipient-PT cells, derived from control (pLKO.1), shRNA-CD36, shRNA-EGR2 or shRNA-GAL SCC-25 cell. Each dot in the plot represents a neural gene. f, g, GSEA showing the negative enrichment of biological processes in neural-mouse stroma associated with 2ary PTs from shEGR2/PALM or shGAL/PALM vs pLKO.1/control diet conditions. Fig. 12 | Schwann cell development and neuron projection regeneration are processes enhanced in the stroma of palm-oil derived tumours. a, Functional network interaction between tumour-mouse stromal palm oil (PALM)-controlled genes and the human OSCC primary tumours with a palm memory. The network graph shows the co-regulated functional nodes of interaction between the two compartments (fold-change, FC > 1.5; P value < 0.05). b, c, 10X single-cell (sc) RNA-seq clustering analysis of the tumour-associated stroma purified from secondary (2ary) oral SCC-25 PTs ±palm oil memory.

Extended Data
The principal component UMAP plot is shown in which the specific cell types have been annotated to each respective cluster. d, Overlap analysis represented by Venn diagrams showing the intersection between bulk-stroma palm memory-controlled genes and 10X single-cell (sc) RNA-seq clusters from SCC-25-derived stromal cells in 2ary recipients. Representation factor (RF) and P values are shown for the overlap. Hypergeometric test; estimated number of protein-encoding genes = 25,000. Fig. 14 | Spatial transcriptomic analysis of secondary recipient OSCC primary tumours. a, Images showing Haematoxylin-Eosin (H&E) staining of primary tumours (PTs) from secondary (2ary) recipients injected with control (pLKO.1) and shCD36 SCC-25 cells derived from primary recipient mice fed a normal diet or a palm oil-rich diet. e, Spatial transcriptomics analysis of tumours in a) showing the proportion content and spatial distribution of mouse and human transcriptome per spot, the analysis stratification of the tissue as healthy/tumour invasive front/tumour and the proportion content and spatial distribution of the tumour-associated Schwann cells. c, Quantitative analysis of the proportional content of tumour-associated Schwann cells within each compartment (healthy/tumour front/tumour) of the PTs from 2ary recipients injected with control (pLKO.1) and shCD36 SCC-25 cells derived from primary recipient mice fed a normal diet or a palm oil-rich diet. Fig. 15 | Immunofluorescent analysis of secondary recipient OSCC primary tumours. a, b. e, Immunofluorescence analysis and quantification of primary tumours (PTs) from secondary (2ary) recipient mice injected with VDH-15-or SCC-25-derived cells, as indicated, derived from primary recipient mice fed a normal diet, a palm oil-rich or an olive oil-rich diet. Cytokeratin-14 (KT-14) is shown as an epithelial marker of OSCC, and the glial/Schwann markers s100 or GAP43 as a tumour stroma marker. Nuclei are stained with DAPI. Images are representative of n = 3 independent experiments. c, Graphs showing the values of integrated density of (a, b). (n = 3 biological replicates per group; *P = 0.05, ****P < 0.0001; two-tailed t-test). f, Immunofluorescence analysis of PTs from M15-treated or vehicle-treated mice ± palm oil memory. The expression of cytokeratin-14 (KT-14) is shown as an epithelial marker of OSCC, and the glial Schwann cell marker s100 is shown in the tumour-associated stroma. Nuclei are stained with DAPI. Note that the increased expression of s100 in the stroma of palm memory tumours is slightly reduced in the case of M15-treated palm condition.

Corresponding author(s): Benitah
Last updated by author(s): Sep 7, 2021 Reporting Summary Nature Research wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Research policies, see Authors & Referees and the Editorial Policy Checklist.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted Give P values as exact values whenever suitable.

For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings
For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes Estimates of effect sizes (e.g. Cohen's d, Pearson's r), indicating how they were calculated Our web collection on statistics for biologists contains articles on many of the points above.

Software and code
Policy information about availability of computer code

Data collection
Provide a description of all commercial, open source and custom code used to collect the data in this study, specifying the version used OR state that no software was used.

Data analysis
Provide a description of all commercial, open source and custom code used to analyse the data in this study, specifying the version used OR state that no software was used.
For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors/reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information.

Data
Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A list of figures that have associated raw data -A description of any restrictions on data availability nature research | reporting summary Mycoplasma contamination All cell lines tested negative for mycoplasma contamination.

Commonly misidentified lines (See ICLAC register)
Name any commonly misidentified cell lines used in the study and provide a rationale for their use.

Palaeontology Specimen provenance
Provide provenance information for specimens and describe permits that were obtained for the work (including the name of the issuing authority, the date of issue, and any identifying information).

Specimen deposition
Indicate where the specimens have been deposited to permit free access by other researchers.

Dating methods
If new dates are provided, describe how they were obtained (e.g. collection, storage, sample pretreatment and measurement), where they were obtained (i.e. lab name), the calibration program and the protocol for quality assurance OR state that no new dates are provided.
Tick this box to confirm that the raw and calibrated dates are available in the paper or in Supplementary Information.

Animals and other organisms
Policy information about studies involving animals; ARRIVE guidelines recommended for reporting animal research Laboratory animals NOD scid gamma (NSG) (NOD.Cg-Prkdcscid Il2rgtm1Wjl/SzJ) mice were purchased from Charles River and crossed in-house.

Wild animals
Provide details on animals observed in or captured in the field; report species, sex and age where possible. Describe how animals were caught and transported and what happened to captive animals after the study (if killed, explain why and describe method; if released, say where and when) OR state that the study did not involve wild animals.

Field-collected samples
All mice (a mix of male and female) were housed under a regimen of 12 h light/12 h dark cycles and SPF conditions. All tissue samples were collected before ethical final point.

Ethics oversight
All procedures were evaluated and approved by the CEEA (Ethical Committee for Animal Experimentation) from the Government of Catalunya.
Note that full information on the approval of the study protocol must also be provided in the manuscript.

Human research participants
Policy information about studies involving human research participants

Recruitment
Describe how participants were recruited. Outline any potential self-selection bias or other biases that may be present and how these are likely to impact results.

Ethics oversight
The study followed the guidelines of the Declaration of Helsinki, and patient identity and pathological specimens remained anonymous (in the study and elsewise).
Note that full information on the approval of the study protocol must also be provided in the manuscript.

Clinical data Policy information about clinical studies
All manuscripts should comply with the ICMJE guidelines for publication of clinical research and a completed CONSORT checklist must be included with all submissions.

Clinical trial registration
Provide the trial registration number from ClinicalTrials.gov or an equivalent agency.

Study protocol
Note where the full trial protocol can be accessed OR if not available, explain why.

Data collection
Describe the settings and locales of data collection, noting the time periods of recruitment and data collection.

Outcomes
Describe how you pre-defined primary and secondary outcome measures and how you assessed these measures.