Single-cell transcriptome profiling of the stepwise progression of head and neck cancer

Next, we assigned the epithelial cells in CA and LN tissues as malignant cells (n = 5,113). To assess the transcriptional heterogeneity, the malignant cells were sub-classified into six clusters (CC0 to CC5, Fig. 3a, top, see Supplementary Methods). We observed that each pair-matched primary and metastatic tumor (CA and LN) was clustered together as demonstrated previously21 (Fig. 3a, bottom). In addition, we observed that the cell composition of the malignant cell clusters of the metastasized tumors was very similar to that of the pair-matched primary tumors, revealing that each pair of the primary and metastatic tumors was clustered together in the hierarchical cluster analysis (Fig. 3b). This result suggests that the transcriptional diversity of the primary tumors is conserved in metastatic tumors.

Next, we performed deconvolution analyses to evaluate the cell type proportions in bulk RNA-seq data of TCGA-HNSCC (n = 500, Supplementary Methods). Unlike scRNA-seq data, deconvolution analysis revealed that the malignant cells were the most abundant cell type (79.17%) followed by fibroblasts (10.07%) (Supplementary Fig. 4a). Among the malignant cells, CC0 and CC1 were the most abundant malignant cell clusters (Fig. 3c). Similar results on malignant cell composition were observed in independent datasets of GSE41613 (n = 97)22, GSE42743 (n = 74)22, and GSE65858 (n = 270)23, demonstrating the robustness of our result (Supplementary Fig. 4b). Interestingly, we found that the malignant cell clusters were closely located according to the HPV infection status (Fig. 3d, top). CC0 and CC4 cells were mostly HPV-positive tumors, whereas the other clusters were HPV-negative tumors (Fig. 3d, bottom). This finding indicates that our malignant cell clusters harbor distinct transcriptomic programs related to HPV infection. In addition, the tumor subsites (e.g., oropharynx and oral cavity) have been shown to associate with HPV infection status;24 however, we could not observe significant associations between the malignant cell clusters and tumor subsites (Fig. 3e).

We also evaluated whether our malignant cell clusters could be found in independent scRNA-seq datasets. Subtyping of the malignant cell using nearest template prediction (NTP) analysis could re-identify our malignant cell clusters in GSE103322 (77.16%) and GSE164690 (99.10%), respectively (Supplementary Fig. 5); however, our cell subtypes did not cover all the heterogeneity of the malignant cells, which might be due to a small number of patients being sampled. This may support that our malignant cell clusters represent most of the malignant cell types in HNSCC. In addition, we examined the expression of the previous single-cell malignant cell programs of HNSCC in each of the malignant cell clusters21. The CC0 expressed hypoxia and epithelial differentiation-related genes while the CC1 expressed the partial epithelial-mesenchymal transition (p-EMT)-related genes (Fig. 3f). The p-EMT has been suggested to localize in the leading edge in proximity to cancer-associated fibroblasts (CAFs) and are promoted through paracrine interactions between CAFs and malignant cells, resulting in aggressive progression of HNSCC. Supporting this, we could demonstrate that the CC1-high group (>70 percentile of the CC1), compared to the CC1-low group, had more unfavorable clinical outcomes of overall survival (OS; hazard ratio [HR] = 1.38, P = 0.02) and recurrence-free survival (RFS; HR = 1.69, P = 0.006) in TCGA-HNSCC, although the other malignant cell clusters did not (Fig. 3g, top). This finding could be validated using the independent datasets of GSE41613 (OS, HR = 2.30, P = 0.003) and GSE42743 (OS, HR = 2.16, P = 0.016) (Fig. 3g, bottom). Multivariate analyses also revealed the significance of the associations of CC1 with worse prognostic outcomes of HNSCC (GSE41613, OS, HR = 2.68, P = 0.001; GSE42743, OS, HR = 2.33, P = 0.011, Supplementary Table 2). Taken together, we suggest that the abundance of CC1 is associated with the expression of HPV-negative tumor-like features acquiring an aggressive phenotype.

