Mesenchymal-epithelial transition in lymph node metastases of oral squamous cell carcinoma is accompanied by ZEB1 expression

Abstract

Background: Oral squamous cell carcinoma (OSCC), an HPV-negative head and neck cancer, frequently metastasizes to the regional lymph nodes but only occasionally beyond. Initial phases of metastasis are associated with an epithelial-mesenchymal transition (EMT), while the consolidation phase is associated with mesenchymal-epithelial transition (MET). This dynamic is referred to as epithelial-mesenchymal plasticity (EMP). While it is known that EMP is essential for cancer cell invasion and metastatic spread, less is known about the heterogeneity of EMP states and even less about the heterogeneity between primary and metastatic lesions.

Methods: To assess both the heterogeneity of EMP states in OSCC cells and their effects on stromal cells, we performed single-cell RNA sequencing (scRNAseq) of 5 primary tumors, 9 matching metastatic and 5 tumor-free lymph nodes and re-analyzed publicly available scRNAseq data of 9 additional primary tumors. For examining the cell type composition, we performed bulk transcriptome sequencing. Protein expression of selected genes were confirmed by immunohistochemistry.

Results: From the 23 OSCC lesions, the single cell transcriptomes of a total of 7263 carcinoma cells were available for in-depth analyses. We initially focused on one lesion to avoid confounding inter-patient heterogeneity and identified OSCC cells expressing genes characteristic of different epithelial and partial EMT stages. RNA velocity and the increase in inferred copy number variations indicated a progressive trajectory towards epithelial differentiation in this metastatic lesion, i.e., cells likely underwent MET. Extension to all samples revealed a less stringent but essentially similar pattern. Interestingly, MET cells show increased activity of the EMT-activator ZEB1. Immunohistochemistry confirmed that ZEB1 was co-expressed with the epithelial marker cornifin B in individual tumor cells. The lack of E-cadherin mRNA expression suggests this is a partial MET. Within the tumor microenvironment we found immunomodulating fibroblasts that were maintained in primary and metastatic OSCC.

Conclusions: This study reveals that EMP enables different partial EMT and epithelial phenotypes of OSCC cells, which are endowed with capabilities essential for the different stages of the metastatic process, including maintenance of cellular integrity. During MET, ZEB1 appears to be functionally active, indicating a more complex role of ZEB1 than mere induction of EMT.

Keywords: EMT; Epithelial–mesenchymal plasticity; Heterogeneity; MET; Oral cavity; Partial EMT; Single cell RNA; Squamous cell carcinoma; ZEB1.

Fig. 1

Single-cell gene expression signatures in OSCC cells from a single metastasis reveal predominant functional phenotypes. A UMAP based on scRNAseq data of 4076 cells isolated from a metachronous lymph node metastasis. Cells are annotated and summarized according to the presumed cell type. B UMAP of 1906 OSCC cells depicted in A. Cells are annotated according to predominant functional phenotype. C Heatmap for scaled, log-normalized gene expression of tumor cells (columns) split by respective phenotype and the top 10 differentially expressed genes (DEGs) (rows) of the respective phenotype against all other tumor cells. DEGs are sorted from highest to lowest log2 foldchange. Row sections are ordered like column sections. D Top 5 enriched gene sets from log2 foldchanges of respective tumor phenotypes by normalized enrichment scores (x-axis). Gene sets of respective phenotypes are sorted from highest to lowest enrichment. Bars are colored by the negative decadic logarithm of the Benjamini–Hochberg adjusted p-value (padj). DCs: dendritic cells. ECs: endothelial cells

Fig. 2

