A unique subset of glycolytic tumor propagating cells drives squamous cell carcinoma

Head and Neck Squamous Cell Carcinoma (HNSCC) remains among the most aggressive human cancers. Tumor progression and aggressiveness in SCC are largely driven by Tumor Propagating Cells (TPCs). Aerobic glycolysis, also known as the Warburg Effect, represents a characteristic of many cancers, yet whether this adaptation is functionally important in SCC, and at which stage, remains poorly understood. Here, we show that the NAD+-dependent histone deacetylase Sirtuin 6 (SIRT6) is a robust tumor suppressor in SCC, acting as a modulator of glycolysis in these tumors. Remarkably, rather than a late adaptation, we find enhanced glycolysis specifically in TPCs. More importantly, using single cell RNA sequencing of TPCs, we identify a subset of TPCs with higher glycolysis and enhanced pentose phosphate pathway and glutathione metabolism, characteristics that are strongly associated with a better antioxidant response. Altogether, our studies uncover enhanced glycolysis as a main driver in SCC, and, more importantly, identify a subset of TPCs as the cell-of-origin for the Warburg effect, defining metabolism as a key feature of intra-tumor heterogeneity.

T PCs (or cancer stem cells) in SCC are responsible for sustaining primary tumours and are able to re-populate the entire tumour after transplantation, owing to their self-renewal and differentiation capacity 1 . As such, TPCs have emerged as attractive therapeutic targets 2 . Although genetic drivers, such as the surface marker CD34 and the transcription factor SOX2 (refs. [3][4][5][6], have been identified in TPCs, the specific metabolic characteristics of these cells remain poorly investigated. Metabolic reprogramming has emerged as a critical hallmark of cancer 7,8 . In particular, increased glycolysis and lactate production under normoxia (aerobic glycolysis) are among the best-described characteristics of many tumours. Such adaptation, also known as the Warburg effect, provides transformed cells with intermediate metabolites for biomass while balancing cellular redox status for continuous proliferation 9,10 . Sirtuin 6 (SIRT6), a member of the NAD + -dependent protein deacetylases, negatively regulates HIF-1α-dependent glycolysis gene expression (for example, the glucose transporter GLUT1, pyruvate dehydrogenase kinase 1 (PDK1), and lactate dehydrogenase-A (LDHA) as an H3K9/H3K56 deacetylase, affecting glucose homeostasis 11,12 . SIRT6-deficient cells exhibit aggressive tumour formation through enhanced aerobic glycolysis in colon cancer 13 , and increased expression of oncofetal proteins in pancreatic cancer 14 , emphasizing a pivotal role for SIRT6 in glucose metabolism and tumorigenesis.
In this study, we find that enhanced glycolysis in a model of SIRT6 loss enriches for CD34 + TPCs, in turn resulting in a much more aggressive tumorigenic phenotype. Mechanistically, highly glycolytic CD34 + TPCs present a distinct gene signature associated with glutathione (GSH) metabolism and stemness, thereby providing a defence against oxidative stress, which is robustly enhanced on Sirt6 loss. Using metabolite profiling analysis, we further demonstrate that generation of antioxidants and nucleotides through the oxidative phase of the pentose phosphate pathway (oxPPP) is largely responsible for the aggressive tumorigenic phenotype. Remarkably, direct metabolite measurement from in vivo tumour samples in cellular spatial resolution by matrix-assisted laser desorption/ ionization (MALDI) mass spectrometry imaging (MSI) indicates higher glycolysis and more reduced GSH in CD34 + TPCs compared to CD34 − tumour cells. Further, single-cell RNA-sequencing (scRNA-seq) analysis defines a subset of TPCs with such characteristics that are functionally crucial for TPC enrichment and tumourigenic potential. To our knowledge, our studies provide A unique subset of glycolytic tumour-propagating cells drives squamous cell carcinoma the first in-depth characterization of the metabolic adaptations in SCC, identifying a previously unrecognized metabolic heterogeneity within TPCs, a key feature to support antioxidant protection and nucleotide synthesis in these unique tumour-driving cells.

SIRT6 acts as a tumour suppressor in squamous cell carcinoma by modulating aerobic glycolysis.
We sought to define an in vivo model of SCC to closely examine the effect of increased glycolysis in tumourigenesis and its specific subpopulations including TPCs. We reasoned that SIRT6 might act as a tumour suppressor in SCC, one of the major types of epithelial cancers, via modulation of glycolysis. Notably, SIRT6 copy number loss is associated with shorter overall survival when analysed in samples from individuals with head and neck squamous cell carcinoma (HNSCC) in The Cancer Genome Atlas (TCGA; Extended Data Fig. 1a). Further, both SIRT6 copy number and expression were significantly decreased in HNSCC samples compared to matched normal tissue in either the Oncomine platform or the TCGA (Extended Data Fig. 1b,c). SIRT6 expression was already observed in early-stage tumours, thus suggesting that SIRT6 loss may be functionally important in both initiation and maintenance of SCCs (Extended Data Fig. 1b). We next analysed SIRT6 protein expression in human samples from individuals with HNSCC and normal skin tissues by immunohistochemistry and found that less differentiated tumours tended to have less SIRT6 expression (Extended Data Fig. 1d). Last, among available human cancer cell lines listed in the Cancer Cell Line Encyclopedia (CCLE), almost all HNSCC cell lines (denoted as upper aerodigestive tract) exhibited SIRT6 loss (Extended Data Fig. 1e). Altogether, these analyses suggest a potential tumour suppressive role for SIRT6 in SCC.
To better define the roles of Sirt6 in squamous cell carcinogenesis, we generated an in vivo Sirt6 conditional knockout (cKO) mouse model (Sirt6 F/F ;K14-cre + ) that specifically deletes Sirt6 in the skin epithelium. These mice, along with wild-type (WT) animals (Sirt6 +/+ ;K14-cre + ), were treated with 7,12-dimethylbenz(a) anthracene (DMBA), a known carcinogen, followed by repetitive 12-O-tetradecanoylphorbol-13-acetate (TPA) treatment for 14 weeks, a well-established protocol to recapitulate SCC in vivo (Fig. 1a). Remarkably, Sirt6 cKO animals showed an earlier onset of tumours (Fig. 1b) and significantly larger tumours at 14 weeks after DMBA treatment (Fig. 1c). The C57BL/6 strain is known to be highly resistant to tumourigenesis in vivo, and without continuous TPA treatment, existing skin tumours tend to regress 15 . Consistently, most of the WT tumours became smaller and regressed after discontinuing TPA treatment for more than 7 weeks. In contrast, multiple Sirt6-deleted tumours remained and even grew larger (Fig. 1d). Notably, fully transformed SCC were exclusively formed in Sirt6 cKO animals (Extended Data Fig. 2a). We next assessed tumour cell proliferation and detected a major increase in cells positively staining for proliferating cell nuclear antigen (PCNA + ) in Sirt6 cKO tumours (Fig. 1e). Overall, these data suggest that loss of Sirt6 pro-motes tumour cell proliferation, resulting in enhanced tumour progression and maintenance.
To directly determine whether increased glycolysis plays a functional role in driving tumour progression, we next administered an inhibitor of glycolysis, dichloroacetate (DCA), in drinking water during the DMBA/TPA treatment in vivo (Fig. 1f). Inhibition of glycolysis significantly reduced tumour size in Sirt6 cKO animals (Fig. 1g), while continuous inhibition of glycolysis following TPA withdrawal completely impaired the progression of the existing tumours in Sirt6-deleted animals (Fig. 1h). These findings emphasize a critical role for aerobic glycolysis in SCC growth and maintenance. Molecularly, dysplastic proliferating cells started to express high levels of GLUT1, while normal proliferating (PCNA + ) keratinocytes and hair follicular stem cells barely expressed GLUT1, suggesting that increased glycolytic metabolism is an important adaptation of transformed cells, not normal proliferating cells (Extended Data Fig. 2b). In line with these results, RNA expression of Glut1 and another glycolytic gene, Ldha, was also increased in tumours compared to adjacent normal skin (Extended Data Fig. 2c). Further, GLUT1 appeared to be specifically expressed in tumours, particularly in tumour basal layers, compared to normal adjacent skin ( Fig. 1i and Extended Data Fig. 2d). This was also true for the PDK-dependent phosphorylated form of pyruvate dehydrogenase complex (phospho-PDH; Fig. 1i and Extended Data Fig. 2d), the rate-limiting enzyme that converts pyruvate into acetyl-CoA for usage in the tricarboxylic acid cycle. Phosphorylation of PDH inactivates the enzyme, forcing production of lactate instead. Mitochondria pyruvate carrier 1 (MPC1) expression was higher in more differentiated tumour cells compared to basal tumour cells (Extended Data Fig. 2d), suggesting less pyruvate entry into mitochondria in these cells. Together, these results indicate that tumour basal cells are exquisitely glycolytic.
We next analysed expression of SIRT6 and glycolytic genes in samples of SCCs from individuals. Several glycolytic genes were inversely correlated with the level of SIRT6, a phenotype observed even at early stages of carcinogenesis (Extended Data Fig. 2e). Last, SIRT6 low, less differentiated tumour cells tend to express ubiquitously higher levels of GLUT1 (Extended Data Fig. 2f). This further strengthens our mouse data and validates this in vivo model as highly relevant to human HNSCC.
Increased glycolysis from Sirt6 loss enriches the number of tumour-propagating cells. Next, we sought to follow these GLUT1 + tumour cells and to assess their functional role in tumourigenesis. First, we analysed the differentiation state of GLUT1 + cells in the tumours by using keratin 5 and keratin 10, markers of basal progenitors and differentiated cells in skin, respectively. Most of the GLUT1 + cells co-stained with keratin 5 and were mutually exclusive with keratin 10 (Extended Data Fig. 3a), suggesting that less differentiated cells are particularly glycolytic. We next stained for the surface protein CD34, an established marker of TPCs. As expected, most of CD34 + cells co-expressed keratin 5 and were negative Fig. 1 | Sirt6 acts as a tumour suppressor in squamous cell carcinoma by negatively regulating aerobic glycolysis. a, DMBA/TPA-induced skin carcinogenesis protocol in Sirt6 WT or Sirt6 cKO animals. b, Tumour-free period after starting DMBA treatment in Sirt6 WT or Sirt6 cKO animals. Statistical analysis was done by log-rank test. c, Tumour size was measured at 14 weeks after DMBA treatment. Data are the mean ± standard deviation (s.d.). d, Tumour progression was assessed after stopping TPA treatment (at 14 weeks after DMBA) for at least 7 weeks. Fisher's exact test was performed for statistical analysis (P < 0.0001; two-sided). e, PCNA immunostaining in DMBA/TPA-treated skin tumours from Sirt6 WT or Sirt6 cKO animals. Representative images (scale bars indicate 100 µm) and quantification of PCNA + layers from normal adjacent skin and skin tumours. IHC, immunohistochemistry. Data are the mean ± s.d. f, Schematic of DCA treatment in DMBA/TPA-treated animals. DCA was administered at 7-8 weeks after DMBA treatment, to avoid any confounding effect of DCA on tumour initiation. g, Tumour size was measured at 14 weeks after DMBA treatment. Data are the mean ± s.d. h, Tumour progression was assessed after stopping TPA treatment with continuous DCA treatment. Data from the first two groups in d were used again for comparison. Fisher's exact test was performed for statistical analysis (P < 0.0001; two-sided). i, GLUT1 and phospho-PDH (Ser293) immunostaining in Sirt6-deleted large papilloma samples. Scale bars indicate 100 µm. Statistics, sample sizes (n) and numbers of replications are presented in 'Statistics and reproducibility'. *P < 0.05, **P < 0.01, ***P < 0.001. for keratin 10 (Extended Data Fig. 3b). Remarkably, most of the GLUT1 + cells were positive for CD34 and SOX9, supporting the idea that glycolytic basal cells are putative TPCs (Fig. 2a). Importantly, CD34 + SOX9 + hair follicular stem cells were negative for GLUT1, indicating that specifically the CD34 + tumour cells, rather than normal skin stem cells, benefit from enhanced glucose uptake  (Fig. 2a). This finding is further strengthened by co-staining with GLUT1, CD34 and keratin 10 in whole-tumour samples and subsequent calculation of correlation values (Extended Data Fig. 3c). Similarly, α6-integrin high /CD34 + cells have much higher expression levels of GLUT1 compared to α6-integrin high /CD34 − cells (Extended Data Fig. 3d).
Although markers used to identify TPCs in human HNSCC remain controversial, TPCs in murine cutaneous SCC have been well defined in the past decade [3][4][5][6][16][17][18][19] . We attempted to analyse and prospectively isolate TPCs (α6-integrin high /CD34 + ) from Sirt6 WT or Sirt6-deleted skin tumours (Extended Data Fig. 4a). Strikingly, when we analysed the percentage of TPCs in live lineage-selected tumour cells negative for propidium iodide (PI) and positive for yellow fluorescent protein (YFP; Methods), we found that Sirt6-deleted tumours exhibited a significant increase in TPCs ( Fig. 2b and Extended Data Fig. 4b), especially in tumours that were large (more than 2.5 mm in size), compared to size-matched tumours from WT mice. We next inhibited glycolysis in vivo by treating the animals with DCA, which caused significantly decreased blood lactate and reduced phospho-PDH (Ser293) levels (Extended Data Fig. 4c,d). Both results confirm successful inhibition of glycolysis. More importantly, DCA treatment severely reduced the percentage of TPCs both in WT and Sirt6 cKO tumours ( Fig. 2c and Extended Data Fig. 4b), suggesting that enhanced glycolysis could provide a unique advantage to TPCs, a phenotype exacerbated in the absence of SIRT6. Together, these results indicate that the increased tumour growth and maintenance phenotype observed in Sirt6-deleted tumours could be due to an increase in TPCs, a population uniquely glycolytic.

TPCs exhibited higher glutathione metabolism and better antioxidant response.
To pinpoint mechanistic pathways that could explain the benefit of enhanced glycolysis in TPCs, we performed RNA sequencing (Extended Data Fig. 4e,f). First, we established a 'common TPC gene signature' . In this analysis, we identified 397 commonly upregulated genes and 191 commonly downregulated genes ( Fig. 2d and Supplementary Table 1). DAVID pathway analysis in these upregulated genes from differentially expressed genes (DEGs) of Sirt6 cKO TPCs compared with Sirt6 cKO α6 high / CD34 − cells (Supplementary Table 2) and DEGs of WT TPCs compared with WT α6 high /CD34 − cells (Supplementary Table 3) revealed upregulation of several pathways in TPCs, including GSH metabolism, lipid metabolism, amino acid transport and multicellular organism development ( Fig. 2e and Extended Data Fig. 4g). Some of the commonly upregulated genes are already known to be important in TPCs (for example, Sox2, Ptk2 and Ereg), supporting the quality of our data 5,6,16,17 .
Enhanced GSH metabolism is important for antioxidant defence, and lipid metabolism and amino acid transport are vital for cellular energy and biomass, suggesting that rewiring metabolism is pivotal in TPCs. Consistently, Sirt6 cKO TPCs exhibited higher expression of genes in glycolysis and the PPP, consistent with their aggressive phenotype (Supplementary Table 4 and Fig. 2f). Moreover, genes involved in both GSH metabolism and redox balance, which were enriched in TPCs, were even higher in Sirt6 cKO TPCs (Fig. 2f). Strikingly, genes involved in stemness and carcinogenesis were expressed at a much higher level in Sirt6-deleted TPCs compared to WT TPCs, providing further rationale for the increased aggressiveness in the Sirt6-deleted tumours (Extended Data Fig. 4h).
Increased oxidative pentose phosphate pathway, generation of reduced glutathione and nucleotides in squamous cell carcinoma. To gain further insights into the metabolic adaptations that can be regulated by SIRT6, we took advantage of two human SCC cell lines. In HSC2 cells that barely express SIRT6, we found that ectopic expression of SIRT6 causes repression of glycolytic gene expression, reduced H3K9/K56 acetylation (Extended Data Fig. 5a) and decreased glycolytic reserve capacity (Extended Data Fig. 5b). Using stable-isotope tracing with [U-13 C]glucose, we found diminished glycolytic flux towards glycolytic intermediates including fructose-6-phosphate (F6P), pyruvate and lactate ( Fig. 3a and Extended Data Fig. 5c), and delayed 13 C incorporation into ribose-5-phospohate (R5P; Fig. 3a), at a time when most of the cells (both in WT and SIRT6-H133Y) were alive (Extended Data Fig. 5d). Of note, although citrate labelling did not reach steady state, the flux into citrate from glycolysis via acetyl-CoA (M + 2) showed similar labelling kinetics (Extended Data Fig. 5e). In addition, the flux into α-ketoglutarate from glycolysis (M + 2) that reached a pseudo-steady state was not affected, indicating that cells sustained normal mitochondrial respiration (Extended Data Fig. 5e). Decrease in glycolysis from SIRT6 overexpression was dependent on SIRT6 enzymatic activity, since expression of the SIRT6-H133Y catalytically inactive mutant (HY) did not influence glycolysis (Fig. 3a,b and Extended Data Fig. 5a-c). Although these experiments suggest that the enzymatic activity of SIRT6 is necessary, we cannot rule out additional, nonenzymatic roles for SIRT6 in this phenotype.
The PPP consists of an oxidative phase (critical to regenerate GSH) and a non-oxidative phase, which can be distinguished with [1,2-13 C]glucose (Fig. 3b). Interestingly, the relative fraction of M + 1 R5P was appreciably decreased when SIRT6 was overexpressed, while the relative fraction of M + 2 R5P remained similar, indicating that the oxidative arm of the PPP was strongly affected (Fig. 3b). In addition, sedoheptulose-7-phosphate (S7P), one of the intermediates in the non-oxidative PPP, did not show a notable difference in 13 C incorporation (Fig. 3a,b), further confirming our findings. Last, 13 C incorporation into DNA (M + 5 isotopologue) was significantly reduced in these cells (S6WT), indicating that enhanced PPP in SIRT6-deficient tumours serves as a critical precursor in de novo nucleotide synthesis ( Fig. 3c; pool size data are available in Supplementary Table 5).
As a complementary approach, we used SCC13, a skin SCC cell line that endogenously expresses high levels of SIRT6. We observed that inducible knockdown of SIRT6 caused increased glycolytic gene expression, increased bulked H3K56Ac and increased H3K9 acetylation specifically in glycolytic genes (Extended Data Fig. 6a,b). Consistently, a global increase in H3K56Ac was also observed in tumour samples in vivo from Sirt6 cKO animals compared to WT animals (Extended Data Fig. 6c). Intriguingly, CD34 + TPCs in WT tumour samples showed higher H3K56Ac levels compared to those in CD34 − tumour cells, indicating that SIRT6 activity might be reduced in CD34 + TPCs, likely contributing to the enhanced glycolytic phenotype in TPCs (Extended Data Fig. 6d). Knockdown of SIRT6 also boosted glucose uptake (consumption), lactate secretion and glycolytic capacity (Extended Data Fig. 6e- while enriching the relative abundance of several glycolytic intermediates including pyruvate and lactate (Extended Data Fig. 6l). These results are indicative of a quantitative increase in glucose metabolism and glycolysis, rather than an imbalance between glycolysis and mitochondrial respiration, since the ratios of lactate secretion/glucose consumption and NAD + /NADH did not show     Fig. 6h,j). Notably, highly glycolytic cells (shSIRT6) exhibited more GSH and less GSSG/GSH (oxidized versus reduced forms of GSH) and lower reactive oxygen species ( Fig. 3d and Extended Data Fig. 6k). We also observed increased levels of several nucleotides and their precursors, consistent with what we observed in HSC2 cells (Extended Data Fig. 6m; pool size data are available in Supplementary Table 5). Inhibition of SIRT6 as a driver of metabolic rewiring seems critical for tumour survival  Overlay image (right) of G6P and citrate with a correlation value of distribution between two metabolites. Scale bars indicate 300 µm. f, Representative immunofluorescence images of CD34, GLUT1 and keratin 10, and MALDI-MSI (GSH) from DMBA/TPA-treated skin tumours. Box plots show comparison of GSH abundance between different tumour subpopulations. Scale bars indicate 300 µm. g, Immunohistochemical analysis against MDA, a lipid peroxidation marker and SOX2, a functional TPC marker in the same tumour samples (serially sectioned). Scale bars indicate 100 µm. Statistics, sample sizes (n) and numbers of replications are presented in 'Statistics and reproducibility'. *P < 0.05, **P < 0.01, ***P < 0.001. and growth, since prolonged overexpression of SIRT6 caused apoptosis of these HSC2 tumour cells (Extended Data Fig. 7a), while knockdown of SIRT6 in the SCC13 human line enhanced proliferation in vitro (Extended Data Fig. 7b), a phenotype inhibited by DCA. We next assessed whether SIRT6 inhibition could impair human tumour growth in vivo. For this purpose, we took advantage of an in vivo model where the tumour cells are co-injected with tumour-associated fibroblasts and human primary keratinocytes. Strikingly, SIRT6 knockdown in SCC13 cells significantly increased tumour growth in vivo (Extended Data Fig. 7c-e).
Enhanced antioxidative response in glycolytic CD34 + tumour-propagating cells. We next sought to directly analyse metabolite levels of CD34 + TPCs from in vivo tumours by using MALDI Fourier transform ion cyclotron resonance mass spectrometry imaging (MSI), which provides direct evidence about metabolic characteristics of TPCs 20,21 . We used high spatial (25 μm) and spectral resolution to map the metabolic profiling directly on frozen tumour sections, preserving their structure and histology. A global view of the whole tumours using this cutting-edge technique already distinguished clear metabolic heterogeneity (Extended Data Fig. 8a,b); t-distributed stochastic neighbour embedding (t-SNE) analysis of all MALDI-MSI peaks, along with H&E staining of the same tissue sections, created 'metabolic images' with almost anatomical precision, both in its ability to separate tumours from adjacent tissues, and to depict metabolic heterogeneity within each tumour. To identify TPCs within these tumour samples, neighbouring sections were stained with CD34, as a TPC marker; GLUT1, as a glycolysis marker; and keratin 10, as a differentiated suprabasal cell marker (Fig. 3f). By using a non-linear transformation algorithm to co-register the two images (MALDI-MSI and immunofluorescence) from the sequential sections, we were able to appreciate the relative abundance of specific metabolites in different tumour subpopulations. Glucose-6-phosphate (G6P), one of the earliest glycolysis intermediates, and citrate, one of the mitochondrial tricarboxylic cycle intermediates, both of which yielded robust signals in MALDI-MSI, were used to determine relative glycolytic activity within a tumour ( Fig. 3e and Extended Data Fig. 8c-g). In general, signals from G6P and citrate showed weak correlation value (below 0.3). Importantly, mean abundance (or intensity) of G6P was significantly higher in CD34 + cells, while that of citrate was significantly lower in CD34 + cells than CD34 − cells (Extended Data Fig. 8c-g), strongly indicating increased glycolytic activity in CD34 + TPCs compared to non-TPCs. In this regard, reduced GSH levels were much higher in CD34 + cells (Fig. 3f), consistent with what we observed in highly glycolytic cells by in vitro metabolite profiling, further strengthening our finding that TPCs specifically increase glycolysis and PPP for GSH generation. To functionally assess whether glycolytic TPCs exhibit antioxidant properties, we checked the cellular redox state by using malonyldialdehyde (MDA), a marker of lipid peroxidation, along with SOX2, a TPC marker, in SIRT6-deficient skin tumours. Remarkably, the two markers exhibited a mutually exclusive staining pattern, further indicating that an important reason for the metabolic rewiring in TPCs is to protect against oxidative stress (Fig. 3g).
Single-cell RNA sequencing characterizes a subset of TPCs with higher glycolysis and antioxidant response, responsible for tumour progression. Our t-SNE analysis provided evidence for metabolic heterogeneity within TPCs ( Fig. 4a and Extended Data Fig. 8h). To determine in more detail whether CD34 + TPCs indeed exhibit metabolic heterogeneity, we took advantage of scRNA-seq. Dimensionality reduction analyses using both uniform manifold approximation and projection (UMAP) and t-SNE, and principal-component analysis separated TPCs of the WT2 sample from all the other TPCs, indicating that TPCs of the WT2 sample are qualitatively different from those of the other three tumour samples (Extended Data Fig. 9a-d). Based on tumour size and TPC enrichment, the WT2 sample was the smallest and the least aggressive tumour, suggesting that this tumour was likely regressing (as we observed with most of the WT tumours), and thus we excluded TPCs of WT2 from the analysis.
Dimensionality reduction analyses of TPCs from the other three samples showed that TPCs of each sample mixed together well regardless of sample identity (Extended Data Fig. 9e,f), generating four distinct clusters among TPCs with differentially expressed marker genes (Fig. 4b,c and Supplementary Table 6). Analysis of stemness and pro-differentiation programme (see gene lists in Supplementary Table 7) revealed that two distinct clusters, cluster I and IV, showed the highest stemness score and the lowest pro-differentiation score at comparable levels, implying that these are the most stem-like cells ( Fig. 4d and Extended Data Fig. 9g,h). Based on gene expression, cluster III might represent contamination of differentiating CD34 low cells during the prospective isolation by FACS, showing clearly less stemness markers and higher pro-differentiation scores with Sprr1a/Sprr1b expression (Fig. 4c,d and Extended Data Fig. 9g,h). We then compared cluster I and IV, separated by differences in their whole transcriptome while presenting similar stemness features. Remarkably, even within these two seemingly identical stem cell clusters (based on stem cell markers), cluster I exhibited significantly higher expression of genes involved in glycolysis, PPP, antioxidant response and GSH metabolism compared to that of cluster IV ( Fig. 4e; see gene lists in Supplementary  Table 7), defining, with unprecedented resolution, metabolically distinct subsets within bona fide TPCs. Notably, trajectory analysis in a pseudotime scale strongly suggests that although cells in cluster IV may represent the earliest progenitors in SCC, cells in cluster I, equally stem-like as in the cluster IV, could be the major contributors to aggressive SCC, in part due to their high glycolytic metabolism and ability to generate GSH (Extended Data Fig. 9i,j).

Redox regulation is critical for TPC enrichment and tumour progression.
To define whether increased GSH plays a functional role in TPCs, we administered an antioxidant, N-acetylcysteine (NAC), to animals treated with DCA in vivo. Strikingly, depleted CD34 + TPCs were completely rescued by exogenous antioxidant administration both in WT and Sirt6 cKO tumours ( Fig. 4f and Extended Data Fig. 10a). To confirm these results, we grew in vitro WT and shSIRT6 SCC13 cells in suspension, an approach that selects for cancer stem cell activity, as analysed by tumoursphere formation. Indeed, we saw a fourfold increase in tumourspheres in shSIRT6 cells, a phenotype abolished by DCA treatment, and fully rescued following addition of NAC ( Fig. 4g and Extended Data Fig. 10b).
Last, we took advantage of a dataset from a recent single-cell transcriptomic analyses of ten different human HNSCC samples 22 , and assessed glycolysis and antioxidant response gene expression. Highly glycolytic tumours (as defined by a 'glycolysis score'; Methods) exhibited a robust co-expression signature of antioxidant response genes (as defined by an 'antioxidant gene score'; Fig. 4h and Extended Data Fig. 10c,d). This positive correlation seems more obvious in the classical subtype of HNSCC, but less so in the atypical or basal subtypes of HNSCC that barely express glycolytic genes (Extended Data Fig. 10c,d). Furthermore, at a single-cell transcriptome level within each classical subtype of HNSCC (MEEI20 and MEEI6), the relationship between glycolysis and antioxidant response gene expression showed a significant positive correlation (Fig. 4h). Analyses from these independent human RNA-seq data further strengthen our findings that increased glycolysis and thus enhanced antioxidant response may promote squamous cell carcinogenesis, especially through TPCs.  I  II III IV  I  II III Fig. 2c. g, The number of tumourspheres at day 10 in SCC13 cells in indicated conditions. Data indicate the mean ± s.e.m. h, Scatterplot analysis of glycolysis score and antioxidant gene signature score in single cells from two classical subtypes of HNSCC using linear regression. Pearson correlation coefficients are 2.7 × 10 −48 for MEEI20 and 6.2 × 10 −6 for MEEI6. Statistics, sample sizes (n) and numbers of replications are presented in 'Statistics and reproducibility'. *P < 0.05, **P < 0.01, ***P < 0.001.

Discussion
The Warburg effect is considered a late adaptation of rapidly proliferating advanced tumours, yet several descriptive studies indicated that such adaptation may occur early, with increased glycolysis and its branched metabolic pathways (for example, serine/glycine metabolism) shown to be important specifically in tumour-initiating cells in lung, brain and breast cancer models [23][24][25][26][27] . In this study, using both genetically engineered mouse models and human SCC lines in which glycolytic metabolism was manipulated, we provide direct evidence for a driving role of glycolysis in TPCs, a phenotype regulated at least in part by epigenetic mechanisms. As such, to our knowledge, our data represent the first comprehensive analysis that metabolic reprogramming to increase glycolysis is critical for nucleotide biosynthesis and antioxidative functions in TPCs, thus identifying TPCs as the cells of origin for the Warburg effect in SCC. Importantly, despite exhibiting a Warburg phenotype (increased glucose uptake and increased glycolytic flux towards the PPP and lactate production), they appeared to maintain robust mitochondrial respiration, something that has been appreciated in several other cancer types 28 .
In cutaneous SCC, quiescent, transforming growth factor-βresponding SCC stem cells have an increased expression of genes in GSH metabolism and redox balance, providing cisplatin resistance and tumour recurrence 18 . In our study, we discovered that a subpopulation of highly proliferating CD34 + TPCs showed an augmented antioxidant response and nucleotide biosynthesis via increasing glycolysis, defining a metabolic mechanism to explain the ability of these cells to drive tumourigenesis. Recent studies have shown that both quiescent and proliferative TPCs, regardless of their proliferation states, exhibit elevated gene expression of the antioxidant response genes, including GSH metabolism genes; of note, most observations were in vitro, lacking in vivo relevance and evidence from direct measurements of metabolites 19 . Our study defined a subset of TPCs as exhibiting increased GSH and provided strong in vivo data for the importance of increased glycolysis in generating GSH and defending against oxidative stress. We further discovered a unique 'metabolic heterogeneity' within TPCs, indicating that only a defined (previously unknown) subpopulation of CD34 + TPCs acquires metabolic adaptations that could drive tumourigenesis. Our studies add metabolism to the large list of heterogeneous traits that have emerged in recent years in cancer cells 29 . These results may lead to a new therapeutic opportunity to target this specific subpopulation of TPCs in SCC by modulating glucose metabolism and, in turn, provide new hope for patients with these aggressive tumours for which treatment remains a challenge.

Methods
Mice and chemical-induced skin carcinogenesis. Mice were housed in pathogen-free facilities. All experiments were conducted under protocol 2019N000111 approved by the Subcommittee on Research Animal Care at Massachusetts General Hospital (MGH). Mice were maintained as a highly pure C57BL/6 background (>96%). Unless indicated, all animals were maintained on a standard diet (Prolab Isopro RMH 3000, 0006972). Data presented include both male and female mice. The Sirt6 F/F conditional strain 13 was crossed with the K14-cre strain (Jackson Laboratory). For flow cytometry analysis, these animals were further crossed with ROSA26-LSL-EYFP (Jackson Laboratory). Two days after shaving of their back hair, 8-week-old mice were subjected to a one-time DMBA (Sigma, 200 nmol in acetone) treatment followed by TPA (Sigma; 20 nmol in acetone) three times a week for 14 weeks. After 14 weeks of DMBA/TPA treatment, some mice were kept without TPA treatment for at least 7 weeks to observe tumour progression. The appearance and the number of tumours were closely monitored twice a week. Any visible mass that was more than 1 mm in size and existed for more than 1 week was considered as a tumour for onset and counting numbers of tumours. Some of the DMBA/TPA-treated animals were administered with DCA (Acros Organics, 5 g l −1 ) and/or NAC (Sigma, 1 g l −1 ) in their drinking water. NAC containing water was changed every week due to stability of this drug in aqueous solution. Blood from tail veins was collected into EDTA-coated tubes. Plasma was separated by centrifugation (15,000 r.p.m. for 10 min at 4 °C) for further analysis.
For the analysis and prospective isolation of TPCs, we generated K14-cre + ;Sirt6 F/F ;ROSA26-LSL-YFP (Sirt6 cKO) animals, as well as Sirt6 WT (K14-cre + ;Sirt6 +/+ ;ROSA26-LSL-YFP) animals to specifically isolate YFP + epithelial cells by FACS. In addition, because we had difficulty obtaining sizeable WT tumours with the previous DMBA/TPA treatment protocol, we modified the protocol by treating mice with TPA for more than 24 weeks to obtain sufficient number of cells to perform subpopulation analysis.
Humane end point: for all tumour assays, animals were euthanized according to the Institutional Animal Care and Use Committee protocol (tumours reached 20 mm, ulcerated mass or loss of 15% body weight).
Human datasets. SIRT6 expression and copy number, and glycolytic gene expression data were obtained from the Oncomine 30 and the CCLE 31 . SIRT6 copy number and the corresponding survival data of each participant were obtained from the TCGA 32 . A Kaplan-Meier plot was made and the log-rank P value was calculated for the SIRT6 copy number in the TCGA HNSC samples. TCGA HNSC clinical information was downloaded from the TCGA data matrix access portal (http://cancergenome.nih.gov/) during December 2015 and copy number data processed by Genomic Identification of Significant Targets in Cancer 2.0 (GISTIC2) for TCGA HNSC were downloaded from the Broad Institute (http:// gdac.broadinstitute.org/) during April 2015. Follow-up clinical data files were merged with the original clinical data file to ensure that the most up-to-date patient follow-up information was used for survival analysis. The Kaplan-Meier survival plot and the log-rank P value were generated using R using the 522 TCGA HNSC primary tumour samples of which both clinical data with a death event or at least 6 months of follow-up and GISTIC2 copy number data were available. Kaplan-Meier plots used overall survival with death from any cause as the end point, and individuals still alive at last follow-up were censored at last follow-up time. Samples were split into five groups based on the GISTIC2 data for SIRT6 in the all_thresholded.by_genes.txt results file (a gene-level table of discrete amplification and deletion indicators at for all samples). Within the GISTIC2 table, values of 0 represent no amplification or deletion above the threshold (0.1), while positive numbers represent amplifications, and negative numbers represent deletions (1 denotes amplification above the amplification threshold, 2 denotes amplifications larger than any arm-level amplifications observed for the sample, −1 represents deletion beyond the threshold, and −2 represents deletions greater than the minimum arm-level deletion observed for the sample).
Human tumour sample analysis. Primary HNSCC tumours and newborn foreskins (used as normal tissue) analysed in this study were collected following institutional review board (IRB) approval. The collection and use of discarded, de-identified tissue was reviewed and approved by the Dana-Farber/Harvard Cancer Center (DF/HCC) IRB (protocol no. 03-204). Briefly, tissues were fixed in formalin, followed by ethanol and paraffin embedding and then serial sectioning onto slides for immunohistochemistry.
Real-time quantitative PCR with reverse transcription analysis. Total RNA was extracted with the TriPure isolation reagent (Roche) as described by the manufacturer. For RNA isolation from mouse skin or tumour samples, an additional RNA clean-up procedure was performed with the RNeasy Protect Mini kit (Qiagen). Complementary DNA (cDNA) synthesis and real-time PCR were performed as previously described 14 . In brief, isolated RNA was reverse transcribed by using the QuantiTect Reverse Transcription kit (Qiagen). Real-time PCR was performed using SYBR green master mix (Roche) with the final volume of 12.5 μl per reaction in LightCycler 480 detection system (Roche). Data were presented as relative mRNA levels normalized to the β-actin expression level in each sample. The primer sequences are listed in Supplementary Table 8.

Mouse tumour cell isolation and fluorescence-activated cell sorting.
Mice bearing skin tumours were euthanized and the tumours were collected on ice. Each tumour was cut into small pieces and incubated with 0.5% trypsin (diluted in keratinocyte serum-free medium; Gibco) on a horizontal shaker at 37 °C for 1.5 h. With an 18-gauge syringe, digested tumour cells were physically isolated into a single-cell suspension. The trypsin was inactivated by adding chelexed FBS. After serial filtering with 70-μm and 40-μm strainers (BD Biosciences), tumour cells were centrifuged at 1,200 r.p.m. at 4 °C for 10 min. Cell pellets were resuspended with PBS containing 4% chelexed FBS and then transferred into FACS tubes with a 40-μm filter. The following fluorophore-conjugated antibodies and their dilutions were used: anti-CD34-BV421 (BD Biosciences, 562608; 1:50) and anti-CD49f-PE (eBiosciences, 12-0495-81; 1:200) and anti-GLUT1-A647 (Abcam, ab195020; 1:100). PI (Sigma, P4864; 1:1,000) or Zombie NIR fixable viability dye (BioLegend, 423105; 1:100) were used to negatively select live cells. Proper isotype controls, single-colour controls and fluorescence-minus-one controls were used in every experiment to set up optimal compensation and gates. Cells were analysed and sorted using a FACSAria II (BD Biosciences). Obtained data were analysed by FlowJo. An exemplary gating strategy is described in Extended Data Fig. 4a.
RNA preparation and RNA sequencing. Following FACS isolation of three different tumour subpopulations (α6 high /CD34 + , α6 high /CD34 − and α6 low / CD34 low ) from at least two independent skin tumours obtained from Sirt6 WT and Sirt6-deleted animals, sorted tumour cells were directly collected into RNA isolation buffer provided in the kit (Clontech, 740902.50). RNA isolation was conducted according to the manufacturer's instructions. Library construction was performed using the SMART-Seq v4 Ultra-Low Input RNA kit to produce cDNA (Clontech, 634888). The total RNA input amount for this kit was 10 ng. Eight cycles of PCR were performed for PCR amplification. After cDNA construction, the samples were validated using an Agilent 2100 Bioanalyzer and Agilent's High-Sensitivity DNA kit. Before generating the final library for Illumina sequencing, the Covaris AFA system was used to shear cDNA resulting in a 200-500 bp size range. Sheared libraries were validated using Agilent 2100 Bioanalyzer and Agilent's High-Sensitivity DNA kit. Quantification was completed by using a Qubit 4 fluorometer (Invitrogen) using the Qubit RNA HS Assay kit. Generation of the final library was completed by using a Low Input Library Prep Kit v2 (Clontech, 634899). The total cDNA input amount was 10 ng. Seven cycles of PCR were performed for PCR amplification. After library construction, the samples were validated using the 2200 TapeStation System and High-Sensitivity D1000 ScreenTape kit. Libraries were quantified using the Library Quantification kit (Kapa Biosystems, KK4828) and the Bio-Rad CFX96 instrument. Each lane of sequencing was pooled in 6-plex (six samples per lane) with unique barcodes. Pooled libraries were also quantified using the Kapa Biosystems Library Quantification kit (KK4828) and the Bio-Rad CFX96 instrument. These pools were then denatured to 16 pM with 1% PhiX and sequenced on the Illumina HiSeq 2000 instrument, producing approximately 30 million paired-end 50-bp reads per sample.
RNA-sequencing analysis. STAR aligner 34 was used to map sequencing reads to the mouse reference transcriptome (mm9 assembly). Read counts over transcripts were calculated using HTSeq (v.0.6.0) 35 based on a current Ensembl annotation file for NCBI37/mm9 assembly. Differential expression analysis was performed using EdgeR 36 , and genes were classified as differentially expressed based on cut-offs of at least a twofold change. Analysis of enriched functional categories among detected genes was performed using DAVID 37 .
Cell lines and cell culture. SCC13 cells were grown in keratinocyte serum-free medium (Gibco) supplemented with epidermal growth factor and bovine pituitary extract according to the manufacturer's instructions. HSC2 cells were grown in DMEM/F-12 medium with 10% Tet system approved FBS (Clontech) and 1% penicillin-streptomycin (both 100 U ml −1 ; Gibco). Human primary keratinocytes were obtained from CellnTec and were grown in CnT-PR medium (CellnTec) supplemented with 1% penicillin-streptomycin (both 100 U ml −1 ; Gibco). Human tumour-associated fibroblasts were grown in DMEM supplemented with 10% FBS (Sigma), insulin-transferrin-selenium reagent (Gibco) and 1% penicillinstreptomycin (both 100 U ml −1 ; Gibco). All cells were cultured and maintained at 37 °C under 5% CO 2 . Human primary keratinocytes (CellIn Tec) were passaged with accutase (Gibco) and were used within four passages. All other cell lines were passaged by trypsinization.
Constructs and viral infection. Human pTripZ-shSIRT6 construct (Dharmacon, RHS4740), negative control shRNA vector and pMSCV-luc-PGK-Neo-IRES-eGFP construct were kind gifts. pLVX-Tet-On was obtained from Clontech. Human pRetro-SIRT6 WT and pRetro-SIRT6-H133Y were previously described 17 . Viral particles containing the above-mentioned constructs were generated using either lentiviral (pCMV-d8.9) or retroviral (pCL-ECO) packaging plasmids with pCMV-VSV-G (Addgene) in 293T cells. Virus-containing supernatant was filtered through a 0.45-μm filter and added into target cell lines with 8 μg ml −1 polybrene. For infection of SCC13 cells, filtered virus-containing supernatant was ultra-centrifuged at 20,000g at 4 °C for 2 h to concentrate it into a very small volume (~200 μl), and about 5 μl of virus concentrate was used for infection of SCC13 cells. For efficient infection, six-well plates with SCC13 cells with added virus concentrate and polybrene were centrifuged at 2,250 r.p.m. at 32 °C for 1 h, and virus-containing medium was immediately replaced by regular keratinocyte serum-free medium. The next day, cells were selected in 2 μg ml −1 of puromycin or 1.4 mg ml −1 of neomycin for SCC13 cells, or in 1.5 μg ml −1 puromycin or 0.5 mg ml −1 of neomycin for HSC2 cells for at least a week, and the pooled populations were used for various experiments. SCC13 cells with doxycycline-inducible pTripZ constructs were treated with 1 μg ml −1 doxycycline for at least 3 d. HSC2 cells with doxycycline-inducible pRetro constructs were treated with 100 ng ml −1 doxycycline for 26 h, unless otherwise indicated.
Glucose uptake assay. In SCC13 cells, cells were plated in duplicate on six-well plates (2.5 × 10 5 cells per well) in culture medium 1 d before the experiment. Media containing 2-NBDG (Invitrogen; 100 μM) were added for 2 h. Fluorescence was measured using a FACSAria II. Obtained data were analysed by FlowJo. After proper compensation, geometric mean value was normalized by its negative control (without 2-NBDG) of each group.
Glycolytic capacity. Cells were plated into XFe96 cell culture microplates (Seahorse Bioscience) 1 d before the experiment. Medium was replaced in the Seahorse microplates with assay medium supplemented with 2 mM l-glutamine (Gibco; pH 7.35 ± 0.05) for glycolysis stress tests. The plate was incubated in a CO 2 -free incubator for 1 h at 37 °C. For the glycolysis tests, 10 mM of glucose, 2 μM of oligomycin and 100 mM of deoxyglucose were sequentially injected. Extracellular acidification rate was measured in every well based on the instrument's protocol. Experiments were run using an XFe96 analyser, and raw data were normalized by cell number calculated with a Cyquant cell proliferation assay kit (Thermo Fisher Scientific).
Chromatin immunoprecipitation. Chromatin immunoprecipitation followed by RT-qPCR was performed as previously described 38 with H3K9Ac (Millipore, 07-352). In brief, SCC13 cells were crosslinked with 1% formaldehyde/PBS for 15 min at RT, followed by quenching with 0.125 M glycine. After washing twice with PBS, cells were collected in RIPA buffer and then sonicated to generate DNA fragments of approximately 500 bp in size. About 200 μg of pre-cleared protein extract was used for immunoprecipitation overnight (>12 h) at 4 °C using protein A/G agarose beads (Santa Cruz, sc2003). Washed samples were eluted by incubation at 65 °C for 10 min with 1% SDS, and crosslinking was reversed by incubation at 65 °C for 6 h with 200 mM NaCl. DNA was purified using the QIAquick spin kit (Qiagen) and assessed by real-time qPCR using the LightCycler 480 system (Roche). The primers' sequences are listed in Supplementary Table 9.
Isotope tracing experiment. HSC2 cells were maintained with doxycycline for 26 h, at which point intracellular metabolites were collected after adding labelled medium at different time points. Glucose-free DMEM/F-12 medium (US Biological) was supplemented with 10% dialysed FBS, 17.5 mM of [U-13 C] glucose or [1,2-13 C]glucose (Cambridge Isotope Laboratories), and 15 mM HEPES. SCC13 cells were maintained with doxycycline for 3 d, at which point intracellular metabolites were collected after adding labelled medium at different time points. Keratinocyte serum-free medium (Gibco) was supplemented with 5.8 mM of [U-13 C]glucose (Cambridge Isotope Laboratories). At different time points after replenishing labelled medium, medium from biological triplicates (in six-well plates) was fully aspirated. Each well was washed with 0.9% ice-cold NaCl twice and 1 ml of 80% (vol/vol) methanol was added at dry-ice temperature. After vigorous vortexing, insoluble material in lysates was centrifuged at 16,000g at 4 °C for 10 min. The supernatant was transferred and the solvent was evaporated using a SpeedVac. Samples were stored at −80 °C until analysis.
Metabolite profiling by liquid chromatography-mass spectrometry. Metabolite profiling and isotope tracing liquid chromatography-mass spectrometry (LC-MS) analyses were conducted on a QExactive benchtop orbitrap mass spectrometer equipped with an Ion Max source and a HESI II probe, which was coupled to a Dionex UltiMate 3000 HPLC system (Thermo Fisher Scientific). External mass calibration was performed using the standard calibration mixture every 7 d. Typically, samples were reconstituted in 100 μl water and 2 μl was injected onto a SeQuant ZIC-pHILIC 150 × 2.1 mm analytical column equipped with a 2.1 × 20 mm guard column (both 5-mm particle size; EMD Millipore). Buffer A was 20 mM ammonium carbonate and 0.1% ammonium hydroxide; Buffer B was acetonitrile. The column oven and autosampler tray were held at 25 °C and 4 °C, respectively. The chromatographic gradient was run at a flow rate of 0.150 ml min −1 as follows: 0-20 min: linear gradient from 80-20% B; 20-20.5 min: linear gradient from 20-80% B; 20.5-28 min: hold at 80% B. The mass spectrometer was operated in full-scan, polarity-switching mode, with the spray voltage set to 3.0 kV, the heated capillary held at 275 °C, and the HESI probe held at 350 °C. The sheath gas flow was set to 40 units, the auxiliary gas flow was set to 15 units and the sweep gas flow was set to 1 unit. MS data acquisition was performed in a range of m/z = 70-1,000, with the resolution set at 70,000, the AGC target at 1 × 10 6 , and the maximum injection time at 20 ms. An additional scan (m/z 220-700) in negative mode only was included to enhance detection of nucleotides. Relative quantification of polar metabolites was performed with XCalibur QuanBrowser 2.2 (Thermo Fisher Scientific) using a mass tolerance of 5 ppm and referencing an in-house library of chemical standards. For stable-isotope tracing analyses, data were corrected for natural abundance 39 . Metabolite pool sizes of the above-mentioned metabolites are described in Supplementary Table 5.
Lactate measurement by gas chromatography-mass spectrometry. Ice-cold methanol was added into 5 µl of plasma from tail vein blood, vigorously vortexed at 4 °C for 10 min, followed by centrifugation. The supernatant was transferred and the solvent was evaporated using a SpeedVac. Samples were stored at −80 °C until analysis. Derivatization and measurement by the gas chromatography-mass spectrometry was done as previously described 40 . In brief, polar metabolites were derivatized with 20 mg ml −1 methoxyamine in pyridine for 90 min at 37 °C and subsequently with N-(tert-butyldimethylsilyl)-N-methyl-trifluorosilane and 1% tert-butyldimethylchlorosilane for 60 min at 60 °C. Metabolite levels were then measured with a 5977B GC system (Agilent Technologies). The raw ion chromatograms were extracted to determine metabolite levels using a custom MATLAB M-file 41 .
Reactive oxygen species measurement. Cells were plated in duplicate on 12-well plates (1 × 10 5 cells per well) in culture medium. After 3 d in the presence or absence of doxycycline (1 μg ml −1 ), cells were trypsinized, centrifuged and resuspended in the medium. Cells were incubated with CellROX Deep Red (Invitrogen) for 30 min at 37 °C. Fluorescence was measured using an LSR II (BD Biosciences). Obtained data were analysed by FlowJo.
Proliferation assay. SCC13 cells were pretreated with doxycycline (1 μg ml −1 ) for at least 3 d before plating cells. Cells were plated in triplicate on six-well plates (1 × 10 4 cells per well) in culture medium with doxycycline in the presence or absence of DCA. Adherent cells were trypsinized and counted by trypan-blue exclusion at 2, 3 and 5 d later.
Apoptosis assay. Cells were plated in duplicate on 12-well plates (2.5 × 10 4 cells per well) in culture medium. After 4 d in the presence of doxycycline (100 ng ml −1 ), all the floating and adherent cells were collected and stained with Annexin V-FITC (BD Biosciences) and PI (Sigma) to analyse cell death in Accuri (BD Biosciences). Obtained data were analysed by FlowJo.
Skin xenotransplantation assay and bioluminescence imaging. SCC13 shCtrl or shSIRT6 cells were stably transduced with retrovirus containing pMSCV-luc-PGK-Neo-IRES-eGFP. After neomycin selection (1.4 mg ml −1 ), infected cells were further sorted with green fluorescent protein (GFP) using a FACSAria II (BD Biosciences) and sorted GFP high cells were cultured before the experiment. A total of 4,000 cells of SCC13 shCtrl or shSIRT6 cells, 1,000 cells of human primary keratinocytes and 500 cells of tumour-associated fibroblasts were prepared and mixed immediately before injection. Mixed cells were injected after inserting a silicone chamber (Renner) in the back skin of NSG mice (Jackson Laboratory) under avertin anaesthesia. Doxycycline (200 μg ml −1 ) was administered in the drinking water and was replaced every week due to its light sensitivity. After 8 d, the silicone chamber was removed from the back skin of the mice under avertin anaesthesia, and several sutures were made to aid the wound healing process. All surgical procedures were performed in the aseptic hood of the pathogen-free facility. Every mouse was singly housed due to the small open-wounded area after silicone chamber removal. From day 23, all mice were subject to bioluminescence imaging once a week. Under isoflurane anaesthesia, 300 μl of d-luciferin (15 mg ml −1 ; RR LABS) was injected intraperitoneally, and the mice were imaged every 5 min after injection with a 0.5 s to 60 s exposure time with a binning of 8 and 4 on an Ami X imaging system (Spectral Instruments Imaging) until the total flux and maximal radiance peaked. Total flux (photons/s) and maximal radiance (photons/s/cm 2 /sr) were measured by Spectral Ami X (Spectral Instruments Imaging). A region of interest (ROI) was drawn around the tumour region for each mouse as well as a background ROI outside of the mice, which was subtracted from the total flux and maximum radiance for each mouse.
Humane end point: for all tumour assays, animals were euthanized according to the Institutional Animal Care and Use Committee protocol (tumours reached 20 mm, ulcerated mass or loss of 15% body weight).
Matrix-assisted laser desorption/ionization mass spectrometry in skin tumours. DMBA/TPA-treated skin tumours and adjacent skin samples were collected, snap frozen in liquid nitrogen and stored at −80 °C until analysis. Samples were embedded in a solution of 2% methylcellulose for sectioning in the Microm HM550 cryostat (Thermo Fisher Scientific) at a thickness of 10 µm. The cryo-sections were mounted in indium tin oxide slides pre-coated with 90% CHCA matrix solution (5 mg ml −1 ) in acetonitrile for MALDI-MSI, and consecutive sections were mounted in regular glass slides for immunostaining. The indium tin oxide slides were then coated with 90% 9-aminoacridine (5 mg ml −1 in acetonitrile) using an automated matrix sprayer (TM-sprayer; HTX imaging) under the following conditions: four passes, 0.17 ml min −1 flow rate, 75 °C, 10 p.s.i. nitrogen pressure and 1,000 mm min −1 speed. A 9.4T SolariX XR FT-ICR mass spectrometer (Bruker Daltonics) with a MALDI source was used to acquire spectra in negative-ion mode, with a raster pixel size of 25 µm. The mass range was m/z 60-1,500 with continuous accumulation of selected ions set to m/z 100-520. A tuning mix was used for external calibration in the selected mass range. FlexImaging 5.0 was used for acquiring the images, while SCiLS 2019c was used for data processing.

Multinodal image processing and integration.
1. H&E and immunofluorescence: H&E and immunofluorescence images were pre-processed 42 to exclude the background noise. These two imaging modalities exhibit local deformations (such as compression and tearing) due to the manual nature of the tissue sectioning process and the underlying different staining procedures. The iterative parametric-based image registration framework was used to non-linearly warp the H&E image to be spatially aligned with immunofluorescence image 43 . The non-linear image registration process was first initialized by Affine transformation to capture global deformations (scaling, rotation, translation, and shearing) and map the images to the same coordinate space. The non-linear transformation model of cubic B-spline was used to model the local deformations. The cost function of mutual information has been found efficient as a similarity metric for multimodal registration 44 and it was optimized by the adaptive stochastic gradient descent algorithm 45 . For robustness and faster convergence, the image registration was implemented using a multi-resolution strategy that used eight levels of Gaussian smoothing 46 . The registration algorithm was implemented using the publicly available toolbox elastix 43 , eventually yielding an optimized non-linear transformation model, T μ 1 I . 2. MALDI-MSI and H&E: the complex nature of MSI data poses challenges that hinder direct co-registration with histology. This complexity is described in terms of high dimensionality and the lack of established spatial correspondences between biochemical and anatomical images. We adopted the t-SNE-based non-linear image registration methodology developed by Abdelmoula et al. 42 . Briefly, the t-SNE computes pairwise similarities of the high-dimensional data points (that is, spectra) and non-linearly maps data into a lower dimensional embedding of three dimensions. The t-SNE non-linearity preserves the local similarity of the higher dimensional data points in the embedding space as such similar spectra are projected closely to each other, whereas dissimilar ones are projected further away. The embedding features were spatially mapped to form a t-SNE image that reveals structures in which edges demarcate molecularly distinct regions. The t-SNE image was non-linearly aligned to histology using the elastix toolbox 43 . The registration was initialized with Affine transformation to capture linear deformation and then followed by a cubic B-spline transformation to capture non-linear deformations. The cost function of mutual information and the multi-resolution strategy of four levels of Gaussian smoothing were used. The resultant transformation model, T μ 2 ; I was applied to spatially align a given m/z image to the histological image, and then was aligned to immunofluorescence image using the previously computed model T μ 1 I . 3. ROI and multimodal correlation: we chose ROIs based on matching histology between H&E and immunofluorescence slides because they were not the exact same slides and also excluded edge and folded areas due to strong auto-fluorescence signals. In selected ROIs, signal intensities of MALDI-MSI images were quantified based on the cellular markers defined by immunofluorescence images, using in-house MATLAB codes.
Single-cell RNA-sequencing data generation and processing. FASTQ files were aligned to the NCBI Genome Reference Consortium Mouse Build 38 (mm10) using STAR 34 . Expression levels were computed using the RSEM tool and quantified as both raw counts and transcripts per million (TPM) 48 . For each cell, we used four quality-control metrics. We excluded (1) cells expressing less than 2,500 genes with three or more counts, (2) cells with less than 500,000 RNA molecules detected, (3) cells with mitochondrial gene counts exceeding 10% of total gene expression 49 and (4) cells with an average expression of housekeeping genes 47 (log 2 (TPM + 1) < 6).
Dimensionality reduction. The t-SNE and UMAP methods were used for dimensionality reduction within the R package Seurat 50 . The dimensionality of the dataset was determined by principal-component analysis, and the first five principal components (PCs) were used as inputs to the dimensionality reduction.
Unsupervised clustering of tumour cells. TPCs were clustered using Seurat (v. 3.1.1) 50 . The same PCs as in the dimensionality reduction were used as input to the clustering analysis. The resolution parameter was set at 0.5 to identify the greatest number of clusters with at least 20 significant marker genes. To be considered a marker gene, the gene must be enriched in the given cluster with a Wilcoxon rank-sum test Bonferroni-adjusted P value of less than 0.05 and a minimum log fold change of 0.25 when comparing the cells of the cluster to all other cells. Statistical analysis was performed in R (v.3.6.1).
Pathway score analyses of tumour-propagating cells from single-cell RNA sequencing. For each single cell, we assigned a score for each of several programmes based on the average expression of a selected gene set of interest minus the average of a control gene set, based on the work of Puram et al. 22 . For glycolysis and GSH metabolism, we began with a list of relevant genes and calculated the Pearson correlation coefficient between all pairs of genes in the list. We sequentially removed genes with the strongest negative correlations to other genes until all genes were positively correlated in the list. Then, we removed genes with the weakest positive correlations until all pairs of genes had a Pearson correlation coefficient of at least 0.2. For the PPP and pro-differentiation, we followed a similar process of elimination but stopped eliminating genes once all correlations were positive, as the correlations were weaker. For stemness, we used a list of genes associated with stemness but did not eliminate any genes as the genes on the list were largely uncorrelated. For antioxidant response, we ranked the entire list of analysed genes by their correlation with the NRF2 transcription factor controlling the antioxidant response. We removed unannotated genes to only include genes with known functions, and selected the top 30 positively correlated genes to make the gene list of interest. As described by Puram et al., the gene expression levels of a given programme may be confounded by the overall complexity of a given cell, since cells with high expression complexity would be expected to score well for any given pathway 22 . To control for this, the cell score for a given pathway was defined as SC j (i) = average (Er(G j ,i)) -average (Er(G j cont ,i)), where G j is a given gene set, SC j (i) is the score of each cell i, and Er is average relative expression. The cell score is then the score of the given pathway in that cell minus the score of a control gene set. The control gene set was selected by dividing all analysed genes into 25 bins of equal size according to their overall expression, and for each gene in G j , selecting 100 genes from the same expression bin.

Trajectory analysis of tumour-propagating cells.
To analyse the differentiation trajectory of the TPCs based on the scRNA-seq gene expression data, we used Monocle (v.2.13.0) 51 . Monocle identified a set of 14,321 DEGs between the four clusters established in Seurat. These genes were ordered by ascending q value, and the top 1,000 genes were used as input to Monocle's Reversed Graph Embedding algorithm. Branched expression analysis modelling was used to identify genes with branch-dependent expression, and a q-value cut-off of 1 × 10 −15 was applied.
Single-cell transcriptomic analysis of human HNSCC samples. We restricted the analysis to malignant cells, as defined previously 22 , and scored each cell for each programme, based on the average expression of a set of predefined genes, minus the average of a control gene set (see work by Puram et al. 22 for further description of how the control scores are selected). For glycolysis, we initially examined the gene sets PDK1, G6PD, PGD, PFKM, LDHB and LDHA. Since five of those genes were all positively correlated with one another (across single cancer cells), while the sixth (LDHA) was not correlated and it was removed from the list; the remaining five genes were used to define final scores. For antioxidant gene score, we first calculated the correlation of each gene with NRF2 across all cancer cells. We then selected the top 30 genes (including NRF2) to define antioxidant gene score (NRF2) scores. As expected, this list included many known downstream targets of NRF2 such as GPX2 and genes associated with, for example, detoxification and antioxidants.
Statistics and reproducibility. Figure 1: b, Statistical analysis was done by log-rank test (P = 0.0176). WT, n = 9; cKO, n = 14. c, Student's t-test was performed (two-sided). WT, n = 34 biologically independent tumours; cKO, n = 36 biologically independent tumours. d, Fisher's exact test was performed for statistical analysis (P < 0.0001, two-sided). WT, n = 7; cKO, n = 11. e, Upper: Student's t-test was performed (P < 0.0001, two-sided). At least four different measurements in two different microscopic images (×40) from two biologically different samples were analysed for normal skin samples. At least four different measurements in three different microscopic images (×40) from three biologically different samples were analysed for tumour samples. The number of measurements were as follows; WT normal, n = 16; cKO normal, n = 16; WT tumour, n = 34; cKO tumour, n = 45. Lower: immunostaining against PCNA was performed six times with similar results. g, Student's t-test was performed (two-sided). WT, n = 34 biologically independent tumours; cKO, n = 36 biologically independent tumours; WT with DCA, n = 28 biologically independent tumours; cKO with DCA, n = 43 biologically independent tumours. h, Fisher's exact test was performed for statistical analysis (P < 0.0001, two-sided). WT, n = 7; cKO, n = 11; WT with DCA, n = 4; cKO with DCA, n = 4. i, Immunostaining against GLUT1 and phospho-PDH was performed ten and three times, respectively, with similar results. Figure 2: a, Immunofluorescence staining was performed five times with similar results. e,f, Data are from at least two biological replicates (n = 3 for WT and n = 2 for Sirt6 cKO), presented as the mean. Figure 3: a, Data are from three biological replicates. b, Data are from three biological replicates. c, Student's t-test was performed (two-sided). Data are from three biological replicates. d, Student's t-test was performed (two-sided). Data are from at least two biological replicates (n = 2 for shCtrl and n = 3 for shSIRT6). e, Data are from two biologically independent tumour samples, consisting of more than 100,000-pixel data points per sample. f, Immunofluorescence staining was performed twice with similar results. The non-parametric Wilcoxon rank-sum test (two-sided and 95% significance level) was used after checking normality distribution using the Kolmogorov-Smirnov test. Details of the box plots are listed in the source data. g, Immunostaining was performed three times with similar results. Figure 4: a, Data are always reproducible every time the code was run. The key point for this stability was that the random seed point in the t-SNE algorithm was set to zero, which maintained reproducibility. The t-SNE analysis was done on MATLAB 2018a, which was installed on a workstation operating with Windows 10. d, Right: analysis of variance (ANOVA) test was performed. F values were as follows: stemness, F = 14.47; pro-differentiation, F = 95.23. e, Pairwise comparisons using a t-test were used to calculate P values (P-value adjustment was done by the Benjamini and Hochberg method). g, Data are from three biological replicates. A paired student t-test was performed (two-sided).
Extended Data Fig. 1a, Log-rank test was performed. b and c, left panels, Student t-tests were performed. Details of the box plots are listed in the Source Data file 2. d, Immunostaining was performed two times with similar results.
Extended Data Fig. 2b, Immunostaining was performed ten times for GLUT1 and six times for PCNA with similar results. c, Each dot represents one biologically independent tumour sample. (n = 1, 2, 6, 5, respectively, from left to right in each graph) Student's t-tests were performed (two-sided). d, Immunostaining against GLUT1, p-PDH, and MPC1 was performed ten, three, and two times, respectively, with similar results. e, Definitions of each box plot are listed in Source Data file 3. f, Immunostaining was performed two times with similar results.
Extended Data Fig. 3a,b, Immunostaining was performed three times with similar results. Quantification is done by at least three independent 20x images with hundreds of positive cells. c, Immunostaining was performed two times with similar results.
Extended Data Fig. 4b, Each dot represents one biologically independent tumour sample. (n = 2, 6, 3, 6, 5, 3 respectively, from left to right in each graph) Student's t-tests were performed (two-sided). c, Each dot represents one biologically independent tumour sample. (n = 6, 5, 3 respectively, from left to right in each graph) Student's t-tests were performed (two-sided). g and h, data are from at least two biological replicates (n = 3 for WT and n = 2 for Sirt6 cKO), presented by mean.
Extended Data Fig. 5b, Each dot represents one biologically independent sample, presented by the mean and s.d. (n = 9 for S6HY and n = 16 for S6WT). Student's t-tests were performed (two-sided). c and e, data are from three biological replicates. d, data are from three biological replicates.
Extended Data Fig. 6b, Data are from two independent experiments with two experimental replicates. Student's t-tests were performed (two-sided). c, Immunostaining was performed three times with similar results. d, Immunostaining was performed three times with similar results. e, data are from three independent experiments with two experimental replicates. Student's t-tests were performed (two-sided). f-h, Data are from four biological replicates. Student's t-tests were performed (two-sided). i, Each dot represents one biologically independent sample, presented by the mean and s.d. (n = 14 for shCtrl and n = 14 for shSIRT6). Student's t-tests were performed (two-sided). j, data are from three biological replicates. k, data are from two independent experiments with two experimental replicates. Student's t-tests were performed (two-sided). l and m, data are from at least two biological replicates (n = 2 for shCtrl and n = 3 for shSIRT6).
Extended Data Fig. 7a, data are from three independent experiments (n = 10 each sample). Student's t-tests were performed (two-sided). b, data are from three independent experiments. Two-way ANOVA test was performed. d, Two-way ANOVA test was performed. e, Immunostaining was performed two times with similar results.
Extended Data Fig. 8a, data are always reproducible every time the code was run. The key point for this stability is the random seed point in the t-SNE algorithm is set to zero and that maintained reproducibility. The t-SNE analysis was done on MATLAB 2018a that was installed on a workstation operating with Windows 10. c and e, data are from two biologically independent tumour samples, consisting of more than hundred thousands of pixel datapoints per sample. d and f, The non-parametric Wilcoxon rank-sum test (two-sided and 95% significance level) was used after checking normality distribution using Kolmogorov-Smirnov test. Details of the box plots are listed in the Source Data file 4. g, data are from two biologically independent tumour samples, consisting of more than hundred thousands of pixel datapoints per sample. h, data are always reproducible every time the code was run. The key point for this stability is the random seed point in the t-SNE algorithm is set to zero and that maintained reproducibility. The t-SNE analysis was done on MATLAB 2018a that was installed on a workstation operating with Windows 10.
Extended Data Fig. 10a, Each dot represents one biologically independent tumour sample. (n = 5, 3, 10, 7, respectively, from left to right in each graph) Student's t-tests were performed (two-sided). b, Bright field imaging was performed two times with similar results.
Supplementary Table 6: To be considered a marker gene, the gene must be enriched in the given cluster with a Wilcoxon rank-sum test Bonferroni-adjusted P value of less than 0.05 and a minimum log fold change of 0.25 when comparing the cells of the cluster to all other cells. Statistical analysis was performed in R (v.3.6.1).
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The RNA-seq data of the tumour subpopulations from mouse cutaneous tumours of Sirt6 WT or cKO animals have been submitted to the Gene Expression Omnibus (GEO) database under accession number GSE115953 (related to Fig. 2d-f and Extended Data Fig. 4e-h). The scRNA-seq data of TPCs from mouse cutaneous tumours of Sirt6 WT or cKO animals have been submitted to the GEO under accession GSE147031 (related to Fig. 4b-e and Extended Data Fig. 9a-j). There is no restriction on data availability. Human HNSCC scRNA-seq data from Puram et al. 22 is available under accession GSE103322 (related to Fig. 4h and Extended Data Fig. 10c,d). Oncomine and TCGA datasets (related to Extended Data Figs. 1a-c and 2e) are available in cBioportal (https://www.cbioportal.org/). CCLE data (related to Extended Data Fig. 1e) are available in https://portals. broadinstitute.org/ccle/. Source data are provided with this paper.

Code availability
In-house codes were previously used in the published works. Appropriate references to the original works are provided.