Next, we identified the differentially expressed genes that may represent each cluster (Wilcoxon Rank Sum test P < 0.001, fold difference > 1, Supplementary Fig. 6). We found that the CC1 prominently expressed LGALS7B (Fig. 3h), which could be validated by showing the expression of LGALS7B in the CC1-high samples (i.e., C04, C30, C51, C26, C07, C06, C31, and C15) but not in the CC1-low samples (C38 and C09) (Fig. 3i and Supplementary Fig. 7). Galectin-7 has a tumor-promoting function in diverse cancer types, including HNSCC25,26,27, therefore, it is plausible that the aggressive feature of CC1 might be associated with the expression of p-EMT and LGALS7B.

Fibroblasts-derived COL1A1 expression interacts with CD44 in malignant cells

Fibroblasts are the second most common cell type in the deconvolution analysis (see Supplementary Fig. 4a). We could identify the stepwisely expressed genes during progression from NL to LP and CA tissue in each cell type of epithelial/malignant cells (n = 17, e.g., KRT6A) and fibroblasts (n = 41, e.g., POSTN) (Supplementary Fig. 8), and which could be demonstrated more clearly by a pseudo-time trajectory analysis28 (Fig. 4a). Interestingly, we observed that the pseudo-time of LP was closer to that of CA or LN rather than that of NL, which may indicate that the premalignant LP fibroblasts already acquired a CAF-like feature (Fig. 4b). To further delineate the effects of the transcriptional alteration of fibroblasts on cancer progression, we analyzed the ligand-receptor (LR) interaction between fibroblasts and malignant cells (for details, see Supplementary Methods). We aimed to identify interdependent LR pairs between malignant cells and fibroblasts, which exhibit receptor expression in malignant cells, but whose corresponding ligands are expressed only in fibroblasts and not in malignant cells. We could identify the seven LR pairs showing interdependent and stepwise expression of the ligands in fibroblasts (i.e., COL1A1, COL1A2, COL6A3, THBS1, THBS2, TNC, and LAMA4) and their corresponding putative receptors in malignant cells (i.e., CD44, ITGB1, SDC4, CD47, and ITGA6, Fig. 4c). Remarkably, the COL1A1-CD44 ligand-receptor pair showed the most prominent and interdependent expression (Fig. 4d, top), and which could be validated by Puram’s scRNA-seq data (Fig. 4d, bottom). Moreover, COL1A1 expression was most significantly correlated with the pseudo-time of the fibroblasts (r = 0.57, P < 2.2 × 10−16, Supplementary Fig. 9a). We could successfully demonstrate the correlated expression between COL1A1 in fibroblasts and CD44 in epithelial cells (r = 0.43, P = 7.47 × 10−3, Fig. 4e, and Supplementary Fig. 9b). However, CD44 expression in fibroblasts was not correlated with COL1A1 expression (Supplementary Fig. 9c). These results suggest that the CD44-COL1A1 interaction occurs between fibroblasts and malignant cells but not in an autocrine manner. Immunohistochemical analysis also demonstrated the prominent COL1A1 expression in fibroblasts near the CD44-expressing malignant cells (Fig. 4f and Supplementary Fig. 9d).

a Tissue types of NL, LP, CA, and LN (left) and the pseudo-times (right) of the cells are shown in UMAP plots. b Distributions of the pseudo-time in the fibroblasts across the tissue types are shown. c A heatmap shows the interdependent expression of the ligand-receptor pairs with stepwise expression during HNSCC progression between fibroblasts and malignant cells. d The cell proportions for expressed cells of COL1A1 and CD44 are shown in fibroblasts (F-) and malignant cells (C-) across the tissue types (e.g., F-NL, fibroblasts in NL; C-NL, malignant cells in NL) from our data (top) and that of GSE103322 (bottom). e The correlation plot shows COL1A1 expression levels in fibroblasts and CD44 expression levels in epithelial cells in HNSCC patients. The gray shading represents 95% confidence interval (CI). f An immunohistochemical image shows the expression of collagen and CD44 in fibroblasts (red) and neighbored malignant cells (green), respectively. n = 3. g Invasion of the MSKQLL1 and SCCQLL1 cells treated with siRNAs targeting CD44 (siCD44) or non-target control (NC) in the presence or absence of collagen (1 µg/mL). n = 3. Scale bar, 400 µm. h MSKQLL1 cells and CAF cells are co-cultured with the treatment of NC, siCD44, or siCOL1A1, and the images of the migrated cells (left) and the measured cell population (right) are shown. n = 3 biologically independent experiments. Scale bar, 100 μm. Data are shown as mean ± SD. *,P  <  0.05; **,P  <  0.01; ***,P  <  0.001; ns, nonsignificant by the Student two-tailed unpaired t-test. Source data are provided in a Source Data file.