A progressive epithelial differentiation, but no strong uniform direction of development in pEMT clusters. A UMAP of 1906 OSCC cells annotated based on SNN clustering, defining 4 pEMT (pEMT-1 to 4), 4 epithelial differentiated (epi-1 to 4) and one mixed (mix) cluster; clusters are numbered by size. B Heatmap for scaled, log-normalized gene expression in EMP-associated tumor cell phenotypes (columns) split by EMP cluster and their top 5 DEGs (rows) against all other EMP-related tumor cell phenotypes. DEGs are sorted from highest to lowest log2 foldchange. Row sections are ordered like column sections. C Projection of RNA velocity on the UMAP depicted in A. Arrows indicate the extrapolated direction of development; arrow length indicates strength of future development. D First two principal components of OSCC cells with the three EMP-related principal curves that are derived from trajectory inference. Graph on top visualizes the relationship between EMP clusters described by the three principal curves forming a branching trajectory. E Log-normalized expression (y-axis) of MMP1, VIM, SPRR1B and KLK7 across pseudotime values (x-axis) of curve 2, color-coded by clusters. Red lines indicate smoothed expression values over the trajectory generated with a general additive model; 95% confidence intervals are shaded gray. F Inferred CNVs across EMP-related tumor cells (rows) for all chromosomes (columns). Red indicates copy number gains, white diploid copy number and blue copy number loss. Columns show genes categorized in chromosomes and ordered by genome position; hence the size of the chromosome reflects the number of detected genes and not its nucleotide length. Mitochondrial genes were excluded

Fig. 3

Intra-tumoral heterogeneity of OSCC is driven by EMP. A UMAP based on scRNAseq data of 7263 cancer cells from 16 different patients annotated by patient. B Heatmap for scaled, log-normalized gene expression of tumor cells split by patients and their top 5 DEGs (rows) against all other tumor cells. All patients with less than 100 cells are summarized in the ‘other’ column. DEGs are sorted from highest to lowest log2 foldchange. Row sections are ordered like column sections. C UMAP based on scRNAseq data depicted in A with PCs corrected for patient-specific effects using harmony. Cells are annotated according to their predominant phenotype. D Relative distribution of tumor cell phenotypes (left) and cancer cell abundance (right) across patients. The label on the y-axis shows the sample identification and tumor localization (primary tumor [PT] or lymph node metastasis [MET]). E UMAP based on scRNAseq data of 2948 OSCC cells from patient HN01. Cells were annotated based on SNN clustering and the predominant phenotype. F Triangle heatmap of cosine similarity comparing the intratumoral heterogeneity across all patients. Cosine similarity is calculated between log2 fold changes from patient-specific clusters against all other tumor cells within the respective patient. Left side annotated are patient-specific clusters from patient #1 depicted in Fig. 2A and right side from patient HN01 depicted in E. We included only patients with more than 50 tumor cells

Fig. 4

ZEB1 is highly active in metastatic epithelial differentiated OSCC cells. A Percentage of tumor cells with detectable mRNA expression from scRNAseq (more than one UMI) encoding the indicted EMP-related transcription factors. B Mean inferred activity based on the target genes of the indicated transcription factors across tumor phenotypes from EMP-related patient-specific clusters. On top the log-normalized expression of CDH1 and cornifin-B (SPRR1B) is shown, on bottom the localization (primary tumor [PT] or lymph node metastasis [MET]) and respective EMP phenotype of the cluster. C Mean activity of ZEB1 for epithelial differentiated and pEMT clusters of each patient, respectively. Connecting lines show dots belonging to the same patient. ZEB1 (D) and cornifin-B (E) protein expression detected in serial sections by IHC of the primary tumor from patient #2; comparable areas are depicted. Scale bars equal 200 µm in overview and 100 µm in zoomed image. E Colocalized expression of ZEB1 (green) and cornifin-B (red) detected by double staining in the lymph node metastasis of patient #1. Nuclei are stained in blue (DAPI), Scale bars equal 10 µm

Fig. 5

The OSCC microenvironment is composed of heterogenous fibroblasts from which immunomodulatory cells are present across the metastatic cascade. A Overview of types of analyzed OSCC samples and their localization within the head and neck area. B Number of samples across patients colored by their respective tissue origin; for patient #4, the primary tumor could not be analyzed due to incorrect specimen processing; for patients #6 and #7 two regions of the primary tumor were analyzed, denoted as sample #6.1 and #6.4 and #7.1 and #7.4, respectively. C UMAP of 41,284 cells based on OSCC scRNAseq data from our cohort and colored by cell type. D UMAP of 21,037 cells based on CD45-negative and HPV-negative primary HNSCC from Kürten et al. and colored by cell type. E UMAP of 1,595 fibroblasts and 551 pericytes from C colored by the respective phenotypes derived from shared-nearest neighbor clusters. F UMAP of 2,920 fibroblasts and 683 pericytes from D colored by the respective phenotypes derived from shared-nearest neighbor clusters