Please refer to prior study for details (11). Briefly, CGHNC8, C9, K2, K6 (32), OC3 (33), OEC-M1 (34), and TW2.6 (35) were kindly provided by researchers at distinct institutions in Taiwan. OC3 and TW2.6 were selected for two independent animal studies. We measured the tumor size and mouse weight twice a week. For Exp1, all tumors were collected on day 68. For Exp2, to obtain tumors with ~500 mm³ in size, tumors of TW2.6-NOG and OC3-NOG were collected on day 32 and day 81, respectively. Total RNAs of each cell line and xenograft tissue were extracted by TRIzol^® reagent (Invitrogen Life Technologies, Carlsbad, CA, USA), cleaned up by RNeasy column (Qiagen, Hilden, Germany), and subjected to an Agilent Bioanalyzer (Agilent, Santa Clara, CA, USA) for RNA Integrity Number (RIN) assessment. Only samples that had an RIN > 7 were selected for mRNA amplification and sequencing (stranded paired-end, Illumina platform, San Diego, CA, USA).

Anti-Pecam1 (ab28364) was purchased from Abcam (Cambridge, UK); anti-Adgre1 (#70076S) was from Cell Signaling (Danvers, MA, USA).

Sections were dewaxed, rehydrated, and incubated with Trilogy™ (Cell Marque, Rocklin, CA, USA) (10 mM citrate buffer, pH 6.0 for Adgre1) at 121°C for 10 min to unmask antigens. At room temperature (RT), the slides were immersed in 3% hydrogen peroxide for 15 min to quench endogenous peroxidase activity followed by 1% bovine serum albumin for 60 min to block nonspecific antigenic sites. Slides were incubated with indicated primary antibodies at 4°C overnight (Pecam1, 1:100, Adgre1, 1:1500). After washing with 1X TBS containing 0.05% Tween 20, slides were incubated with horseradish peroxidase-conjugated secondary antibodies and developed by chromogen diaminobenzidine using the DakoReal™ EnVision™ kit (#K5007, DAKO, Glostrup, Denmark). All slides were counterstained with Mayer’s hematoxylin and scanned by a Pannoramic MIDI scanner (3DHISTECH, Budapest, Hungary).

Alcian blue staining kit (ab150662, Abcam, Cambridge, UK) was used to detect hyaluronan deposition in the xenograft tissues according to the manufacturer’s instruction. Briefly, dewaxed and rehydrated FFPE tissue sections were incubated with Alcian Blue Solution for 30 min at RT, slides were counterstained with Nuclear Fast Red Solution for 5 min at RT and scanned by a Pannoramic MIDI scanner (3DHISTECH).

To quantitate the immunostaining of Pecam1, Adgre1, and Alcian blue staining in the xenograft tissues, scanned images at 100× magnification were digitized and quantitated by using ImageJ plugged-in with the Immunohistochemistry Image Analysis Toolbox (v1.40p) (NIH, Bethesda, MD, USA). Quantitation results were visualized by R package ggplot2 (v3.3.3). Statistical differences between OC3 and TW2.6 groups were evaluated by two-sample t-test of means (compare_means), as denoted in each plot.

Cell line dataset (GSE150469): human reference genome (hg19) aligned reads for individual genes (n=21,916) were used to compute expression values in transcripts per kilobase million (TPM) by using Cufflinks (v2.1.1). R pheatmap (v1.0.12) was used to visualize the relative expression levels of indicated genes in each cell line. Xenograft tissue dataset (GSE149496): as described previously (11), R XenofilteR (v1.8) processed human (hg38) and mouse (mm10) aligned reads for individual genes were used to perform gene quantification in TPM values. The resulting gene numbers for human and mouse are 17,759 and 16,374, respectively. In the exploratory analysis, R package limma (v3.40.6) was used to compute differentially expressed genes followed by Volcano plots for visualization. To identify biological pathways enriched in the OC3- and TW2.6 tumor stroma, the mouse expression matrix was subjected to gene set enrichment analysis (GSEA) (36). Enrichment was considered significant when false discovery rate (FDR) was less than 5%. R ggcorrplot (v0.1.3) was used to calculate and visualize the correlation matrix comprising genes of interest. The Pearson correlation coefficient was computed; p-value < 0.05 was considered significant. R ggplot2 (v3.3.3) was used to visualize the relative expression levels of indicated genes in the OC3- and TW2.6 stroma. Unpaired two-sample t-test of means was applied to evaluate statistical differences; p-value < 0.05 was considered significant. R pheatmap (v1.0.12) was used to visualize the relative expression levels of indicated genes in each tissue. In silico enumeration of cell fractions was conducted by using CIBERSORTx according to its online documentation of which the ‘single cell RNA-seq HNSCC’ and ‘LM 22’ were used as signature matrices, respectively (27). Puram scRNA-seq dataset (GSE103322): the raw expression matrix was acquired from the UCSC Cell Browser portal. Expression values (TP100K) were normalized, scaled, and log-transformed by using R Seurat (v4.0.3) (37). DotPlot was used to visualize the average expression and fraction of indicated genes across eleven cell types comprising 5,902 cells. NCKU-OrCA-40TN dataset (GSE37991): this normalized microarray dataset comprises 18,047 genes for further analysis. R ggcorrplot was used to calculate and visualize the correlation matrix composed of 20 selected stromal genes. The Pearson correlation coefficient was computed; p-value < 0.05 was considered significant. R survival (v3.2-7) was used to assess the univariate and multivariate Cox proportional hazards of indicated clinical features and genes. The Cancer Genome Atlas (TCGA) datasets: for each indicated cancer type, the expression matrix of 20 stromal genes and associated clinical information were acquired from the UCSC Xena platform (38), followed by univariate Cox proportional hazards assessment and dichotomized Kaplan–Meier overall survival curves prediction by using R survival (v3.2-7).

A recent consensus statement of epithelial-mesenchymal transition (EMT) research reiterated the importance of associating cellular characteristics, rather than the expression of a single or a small set of molecular markers, with EMT phenotypes (6). To comply with these guidelines, we inspected the expression levels of core EMT transcription factors, p-EMT hallmarks, and canonical epithelial markers and mesenchymal regulators (6, 9, 10, 39, 40) in a set of seven Taiwanese cell lines derived from the oral cavity. Among these, while the oral cancer cell line OC3 has the highest mesenchymal propensity, TW2.6 and OEC-M1 are hypothetically p-EMT cells since they maintain both epithelial and mesenchymal genes ( Figures 1A, B ). Indeed, we demonstrated that TW2.6, but not OC3 or OEC-M1, displayed p-EMT multicellular characteristics in vitro and in vivo (10). In addition, mesenchymal OC3 repeatedly grew slower and smaller than its p-EMT TW2.6 counterpart in vivo [ Figure 1C and ref (11)], a phenomenon consistent with one prior study in that head and neck cancer tissues of the inflammatory mesenchymal subtype had a better prognosis (24).

Next, exploratory analysis of differentially expressed genes (DEGs) in the OC3 and TW2.6 xenograft tissues revealed that (1) compared to the OC3 cells expressing various innate immunity responsive genes, TW2.6 cells expressed MYC and E2F targets involved in cell proliferation and cell adhesion ( Supplementary Figure 1A ); (2) compared to the OC3 stroma harboring various extracellular matrix (ECM)- and TGFβ axis-related transcripts, the TW2.6 stroma was characterized by proangiogenic factors and immune-related genes. Gene set enrichment analysis (GSEA) also confirmed that while the epithelial-mesenchymal transition hallmark is uniquely enriched in the OC3 stroma, angiogenesis and immune-related molecular processes are recurrently detected in the TW2.6 stromal compartment ( Supplementary Figure 1B ).

To validate that the p-EMT TW2.6 stroma had higher angiogenesis processes than the OC3 group, we first performed correlational analysis of stromal Pecam1 expression, a microvessel density surrogate, with that of each well-established proangiogenic and endothelial index gene recently established in > 10,000 human tumors (41). The results indicated that the expression of tumor cell-derived VEGFA, IGFBP3, EFNA1, EFNB1, IGF2, PDGFA and stroma-derived Angpt2, Cdh5, Esam, Esm1, Icam2, and Tie1 was statistically correlated with that of Pecam1 ( Figure 2A ) and elevated in p-EMT TW2.6 tumors ( Figure 2B ). In parallel, immunohistochemistry using a Pecam1 antibody not only verified significantly increased vascular densities but also prominent angiolymphatic invasions in the p-EMT TW2.6 tissues compared to their mesenchymal OC3 counterparts ( Figure 2C ). Taken together, these results strongly suggest that a higher blood supply from the tumor stroma might contribute to fostering better growth of p-EMT TW2.6 cells in vivo.

The other prominent expression feature enriched in the TW2.6 stroma is immune-related molecular signatures, including various cytokines and chemokines (e.g., Il6, Cxcl9, Cxcl10, Cxcl12). At first, this seemed perplexing since prior results showed that stronger innate immunity was present in OC3 tumor cells (11). To resolve this confusion, we performed in silico enumeration of cell fractions available at the CIBERSORTx portal (27). We first employed ‘scRNA-seq HNSCC’ as a signature matrix, since both OC3 and TW2.6 cells are derived from the oral cavity. As expected, while no significant difference was detected in the tumor compartment (primarily consisting of ‘malignant’ cells), the ‘fibroblast’ cell fraction was significantly higher in the OC3 stroma than in the TW2.6 stroma, which is consistent with our prior results (11) ( Supplementary Figures 2A, B ). Importantly, the constituents of other immune cell populations did not differ significantly between the two groups. Alternatively, using LM22, the default signature matrix comprising cell type-specific genes derived from 22 leukocytes revealed two immune populations, i.e., ‘Macrophage-M1’ and ‘Tγδ‘, were statistically elevated in the TW2.6 stroma. Further delineation of hallmark genes comprising each LM22 cell type revealed that instead of authentic immune cell markers, increased expression of stromal Cxcl9, Cxcl10, Cxcl11 (Macrophage-M1) and Ccl5 (Tγδ) were likely to be the leading genes contributing to the statistical enrichments ( Supplementary Figures 2C, D ).

Given that ‘macrophage’ represents the major (57.2 ± 6.7%) immune cell type residing in the host stroma by CIBERSORTx analysis, immunohistochemistry using an antibody against the panmacrophage marker Adgre1 (F4/80) was performed to visualize the infiltration of macrophages in each tumor section. Unexpectedly, the spatial locations of macrophages were significantly different between the OC3 and TW2.6 groups. In the majority of OC3 tissues, the tumor cells are interdigitated with irregular Adgre1-positive macrophages and are accompanied by high staining backgrounds. By contrast, in the TW2.6 group, Adgre1-positive prototypical macrophages were frequently detected in the tumor margins, a phenomenon reminiscent of ‘immune privileged’ sites detected in clinical head and neck cancer tissues (42). Note that no significant difference in macrophage staining was noticed between the OC3 and TW2.6 groups after background subtractions ( Figures 3A, B ), which is consistent with the CIBERSORTx results. Intriguingly, as depicted in Figure 3C , the simultaneous upregulation of murine signature genes of tumor-associated macrophages (43) and downregulation of major histocompatibility class II molecules required for antigen presentation in TW2.6 tissues reinforces the accumulation of anergic tumor-associated macrophages surrounding p-EMT TW2.6 tumor cells.

Recent scRNA-seq studies have shown that the inflammatory subtype of CAFs (iCAFs) secretes a variety of inflammatory mediators (e.g., IL6 and CXCL12), whereas the myofibroblastic counterpart (myCAFs) primarily expresses ECM molecules (e.g., FN1 and type I collagen) (18–23). Of special interest, by microarray-based transcriptomics, Costea et al. identified two distinct subtypes of CAFs from clinical oral cancer tissues, designated CAF-N (normal-like) and CAF-D (divergent) (44). CAF-N is intrinsically motile, secretory, proangiogenic, and hyaluronan-rich. CAF-D is less migratory and secretes high levels of TGFβ but is unresponsive to it. Compared to CAF-D, the CAF-N population displayed a greater (85.33% vs. 50%) and faster (7 vs. 14 days) tumor-promoting incidence of otherwise nontumorigenic oral dysplastic cells in the immunodeficient mice (NSG) and significantly deeper invasion of malignant cells in the 3D biometrices constructed in vitro. The authors proposed that a switch from an earlier secretory CAF-N to a later TGFβ^high CAF-D occurred during oral cancer progression.

To test the hypothesis that inflammatory cytokines detected in the TW2.6 stroma might come from the iCAF cell population, in GSEA we applied customized gene matrix transposed (gmt) file of each newly identified stromal cell population as summarized in Table 1 . Interestingly, this approach explicitly assigned enrichment of myCAF and CAF-D to the OC3 stroma and iCAF, tumor-associated macrophages, endothelial cells, and perivascular-like cells to the TW2.6 stroma ( Figure 4A ). In addition, a great number of inflammatory mediators are overlapping genes of the iCAF and immune-related pathways (e.g., IL6/JAK/STAT3 and inflammatory responses), which supports our hypothesis that CAFs of the TW2.6 stroma might be the principle source for the immune-related signatures ( Figure 4B ).

We were aware that two caveats are inherent to the current study: (1) immunodeficient NOG xenografts were used, which retained only limited innate immunity, and (2) human orthologous genes might not be entirely functional in mice. With these concerns, we performed cell type-gene expression analysis using the scRNA-seq dataset of oral cancer tissues in which the presence of CAF subtypes, including myCAFs and iCAFs, has been previously noted (8, 19, 21, 23). The results showed that except for IL6, which was enriched in the iCAF population of xenograft tissues but myofibroblasts in the clinical samples, hallmark genes for myCAFs, iCAFs, and myofibroblasts were fairly consistent between the xenograft and clinical oral cancer tissues, suggesting a certain extent of conservation in non-immune stromal cell populations preserved in the current study ( Figure 4C ).

Another notable feature repeatedly linked to the secretory iCAF population is increased levels of hyaluronan synthases, including HAS1 and HAS2 (18, 19, 44). Interestingly, in our system, significantly increased expression of synthases (HAS2, HAS3) and hyaluronidases (HYAL1, HYAL3, HYAL4) was detected in the tumor cells but not in the stromal portion of TW2.6 tumors ( Supplementary Figure 3A ). Subsequent experiments using Alcian blue staining of OC3 and TW2.6 xenograft tissues confirmed that prominent hyaluronan staining was detected in five out of six TW2.6 tissues compared to sporadic Alcian blue-positive mast cells revealed in the OC3 tissues ( Supplementary Figure 3B ). Thus, our data partly support prior scRNA-seq studies in that iCAFs are frequently present in a hyaluronan-rich tumor microenvironment.

Thus far, our results provide evidence that certain stromal cell populations are preferentially associated with the most invasive yet uncommon EMT (OC3) or p-EMT (TW2.6) tumor cells, which is very different from a real clinical specimen that comprises cells at variable EMT states, e.g., only tumor cells at the invasion front display a p-EMT phenotype (6). A complementary experiment is to correlate stromal genes of interest with transcriptomic datasets from clinical bulk tissues, an approach we employed to identify the CAF index (summed expression of TGFBI, TGFB2, TGFBR2, FN1) (11). Specifically, the Puram oral cancer scRNA-seq (GSE103322) and NCKU-OrCA-40TN (GSE37991) datasets were chosen for cell type mapping and survival analysis, respectively. Of note, the latter is a microarray dataset comprising 40 matched pairs of betel quid-associated oral squamous cell carcinoma and adjacent normal tissues. Due to treatment-related death, only 38 cases were included for survival analysis. In addition, except for recurrence, none of the other clinical features were statistical covariates for overall survival in univariate Cox proportional hazards model assessment ( Table 2 ), an inherent limitation of a small cohort. Further TCGA collections with larger sample sizes will be included as validation datasets to complement this caveat.

First, genes identified in the xenograft experiments were validated for their primary origins of expression by cell typing of the scRNA-seq database, including myofibroblasts (IL6, ACTA2, MYL9, TAGLN), myCAFs (TGFBI, TGFB2, FN1, THY1, COL12A1, FAP, HAS1, HAS2), iCAFs (FAP, CXCL12, GSN), endothelial cells (PECAM1, TGFBR2, HYAL1, HYAL2), mast cells (HYAL3), T cells (HYAL4) and tumor cells (HAS3) ( Figures 4C , 5A ). Next, in the 40 tumor tissues of NCKU-OrCA-40TN, the Pearson correlation coefficient of each paired gene was computed and clustered by cell type to inspect their expression consistency. As depicted in Figure 5B , (1) the expression of myCAF genes correlated with each other the best. (2) Myofibroblast genes share partial similarities to iCAFs and myCAFs, suggesting myofibroblasts might be progenitors of both CAF subtypes. (3) The expression of PECAM1, a microvessel density proxy, was inversely correlated with myCAF-HAS2 and TGFBI but significantly correlated with iCAF-CXCL12 and myofibroblast-ACTA2. (4) The expression of HYAL1 is inversely related to myCAF genes but positively related to tumor cell HAS3.

Third, a univariate Cox proportional hazards model was performed to compute the hazard ratios (HRs) of overall survival associated with each stroma-originating gene in the NCKU-OrCA-40TN cohort. Accordingly, 6 out of 20 selected genes revealed statistical significance, including TGFBI (HR 11, 95% CI 2.4–50, p=0.0021), TGFB2, COL12A1, IL6, HAS2, and HYAL1 (HR 0.15, 95% CI 0.042–0.55, p=0.0039) ( Table 2 ). Multivariable Cox analysis revealed these 6 genes are confounding covariates ( Figure 5C upper panel), which was not unexpected since 4 of them (TGFBI, TGFB2, COL12A1 and HAS2) are derived from myCAFs. As cancer recurrence is a known prognostic factor for overall survival, we included it in multivariate Cox analysis. The results showed that TGFBI and HYAL1, respectively, remained statistically significant in multivariate Cox models ( Figure 5C lower panel). These data suggested that both genes could serve as robust prognostic factors in betel quid-oral cancer. It is worth noting that TGFBI and HYAL1 should act independently rather than combined, as both lost partial significance in bivariate Cox assessment (TGFBI 8.23 95% CI 1.66–40.79, p=0.01, HYAL 0.23, 95% CI 0.06–0.88, p=0.032).

To validate our findings in larger TCGA cancer datasets, we first performed univariable Cox analysis and Kaplan–Meier curves stratified by low- and high-expression groups for each of the 20 stromal genes in the head and neck cancer cohort (HNSC, n=519). The results indicated that while high expression levels of TGFBI, FAP and IL6 were statistically prognostic for poorer survival; the expression of HYAL1 did not reach statistical significance for better survival ( Figure 6A ). It is worth noting that HNSC comprises 11 anatomic subsites including tonsil and larynx, such heterogeneity in tissue source might interfere with precisely identifying stromal biomarkers specific to each HNSC subtype. Next, same approach was used to evaluate TGFBI and HYALI in the other 32 TCGA cancer types. Interestingly, kidney renal clear cell carcinoma (KIRC, n=533) and uveal melanoma (UVM, n=80) concurrently displayed TGFBI and HYAL1 as poor and good biomarkers, respectively ( Figure 6B ). In addition, the expression of TGFBI was also a poor indicator for bladder urothelial carcinoma (BLCA, n=406), breast invasive carcinoma (BRCA, n=841), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC, n=302), and glioblastoma multiforme (GBM, n=153); and the expression of HYAL1 was a statistical protective marker for kidney renal papillary cell carcinoma (KIRP, n=287) and pheochromocytoma (n=148) ( Figure 6C ). Taken together, the prognostic values of TGFBI and HYAL1 identified from betel quid-associated oral cancer tissues ( Figure 6D ) were recapitulated by seven (n=2834) and four (n=1048) TCGA cancer types, respectively.