In addition, we evaluated the functions of COL1A1-CD44 interaction by performing siRNA-mediated knockdown experiments for COL1A1 (siCOL1A1) and CD44 (siCD44). Two HNSCC cell lines expressing CD44 (MSKQLL1 and SCCQLL1, Supplementary Fig. 10a) and the CAFs expressing COL1A1 (CAF48) were selected (Supplementary Fig. 10b). We observed that the treatment with collagen (type I) induced tumor cell invasion, and which was abolished by the transfection with siCD44 (Fig. 4g). In addition, the co-culture of MSKQLL1 and CAF48 enhanced tumor cell migration but did not the co-culture of the siCD44-transfected MSKQLL1 cells with CAF48 or the MSKQLL1 cells with siCOL1A1-transfected CAF48 cells (Fig. 4h). Taken together, we suggest that COL1A1 in fibroblasts and CD44 epithelial/malignant cells interact interdependently, facilitating tumor progression.

CXCL8-expressing CAFs aggravate HNSCC progression

The CAFs from tumor tissues (CA and LN) were further stratified into five clusters (CF0 to CF4), revealing their representative marker genes (Fig. 5a and Supplementary Fig. 11a). When the CAF clusters were compared with the previous CAFs in Puram’s study21, CF0 and CF4 were similar to the Puram’s subtype CAF1 expressing its marker genes (e.g., CTHRC1, COL1A1, COL3A1, POSTN, and MFAP2, Fig. 5b). However, CF4 was different from CF0, showing higher proliferation activity (Supplementary Fig. 11b). CF2 and CF3 were similar to Puram’s CAF2 and myofibroblasts, expressing their marker genes such as CXCL12 and NDUFA4L2. Interestingly, we identified a CAF subtype, CF1, expressing CXCL8 (IL-8 coding gene). CF1 showed a distinct expression pattern differing from the other CAF clusters in the two-dimensional diffusion map embedding analysis29 (Fig. 5c, left). Pseudo-time projection analysis also revealed that CF1 was the last trajectory, implying that CF1 was the most progressed fibroblast (Fig. 5c, right). Moreover, we demonstrated that the abundance of CF1 was significantly correlated with the CXCL8 expression levels in independent datasets of TCGA-HNSCC (r = 0.38, P = 2.72 × 10−18), GSE41613 (r = 0.74, P = 3.41 × 10−18), GSE42743 (r = 0.62, P = 3.54 × 10−9), and GSE65858 (r = 0.68, P = 9.58 × 10−38, Fig. 5d). When we compared the proportion of the CAF subtypes according to the HPV infection status, most CF1 cells were observed in HPV-negative tumors, but not in HPV-positive tumors (P = 1.33 × 10−7, Supplementary Fig. 11c, left). TCGA data also showed that CF1 cells were more abundant in the HPV-negative tumors (10.15%) than HPV-positive tumors (5.04%) (P = 5.26 × 10−7, Supplementary Fig. 11c, right).

Next, we assessed whether the abundances of CF1 and CC1 are associated with each other because both CF1 and CC1 abundances had aggressive features. As expected, CC1 and CF1 abundances were significantly correlated in the pooled HNSCC data (r = 0.41, P = 1.11 × 10−38, Fig. 5e), implying that CF1 and CC1 may interact together. Indeed, galectin-7 (expressed in CC1) promotes metastasis in various cancer types, including HNSCC30. IL-8 (expressed in CF1) is produced by multiple cell types promoting tumor progression31, 32. To evaluate this hypothesis, we demonstrated that galectin-7 treatment significantly induced CXCL8 expression in various CAF cells (CAF30, CAF57, CAF58, and CAF70) as well as HNSCC cells (SCCQLL1 and SNU1066) (Fig. 5f). By contrast, IL-8 treatment did not induce galectin-7 in HNSCC cells. (Supplementary Fig. 12). These results suggest that galectin-7 activates both malignant cells (CC1) and fibroblasts (CF1) to express IL-8, facilitating the aggressive progression of HNSCC.

LAIR2 expression in Treg cells is associated with HNSCC progression

Among the immune cells, NK/T cells were the most common cell type in the scRNA-seq data. The cells with chromosomal aberrations were re-assigned as the non-determined cell type and excluded in the following analyses (see “Methods”). Then, the NK/T cell cluster was further stratified into six subclusters; naïve T cells, CD8+CCL5+, CD4+FOXP3+(regulatory T cell, Treg), CD4+CD154+ (follicular helper T cell), cycling T cells, and NK cells (Fig. 6a, top). The T cell cluster showing higher proliferation activity was designated as cycling T cells, showing higher expression of CD8 coding genes (i.e., CD8A and CD8B) than CD4 as previously described33 (Fig. 6a, bottom, and Supplementary Fig. 13a). During HNSCC progression, the proportions of cycling T cells and CD4+ T cells were increased, whereas the proportions of naïve T cells and CD8+ T cells were decreased regardless of their HPV infection status (Fig. 6b). Pseudo-time trajectory analysis demonstrated the stepwise transition of the T cells according to the tissue types (Fig. 6c). Notably, Tregs in LP exhibited two peaks of pseudo-time distribution, indicating the alteration of the T cell states in the premalignant lesions. In particular, we observed the prominent expression of leukocyte-associated immunoglobulin-like receptor 2 (LAIR2, CD306) in the Tregs in LP (permuted Student’s t-test P = 2.50 × 10−9, fold difference = 2.32, Fig. 6d, and Supplementary Fig. 13b). Although the functions of LAIR2 are largely unknown, LAIR2 has been known as a soluble collagen receptor that activates pro-inflammatory processes34. Therefore, COL1A1 expression in the fibroblasts may activate the LAIR2-mediated downstream signaling in Tregs. In support of this, we observed that COL1A1 expression levels in fibroblasts were weak but significantly correlated with LAIR2 expression levels (r = 0.48, P = 2.48 × 10−3) and the abundance of Tregs (r = 0.38, P = 2.60 × 10−2, Fig. 6e). In addition, we demonstrated that collagen treatment significantly induced FOXP3 expression in mouse Tregs (Fig. 6f and Methods). In addition, by performing fluorescence-activated cell sorting (FACS) analysis, we confirmed that collagen treatment increased the number of both FOXP3+CD25+ Tregs and LAIR2+FOXP3+CD25+ Tregs obtained from the PBMCs in the HNSCC patients (Fig. 6g, h, and Supplementary Fig. 13c, and “Methods”). Taken together, we suggest that the collagen from the fibroblasts interacts with Tregs, inducing LAIR2 and FOXP3 expression and providing a favorable microenvironment for tumor progression.

Discussion

As the next step of TCGA that elucidates molecular heterogeneity using large-scale, multi-omic cancer profiles, a pre-cancer atlas detailing the sequential alterations would be advantageous for delineating driver events during cancer progression. In this study, we aimed to profile the stepwise progression of HNSCC from NL to precancerous LP, CA, and metastatic LN tissues, and demonstrated the single-cell-level transcriptomic dynamics during the stepwise progression of HNSCC. The key findings regarding the interactions among the malignant cells, fibroblasts, and immune cells are summarized in Fig. 7.

The key findings regarding the interactions among the malignant cells, fibroblasts, and immune cells are summarized.

We demonstrated that the LP tissues had a substantial number of CIS cells which showed DNA copy gains and concomitant transcriptional regulation of TP63 and ATP1B3. Indeed, most CIS cells were identified from the patient P15, showing similar CNAs to those of the malignant cells in CA tissues (Supplementary Fig. 14). It is well-known that severe squamous dysplasia or CIS-causing leukoplakia has a higher rate of progression to invasive carcinoma35. Indeed, patient P15 was diagnosed with epithelial dysplasia in 2019, and scRNA-seq was performed at that time. Unfortunately, the patient was re-diagnosed with CIS in 2021. This clinical course of the patient shows that scRNA-seq analysis can help diagnose lesions more accurately and earlier than conventional pathological tests.

In addition, our analyses could identify the cell-cell interactions that promote HNSCC progression. COL1A1-CD44 signaling mediates the interaction between fibroblasts and malignant cells (see Fig. 4). Moreover, a subset of malignant cells (LGALS7B-expressing CC1) interacted with a subset of CAF (CXCL8-expressing CF1) (see Fig. 5). From these results, we suggest that COL1A1-CD44 mediates fibroblast-to-malignant cell interaction, whereas LGALS7B-CXCL8 mediates malignant cell-to-fibroblast interaction.

Immune cells are essential regulatory components in the tumor microenvironment. We demonstrated stepwise alteration of T cell repertoires during HNSCC progression, implying changes in the immune environment before malignancy. Supporting this, premalignant lesions of HNSCC have shown altered inflammatory T cells and increased infiltration of Tregs, promoting tumor progression36, 37. Notably, we observed that Tregs exhibited transcriptional transition from NL to LP, expressing LAIR2. LAIR2 expression has previously been shown to associate with favorable prognostic outcomes in cancer38, whereas the contradictory results have also shown that LAIR2 is a potential marker for T cell exhaustion and worse prognosis in liver cancer39. LAIR2 has a high affinity for various collagen molecules, competing with LAIR1-mediated signaling;34 however, we observed that LAIR1 expression levels were considerably low in NK/T cells; therefore, they are unlikely to compete with each other (Supplementary Fig. 15). Based on our results, we suggest that the COL1A1-LAIR2 interaction between fibroblasts and Tregs provides a permissive microenvironment in tumor development and progression.

In addition to the Tregs, CD154-expressing CD4+ T cells were increased during tumor progression (see Fig. 6b). CD154 is a CD40 ligand that facilitates cell-cell interactions with the CD40-expressing cells, including B cells, macrophages, and follicular helper T cells40. CD154 induces cell proliferation and invasion and poor prognosis, implying that CD154+T cells provide a tumor-promoting immune environment41. By contrast, cycling T cells build up a dysfunctional program and tumor-infiltrating lymphocyte (TIL) activation33. Consistently, we demonstrated that the expression levels of the cycling T cell markers were significantly correlated with the expression levels of dysfunctional T cells markers (r = 0.47, P = 1.5 × 10−28) and the TIL scores (r = 0.3, P = 3.66 × 10−12, Supplementary Fig. 16). Although we did not further evaluate the functions of CD154+CD4 T cells and cycling T cells, our results imply that the alteration of T cell repertoires plays critical roles in the stepwise progression of HNSCC.

We evaluated the cell proportion in tissues using deconvolution analysis. Cell composition in individual samples can be highly heterogeneous according to the sample preparation and experimental procedures. Therefore, the cell proportion obtained from scRNA-seq may not represent the actual cell proportion in tissues. In this study, we compared the cell composition using the deconvolution analyses of the bulk RNA-Seq data, which could estimate the cellular composition in the tissues. Indeed, we identified a large proportion of malignant cells in deconvolution analyses, although the scRNA-seq data detected only a low proportion of malignant cells. Our analysis shows that in addition to measuring the change in the transcriptional expression of cells, measuring the proportional changes of the cell types can reflect the cancer progression according to the microenvironment alteration.

In conclusion, we provide insights into the single cell-level heterogeneity and interactions during the stepwise progression of HNSCC. Targeting the cell-cell interactions that drive the HNSCC progression would be a crucial strategy for the management of HNSCC.