INHBA⁺ macrophages and Pro-inflammatory CAFs are associated with distinctive immunosuppressive tumor microenvironment in submucous Fibrosis-Derived oral squamous cell carcinoma

Abstract

Transcriptomic and metabolic profiles of tumor cells and stromal cells in oral squamous cell carcinoma (OSCC)-derived from oral submucosal fibrosis (OSF) (ODSCC) have been reported. However, the complex intercellular regulatory network within the tumor immunosuppressive microenvironment (TISME) in ODSCC remains poorly elucidated. Here, we utilized single-cell RNA sequencing (scRNA-seq) and spatial transcriptomics (ST) data from GEO database and multiple immunofluorescence staining (mIF) to reveal distinctive TISME of ODSCC. Results found that compared to OSCC without OSF history (NODSCC), OSCC derived from OSF (ODSCC) showed a significant increase in exhausted CD8⁺T and Treg cells (Ro/e > 1, p < 0.05) and a decrease in cytotoxic T (CTL) (Ro/e < 1). ODSCC enriched in more Inhibin subunit beta A⁺ Macrophages (INHBA⁺Mac) and Proinflammatory Cancer-associated Fibroblast (iCAF) versus NODSCC. INHBA⁺Mac possessed strongest immune-suppressive functions, evidenced by highest immune checkpoint scores, lowest MHC scores and highest expression of SPP1 among macrophages. Moreover, INHBA⁺Mac in ODSCC presented stronger immune-suppressive functions than that in NODSCC. iCAF differentially highly expressed INHBA and enriched in immune-related pathways and collagen/ECM pathways across CAF subsets, and possessed stronger immune-suppressive functions, as shown by up-regulated gene expression of TDO2, IDO1 and DUSP4 in ODSCC versus in NODSCC. Furthermore, INHBA expression was higher in ODSCC than in NODSCC (p < 0.01). The classic OSF-inducing molecule arecoline significantly increases the expression of INHBA (p < 0.0001) in vitro experiments stimulating THP-1 cells. ST analysis revealed a close co-location of INHBA⁺Mac, iCAF and Treg and SpaGene identified INHBA-ACVR1/ACVR2A/ACVR2B interaction regions overlapping with distribution of three types of cells. Collectively, ODSCC shows a more severe TISME and potentially poorer sensitivity to immunotherapy than NODSCC. The increased INHBA⁺Mac and iCAF in ODSCC are associated with the observed more severe TISME. The upregulated INHBA in ODSCC and its interaction with INHBA-ACVR1/ACVR2A/ACVR2B may mediate the modulation effect of INHBA⁺ Mac and iCAF on Treg differentiation and functionality.

Introduction

Oral squamous cell carcinoma (OSCC) is the most common malignant tumor in the oral and maxillofacial region. According to Global Cancer Statistics 2022, there are 389,458 new cases annually worldwide, with 188,230 deaths attributed to this disease [1]. The development of OSCC is affected by a variety of complex factors including heavy use of tobacco, betel quid chewing, consumption of alcoholic beverages, and chronic inflammation [2]. Some oral mucosal lesions such as leukoplakia, erythroplakia, oral submucosal fibrosis (OSF) and lichen planus are regarded as oral potential malignant disorder (OPMD) [3]. It has been demonstrated that in Southeast Asia and in South of China, there is a high prevalence of OSCC associated with OSF [4].

Previous studies have revealed distinct clinicopathological profiles between OSF-derived OSCC (ODSCC) and non-OSF-associated OSCC (NODSCC) [5, 6]. However, comparative prognostic analyses yield contradictory findings. Pankaj et al. [6] reported superior disease-specific survival (DSS) in ODSCC compared to NODSCC, attributing this to its favorable clinicopathological features and improved oncological outcomes. Similarly, Divya et al. [7] observed earlier tumor staging, better differentiation, and enhanced prognosis in ODSCC. In contrast, Feng et al. [5] demonstrated that ODSCC exhibits heightened clinical aggressiveness, increased metastatic potential, and poorer survival rates. These conflicting data underscore the necessity of further studying distinct pathogenesis of ODSCC.

The tumor microenvironment (TME) plays a crucial role in tumor malignancy, immune evasion, and therapy resistance [8]. Tumor-infiltrating T cells are often in a state of exhaustion, reflecting a tumor immunosuppressive microenvironment (TISME). Macrophages, cancer-associated fibroblasts (CAFs), and endothelial cells within the TME secrete various cytokines that shape the immune landscape while promoting tumor proliferation, invasion, and metastasis [9]. Nevertheless, these cells present significant heterogeneity. The integration of single-cell RNA sequencing (scRNA-seq) and spatial transcriptomics (ST) are commonly used to unveil this kind of heterogeneity and interactions of different cell types. Studies by Kurkalang et al. [10] and Zhi et al. [11] through scRNA-seq and ST analysis have reported the transcriptomic and metabolic profiles of tumor cells, CAFs, and immune cells and highlighted the critical roles of the p-EMT process and metabolic reprogramming in ODSCC. Yet, immunological characteristics of TME, especially macrophage and CAF subtypes and their key responsible molecules associated with immunosuppression in ODSCC, remains poorly elucidated.

INHBA, a key subunit of activin A and a member of the TGFβ superfamily has been demonstrated to be overexpressed in multiple solid tumors (e.g., colorectal, gastric, and ovarian cancers) and significantly correlated with tumor invasion, metastasis, and poor prognosis [12,13,14]. Additionally, it induces cancer-associated fibroblasts (CAFs) to secrete IL-6 and VEGF, fostering angiogenesis and immunosuppression [15]. However, whether and how INHBA plays immunosuppressive role in ODSCC keep to be explored.

In this study, we utilized scRNA-seq and ST data from the GEO database and experimental validation to reveal the distinctive TME landscape of ODSCC. By comparing ODSCC with NODSCC at the single-cell level, we identified that the high proportion of inhibin subunit beta A⁺ macrophages (INHBA⁺Mac) and proinflammatory cancer-associated fibroblast (iCAF) that highly expressed INHBA was significantly associated with the formation of a more potent immunosuppressive microenvironment, influencing tumor progression in ODSCC. Furthermore, our findings suggest that INHBA-driven SMAD signaling activation contributes to TISME formation, positioning INHBA as a potential therapeutic target for OSCC, particularly ODSCC.

Materials and methods

Data collection

We downloaded the GSE215403, GSE208253, and GSE220978 datasets from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), containing scRNA-seq data of 12 samples in the GSE215403 dataset, 9 samples from the NODSCC group, and 3 samples from the ODSCC group. In addition, ST data of 4 NODSCC samples from the GSE208253 dataset, and 4 ODSCC samples from the GSE220978 dataset were used for analysis, after performing data consistency processing between two groups.

scRNA-seq data preprocessing and integration

For gene expression sequencing, the downloaded count matrices were imported into the R package Seurat (v4.1.0). Samples in GSE215403 were merged into a single Seurat object for consistent filtering. After quality control, including removing cells with gene counts less than 200 and exceeding 5000 or cells with abnormally low or high UMI counts and high mitochondrial read percentages, genes from red blood cells and any remaining multiplets expressing mutually exclusive marker genes, “NormalizeData”, “FindVariableFeatures”, and “ScaleData” were applied to normalize the scRNA-Seq data.

Cluster annotation and data visualization

Normalized and filtered data were processed using the standard Seurat pipeline (v4.1.0). TSNE dimensionality reduction was used for visualization, and Seurat’s “FindClusters” function (v4.1.0) was used to separate cells into unsupervised clusters. Cell types in clusters were defined using the marker genes from references.

Differentially expressed genes (DEGs) analysis

Genes specific to each cluster or group were identified using the “FindAllMarkers” function, and adjusted p-values were calculated using the Wilcoxon rank-sum test. Volcano plots and heat maps were used to show the fold changes and log-adjusted p-values for DEGs.

Analysis of Single-Cell trajectories

Developmental pseudotemporal ordering of single cells was inferred through the Monocle2 computational framework within the R statistical environment (v4.1.0).The “newCellDataSet”, “estimateSizeFactors”, and “estimateDispersions” were used to perform these analyses. The “detectGenes” was used to filter low quality cells with “min_expr = 0.1”.

Cell–cell communication analysis

The R package CellChat was utilized to analyze cell–cell communication, calculating the total number of ligand–receptor interactions and cell-cell interactions among cell types [16].

Functional enrichment analysis

The GO and KEGG pathways were analyzed using the ClusterProfiler R package. Analysis was performed by GSVA and GSEA algorithm using “c2.cp.kegg.v7.4.symbols.gmt”, “c5.go.bp.v7.4.symbols.gmt” and “h.all.v7.4.symbols.gmt” in MSigDB to get the differences in enrichment pathways between different groups [17].

Tissue distribution of specific cell subtypes

Tissue preference of each cluster was estimated by the STARTRAC-dist index, in which Ro/e denotes the ratio of observed to expected cell number. Re/o indicates whether cells of a certain subcluster are enriched or depleted in a specific tissue [18].

CIBERSORTx algorithm

CIBERSORTx algorithm is employed to digitally “purify” cell-type-specific expression profiles from datasets obtained from GEO and TCGA by utilizing our scRNA-seq-derived reference profiles [19].

10 Survival analysis in TCGA HNSC data set

Prognostic outcome assessment was conducted via the GEPIA2 web-based platform (http://gepia2.cancer-pku.cn/#index), which demonstrated significant differential survival outcomes through log-rank testing and confirmed expression differences using Student’s t-test.

11 Spearman’s correlation analysis

Nonparametric Spearman’s rank correlation was employed to quantify associations between immune cell infiltration levels, with statistically significant relationships defined by an absolute coefficient threshold (|Rs| > 0.3) and Benjamini-Hochberg adjusted p-values < 0.05.The R package ggpubr was used to assess and visualize the correlation of INHBA⁺Mac and Treg in GSE65858.

12 Single-cell copy-number variation (CNV) evaluation

The CNV evaluation of each cell was conducted by infercnv R package. The CNVs of Epithelial cells were calculated and the immune cells were applied as the reference. The inferCNV analysis was performed with parameters including “denoise”, default hidden Markov model (HMM) settings, and a value of 0.1 for “cutoff”.

SCENIC

The transcriptomic factors (TFs) were predicted using single-cell regulatory network inference and clustering (SCENIC) by performed SCENIC R package [20].

Score according to different gene sets

To calculate module scores and the fraction of enrichment for gene expression of specific gene set in single cells, “AddModuleScore” and “AUCell” function were performed. Using the ggstatsplot R package, a violin plot was created to visualize the scoring results. The heatmap R package was utilized to generate a heatmap for visualizing gene expression.

ST data preprocessing, integration and data visualization

We followed the standard Seurat workflow for dimensionality reduction and clustering to create the spatial transcriptomics dataset. The “AddModuleScore” function was employed to score the spatial transcriptomics data based on gene sets representing specific cell subtype characteristics from our annotated single-cell data, and visualization was accomplished using the “SpatialPlot” function. To visualize the expression of specific ligand-receptor pairs in the slices, we used the “plotLR” function from the SpaGene R package.

Cell culture and Arecoline stimulation

THP-1 cells were obtained from the Cell Bank of National Collection of Authenticated Cell Cultures. Culturing THP-1 cells in RPMI-1640 medium with 10% FBS. Prepare a cell suspension with a density of 1 million cells/ml. Add 1 µl of 0.1 mg/ml PMA to each 1 ml of the suspension and inoculate 4 ml of this cell suspension into each 6 mm dish. Incubate in a 37 °C incubator for 24 h, then add different concentrations of arecoline (0, 0.5, 5 µg/ml) for stimulation for 48 h.

Quantitative real-time PCR

Total RNA was extracted using the RNeasy mini kit and was reverse-transcribed to cDNA using the QuantiTect Reverse Transcription Kit. cDNA was then mixed with primers and iQ SYBR Green Supermix in a PCR eight-row tube. qRT-PCR was performed using the iCycler Thermal Cycler (Bio-Rad Laboratories, USA). Relative gene expression levels were calculated via the 2 − ΔΔCT method with GAPDH as the endogenous control. Gene expression data were statistically analyzed using unpaired t-tests in GraphPad Prism 9.0 (GraphPad Software, USA) and visualized through column graphs depicting fold-change values normalized to control groups. The unpaired t-test was applied to compare differences in INHBA and TGFβ expression levels across THP-1 cells treated with varying concentrations of arecoline using GraphPad PRISM (version 8.0; GraphPad Software).

Multiple immunofluorescences staining for clinical samples

Paraffin tissue sections of 4 patients with NODSCC and 4 patients with ODSCC who underwent surgery at Qilu Hospital of Shandong University were selected. All pathological states were confirmed histopathologically by H&E staining. Multiple immunofluorescence (mIF) staining of tissue was performed using Opal Chemistry (PerkinElmer, Waltham, MA, USA). Briefly, the sections were labeled with primary antibodies anti-INHBA (Proteintech, 60352-1-Ig), anti-CD3 (ZSGB-BIO, ZM-0417), anti-CD8 (ZSGB-BIO, ZA-0508), and anti-PD-1 (ZSGB-BIO, ZM-0381), CD4 (Abcam, ab133616), Foxp3 (Abcam, ab20034), followed by HRP-conjugated secondary antibody. Subsequently, the fluorophore-conjugated tyramide amplification system (PerkinElmer) was used for signal amplification, and DAPI was used to counterstain the nuclei. Visualization and quantitation of the different fluorophores were achieved with Tissue FAXS Spectra Systems and Strata Quest analysis software (Tissue Gnostics).

Additional methodological details are provided in the Supplementary Materials (Suppl. Materials. 1).

Result

Global landscape of single-cell transcriptomics and intercellular communication in ODSCC and NODSCC

We conducted single-cell RNA sequencing on the GSE215403 dataset from the GEO database, including samples from 3 ODSCC and 9 NODSCC patients. After data preprocessing using standard quality control metrics (e.g., nCount_RNA < 100000, 7500 > nFeature_RNA > 400), we obtained a total of 30,303 cells, with 9,881 from the ODSCC group and 20,422 from the NODSCC group. Unsupervised clustering with Seurat (resolution = 0.6) identified 15 distinct clusters (Suppl. Figure 1A). t-SNE dimensionality reduction(perplexity = 15) and manual annotation based on classic marker genes categorized these cells into 8 clusters (Fig. 1A): B cells (CD19, CD79A, MS4A1), endothelial cells (PECAM1, VWF), epithelial cells (DSP, KRT18, CDH1, KRT8, EPCAM), fibroblasts (FGF7, MME, ACTA2, DCN, LUM), mast cells (TPSB2, TPSAB1), myeloid cells (C1QA, C1QB, MMP19, FCGR3A, FCN1, S100A12, CD1E, CD1C), plasma cells (IGHG1, MZB1), and T cells (CD3E, CD3D, PTPRC, NKG7) (Fig. 1B). Notable differences in cell type proportions among samples reflect strong tumor heterogeneity (Suppl. Figure 1.B). Ro/e algorithm analysis revealed that in ODSCC the epithelial (Ro/e = 1.28) and plasma cells (Ro/e = 1.39) were enriched, while endothelial cells (Ro/e = 0.44), B cells (Ro/e = 0.58), and mast cells (Ro/e = 0.46) were less prevalent compared to the NODSCC group (Fig. 1C). Cellchat analysis showed a complex intercellular communication, of which fibroblasts had more and stronger communication with other cell types, particularly myeloid cells (communication strength = 2.1, count = 87), endothelial cells (communication strength = 2.6, count = 145), T cells (communication strength = 1.4, count = 38) et al. (Fig. 1D, Suppl. Table 1). Comparative analysis indicated that T cells in the ODSCC group receive more regulatory signals, while epithelial cells had relatively weaker interactions with the other cells, potentially due to obstructive effect of more collagen deposition in submucous fibrosis (Fig. 1E, Suppl. Figure 1C, Suppl. Table 1). Further comparative analysis revealed uniquely activated signaling pathways in the ODSCC group (p < 0.01) (Suppl. Figure 1D), including TWEA (TNFSF12-TNFRSF12A) pathway, involved in myeloid cells activating CAFs and promoting CAF-monocyte interaction [21] and the BMP2-(BMPR1A + BMPR2) receptor pathway, associated with CAF transitioning to a lipid-loaded phenotype, thus promoting tumor metastasis and proliferation [22] (Fig. 1F). Conversely, the ODSCC group lacked an activated FGF pathway (Fig. 1.F), affecting CAF differentiation or interaction with endothelial cells [23]. These findings highlight the distinct intercellular communication, especially between fibroblasts and myeloid cells or T cells, in ODSCC compared to NODSCC.

Fig. 1

Global landscape of single-cell transcriptomics and intercellular communication in ODSCC and NODSCC. (A) t-SNE dimensionality reduction-categorized 8 major cell types (left), 12 samples (middle) and 2 groups (right). (B) Dot plot of the marker genes in each major cell type. (C) Ro/e algorithm analysis of tissue distribution of major cell types in each group. (D) Cellchat analysis of intercellular communication in ODSCC and NODSCC. (E) Chord chart showing cell-cell interaction among the major cell types in ODSCC and NODSCC (up: numbers of interactions, down: strength of interactions). (F) Bubble plot comparing the significant ligand-receptor pairs between ODSCC and NODSCC

The epithelial cells in ODSCC exhibits stronger malignant characteristics

29,261 epithelial cells were classified into six subcluster (C0-C5) (Suppl. Figure 2 A). Analysis of the epithelial cell composition across different samples revealed that clusters C5 (Ro/e = 2.37 vs. 0.02) and C0 (Ro/e = 1.24 vs. 0.83) predominantly presented in the ODSCC group, while clusters C1 (Ro/e = 1.34 vs. 0.53) and C4 (Ro/e = 1.38 vs. 0.47) were primarily found in the NODOSCC (Suppl. Figure 2B). t-SNE (perplexity = 50) was employed to visualize the distribution of benign and malignant cells (Suppl. Figure 2 C) and Copy Number Variation (CNV) analysis was used to differentiate between benign and malignant cells (Suppl. Figure 2D, E). GSVA of Hallmark gene sets across different epithelial subpopulations (Suppl. Figure 2 F) showed that the C4 subcluster, characterized as benign based on CNV analysis, had downregulated proliferation-related pathways, suggesting that it represents normal epithelial cells. The low representation of C4 in ODSCC indicates that even non-malignant epithelium in OSF deviates from the normal epithelial expression profile. In contrast, C5 enriched in pathways related to angiogenesis, cell proliferation, EMT and HIF-A related to hypoxia, representing a unique malignant cell subcluster in ODSCC. The increase in this subpopulation reflects a stronger malignant phenotype and a more pronounced hypoxic microenvironment in ODSCC.

T cells in ODSCC exhibit a more severe immunosuppressive landscape

In the TME, T cells can be activated into effector T cells to kill tumor cells upon antigen stimulation. However, under persistent homologous antigen stimulation, the effector functions and proliferation abilities of T cells will be impaired, a phenomenon known as T cell dysfunction [24]. Analysis of the two groups revealed a decrease in T cells in ODSCC (Fig. 1C). To validate the distribution pattern of T cells, we integrated spatial transcriptomic data from 4 cases of NODSCC (GSE208253) and 4 cases of ODSCC (GSE220978). After scoring and mapping the T cell subpopulations, we found that, compared to NODSCC, T cells in ODSCC predominantly localized to the stromal region, with very few infiltrating the tumor area (Fig. 2A, Suppl. Figure 3 A). The presence of T cell exclusion effects suggests an immune-excluded tumor microenvironment in ODSCC [25, 26]. Via dimensionality reduction clustering, T cells were categorized into eight subclusters: CD4⁺ exhausted T cells, CD4⁺ naive T cells, CD8⁺ exhausted T cells (CD8⁺Tex), CD8⁺ naive T cells (CD8⁺Tn), CTLs, naive T cells, NK cells, and Treg (Fig. 2B, Suppl. Figure 3B). Further analysis of the two groups revealed increased Treg (Ro/e = 1.23), CD8⁺Tex (Ro/e = 1.12) and CD8⁺Tn (Ro/e = 1.32), but a decreased CTLs (Ro/e = 0.91) in ODSCC (Fig. 2C). Further dimensionality reduction clustering of CD8⁺T cells allowed for the re-classification of CD8⁺Tex cells into Tterm (PD1^hiHAVCR2⁺TOX⁺) and Tprog (PD1^intGZMA⁺ITGAE⁺CTLA4⁺) which can be reversed by anti-PD-1 treatment [27] (Fig. 2D, Suppl. Figure 3 C). In ODSCC, Tterm exhibited higher expression of exhaustion markers such as INFG, CXCL13, CCL3, PDCD1, and LAG3 [28] (P = 1.49e-03, 95% CI [-0.53, -0.13], Fig. 2E, F, Suppl. Table 2). Monocle2 analysis of CD8⁺T cell differentiation trajectories confirmed that Tterm represented the terminal differentiation state of CD8⁺T cells, while Tprog represented an intermediate state in the CD8⁺T cell exhaustion process (Suppl. Figure 3D). Although CD8⁺T cells increased in ODSCC, the increase was primarily in CD8⁺Tn (Ro/e = 1.37) and CD8⁺Tex with a notable rise in Tterm (Ro/e = 1.22), while Tprog kept a similar distribution (Ro/e = 1, Fig. 2G). Additionally, exhaustion markers and immune dysfunction-related transcripts such as SOX4, FOXP3 and PRDM1 were significantly upregulated in ODSCC compared to NODSCC (Suppl. Figure 3E, Suppl. Table 2), suggesting a lower sensitivity of ODSCC to anti-PD-1 immunotherapy. Notably, the MHC-I signaling pathway was upregulated in ODSCC, particularly affecting CD8⁺Tex (Fig. 2H, I). This implies that the sustained activation of MHC-I signaling pathway might be one of the reasons for the functional impairment of CD8⁺T cells in ODSCC [29].

To validate the distribution patterns of CD8⁺T cell, we performed CD3, CD8 and PD1 mIF staining on four ODSCC and four NODSCC tissue slices. We found that the proportion of CD8⁺Tex cells (CD3⁺CD8⁺PD1⁺) in ODSCC was significantly higher than in NODSCC (p < 0.05) (Fig. 2J), and this result was revalidated using spatial transcriptomics data (Fig. 2K, Suppl. Figure 3 F), indicating a high exhausted T cell state in ODSCC.

Fig. 2

T cells in ODSCC exhibit a more severe immunosuppressive landscape. (A) Distribution of T cells in ODSCC and NODSCC revealed by spatial transcriptomics data. (B) t-SNE of major subclusters of T cells. (C) Tissue distribution of major T cell subtypes in each group. (D) Further dimensionality reduction and clustering-categorized CD8 + Tex cells into Tterm and Tprog subtypes. (E) Heatmap comparing the expression of marker genes associated with T cell exhaustion in Tterm between ODSCC and NODSCC. (F) Violin plot comparing the expression of marker genes associated with T cell exhaustion in Tterm between ODSCC and NODSCC. (G) Tissue distribution of CD8 + T subsets. (H, I) Heatmap illustrating the differences in MHC-I signaling regulation between ODSCC (H) and NODSCC (I). mIFC demonstrating a higher prevalence of exhausted T cells in ODSCC than NODSCC. (J) mIFC demonstrating a higher prevalence of exhausted T cells in ODSCC than NODSCC. (K) Spatial feature plots of the signatures of CD8 + Tex in ODSCC and NODSCC tissue sections

Tregs, as crucial regulators of the TIME, were more prevalent in ODSCC (Fig. 2C). Immune suppression function analysis showed that the co-inhibitory molecule score of Tregs in ODSCC was higher than in NODSCC (P < 0.001, 95% CI [-0.41, -0.19], Fig. 3A, B, Suppl. Table 2), Further mIF analysis revealed that ODSCC had a higher number of Tregs (p < 0.05 Fig. 3C), further validated via using spatial transcriptomics data (Fig. 3D). These results suggest that the increased proportion and stronger immune-suppressive function of Tregs in ODSCC may be significant factors contributing to the more severe TISME in ODSCC.

Fig. 3

Enrichment of Tregs with stronger immune suppression functions in ODSCC. (A) Heatmap showing a higher co-inhibitory molecule expression of Tregs in ODSCC than in NODSCC. (B) Violin plot showing a higher co-inhibitory molecule score of Tregs in ODSCC than in NODSCC. (C) mIFC showing increased enrichment of Tregs in ODSCC. (D) Spatial feature plots of the signatures of Treg in 4 ODSCC and 4 NODSCC tissue sections

Overall, these findings demonstrate that total T cell infiltration levels are lower in ODSCC, yet, the distribution of T cell subgroups shows a more severe immune-suppressive landscape compared to NODSCC. This may predict a poorer response to anti-PD-1 immunotherapy and potentially worse prognosis for ODSCC.

INHBA⁺Mac regulates the TIME in ODSCC and is associated with a poorer prognosis

Macrophages play a central role in immune regulation in TME [30]. To explore the impact of macrophage in the TIME in ODSCC, macrophages were divided into 6 subclusters. (Fig. 4A). Based on the top 10 marker gene (Suppl. Figure 4A) and GSVA functional enrichment analysis (Suppl. Figure 4B), the six subclusters were respectively defined as IDO1⁺Mac, CCL18⁺Mac, CCL2⁺Mac, S100A2⁺Mac, CXCL10⁺Mac, and INHBA⁺Mac. IDO1⁺Mac was characterized by high expression of CLEC10A, AREG and IDO1, immune evasion-related genes [31,32,33], and enriched in the Th17 differentiation pathway. CCL18⁺Mac exhibited high expression of CCL18, APOE and SLC40A1, which is associated with immunosuppression, pro-cancer and tumor cell metabolism [34,35,36], and enriched in lipid metabolism pathways. CCL2⁺Mac was characterized by the expression of chemokines such as CCL2, CCL8, and CXCL1 and primarily enriched in chemokine pathways related to immune cell infiltration. S100A2⁺Mac showed high expression of genes related to angiogenesis, such as S100A2and LGALS3 [37, 38], and enriched in mucosal innate response and vascular endothelial growth factor-related pathways. CXCL10⁺Mac was marked by high expression of genes such as TNFSF10, LGALS2, and mainly enriched in pathways related to B cell proliferation and immune suppression. INHBA⁺Mac exhibited high expression of matrix remodeling genes e.g. TNFAIP6, SERPINB2 and MMP1 [39,40,41] and primarily enriched in pathways related to angiogenesis.

INHBA⁺Mac was more prevalent in ODSCC (Ro/e = 1.3) compared to NODSCC (Ro/e = 0.88), making it the most significantly different subcluster between the two groups(Fig. 4B). In the TCGA HNSCC data, analysis using deconvolution methods revealed that INHBA⁺Mac was significantly associated with poorer prognosis (HR = 1.37, 95% CI [1.05, 1.8], P = 0.0219, Fig. 4C). The relationship between the top 10 marker genes of INHBA⁺Mac and HNSCC prognosis in TCGA database was evaluated using GEPIA2, revealing that INHBA is most significantly negatively correlated with prognosis (Logrank p = 0.0011, p(HR) = 0.0012, Suppl. Figure 4C). The common immunosuppressive molecule SPP1 is also highly expressed in INHBA⁺Mac (p < 0.0001, FDR < 0.0001) [42] (Suppl. Figure 5 A, Suppl. Table 3). Immune checkpoint scoring showed the highest score for INHBA⁺Mac among macrophages (P = 1.62e-37, 95% CI [0.16, 1.00], Fig. 4D, Suppl. Table 2). Moreover, INHBA⁺Mac in ODSCC significantly increased expression of immunosuppressive molecules such as CD274/PD-L1, ADORA2A, and PVR (Fig. 4E, Suppl. Table 2), while specifically reduced co-stimulatory molecules like CD86, CD40 and TNFSF8 (Suppl. Figure 5B) compared with that in NODSCC. MHC sensitivity scores related to immune therapy showed that ODSCC had a lower score overall versus NODSCC (p < 0.001) (Suppl. Figure 5 C, D), with INHBA⁺Mac in particular having the lowest MHC score among macrophage subclusters [43] (P = 7.95e-21, 95% CI [0.08, 1.00], Fig. 4F, Suppl. Figure 5E, Suppl. Table 2). To show the modulation of INHBA⁺Mac on Treg, we used CIBERSORT to deconvolute the expression matrices of macrophage subclusters, and then analyzed the correlation between INHBA⁺Mac and Treg. We found that INHBA⁺Mac was significantly positively correlated with Treg enrichment (Rs = 0.18, P = 2.4e-03; Suppl. Figure 5 F). mIF further validated distribution correlation of Treg and INHBA⁺Mac (p < 0.05 Fig. 4G).

All the above results suggest that INHBA⁺Mac is a main subset in ODSCC among macrophages and more prevalent than in NODSCC. INHBA⁺Mac in ODSCC exhibits more pronounced immunosuppressive functions and lower sensitivity to immune therapy than that in NODSCC.

Fig. 4

INHBA + Mac regulates the TIME in ODSCC. (A) t-SNE plot showing the classification of macrophage subclusters. (B) Tissue distribution of macrophage subclusters in each group. (C) The relationship between the INHBA⁺Mac subcluster and patient overall survival (OS) in the TCGA database. (D) Violin plot showing the immunosuppression scoring of macrophage subsets. (E) Heatmap showing the expression of immunosuppressive ligand and receptor molecules in macrophage subclusters of ODSCC and NODSCC. (F) Violin plot showing reduced MHC scoring for INHBA⁺Mac among macrophage subclusters. (G) mIFC-confirmed distribution correlation of Tregs (CD4 + FOXP3+) and INHBA⁺Mac (INHBA + CD68+)

The pro-cancer and immunosuppressive functions of iCAF in ODSCC

CAFs are crucial regulators in the TME, particularly in cancer cell proliferation and invasion, neovascularization, inflammation, extracellular matrix (ECM) remodeling [44] and immunosuppression [45]. To explore the distinct subset of CAFs and its special roles in ODSCC, fibroblasts were re-clustered into six distinct clusters (0–5) based on high-variance genes. Cluster 0, characterized by elevated expression of cytokines and chemokines IL6, IL11, CXCL1 and CXCL8, was identified as iCAF. Clusters 1, 2, and 4 were classified as myCAF, marked by ACTA2, MYL9, and MYLK. Cluster3 represented mCAF, identified by POSTN, COL1A1, COL1A2 and COMP. Cluster5 was classified as apCAF, defined by HLA-DRB1, HLA-DRA, and CD74 (Fig. 5A, B, Suppl. Figure 6 A). Analysis of CAF distributions revealed that iCAF level was significantly higher (Ro/e = 1.25), while apCAF level lower (Re/o = 0.5) in ODSCC compared to NODSCC (Fig. 5C). GSVA enrichment analysis indicated that iCAF mainly enriched in immune-related pathways and collagen/ECM pathways, correlating with collagen deposition in the matrix (Fig. 5D). Further analysis of the differential iCAF gene expression between the two groups revealed that genes involved in collagen metabolism and promoting tumor cell invasion and proliferation, such as WNT5A (p < 0.001, FDR < 0.001), CTSK (p < 0.001, FDR < 0.001), LUM (p < 0.001, FDR < 0.001), DCN (p < 0.001, FDR < 0.001) and COL7A1 (p < 0.001, FDR < 0.001) [46, 47], associated with immune-suppression, e.g. TDO2 (p < 0.001, FDR < 0.001) [33, 48], and the T-cell activation pathway inhibitor DUSP4 (p < 0.001, FDR < 0.001) [49], were upregulated in iCAF from the ODSCC group (Fig. 5E, Suppl. Figure 6B, Suppl. Table 4). Further Cellchat analysis revealed that iCAF exerted intimate communication with T cells and macrophages via immunosuppressive receptor-ligand pairs (Fig. 5F, Suppl. Figure 6 C). Pseudotime analysis of CAFs subclusters revealed that iCAF were the starting point of the CAFs pseudotime trajectory, with the endpoint being apCAF and myCAF (Fig. 5G). This further shows out the importance of iCAF, in that myCAF enrich in ACTA2, S100A2 and RGS5 and have been demonstrated to exert pro-cancer [50,51,52] and immune-suppression [53, 54].

All these results indicate that iCAF exhibit stronger pro-cancer and immune-suppressive functions in ODSCC than in OSCC.

Fig. 5

iCAF significantly increased and possessed stronger immune-suppressive functions in ODSCC. (A) t-SNE dimensionality reduction-categorized 4 subtypes of CAFs. (B) Heatmap showing TOP 5 markers in indicated CAF subclusters (C) Tissue distribution of CAF subclusters in each group. (D) Bubble plot of representative GO pathways enrichment in CAF subclusters predicted by GSVA. (E) Volcano plot of differentially expressed genes between NODSCC and ODSCC groups in CAFs. (F) The regulation of T cells and macrophages by iCAF through immune-suppressive signals. (G) Pseudotime analysis of CAFs (cell types, pseudotime)

INHBA is involved in the modulation effect of INHBA⁺Mac and iCAF on Treg through the SMAD pathway in ODSCC

The above presented results clearly show that both INHBA⁺Mac and iCAF were associated with the TISME formation whereas INHBA⁺Mac has the highest INHBA expression among macrophage subsets (Suppl. Figure 4 A) and iCAF, the predominant CAF subtype in ODSCC, also highly express INHBA (Fig. 5E, Suppl. Figure 6B). This pushed us to find the key role of INHBA in inducing TISME. INHBA has been revealed to induce Foxp3 expression and Treg generation by activating SMAD2/3 phosphorylation [55]. Here, gene differential expression analysis revealed that INHBA presented the highest expression in myeloid cells, followed by CAFs and epithelial cells and INHBA expression in myeloid cells and CAFs from ODSCC was higher than from NODSCC (p < 0.001, FDR < 0.001, Fig. 6.A, Suppl. Figure 7 A, Suppl. Table 5). Analysis of INHBA expression in four samples of ODSCC from the GSE220978 spatial transcriptomics data revealed that INHBA expression was primarily concentrated in the tumor region and in the OSF regions (Suppl. Figure 7B), mIF analysis in clinical tissues revealed that the expression of INHBA is higher in ODSCC than in NODSCC (P < 0.01, Fig. 6B). Interestingly, arecoline, a primary alkaloid found in betel nuts and a classic inducer for OSF, significantly increased the mRNA expression of INHBA and TGFβ in in vitro cultured THP-1-derived macrophages (p < 0.001), with a more pronounced increase in INHBA mRNA expression (100 to 10,000 times) (Fig. 6C).

To investigate if INHBA is involved in the modulation effect of INHBA⁺Mac and iCAF on Treg, spatial transcriptomic analysis was performed to observe spatial organization of INHBA⁺Mac, iCAF and Treg in spatial transcriptomics slices. Results revealed a close colocation of three types of cells and more interestingly, using SpaGene we identified the ligand INHBA-receptor ACVR1/ACVR2A/ACVR2B interaction regions overlapping with distribution of three types of cells (Fig. 6D, Suppl. Figure 7B). Moreover, our further validation using mIF revealed that INHBA is in close proximity to Treg, and this finding is statistically significant (p < 0.01) (Fig. 6E, F). These results suggest that INHBA is potentially related to the regulatory effect of INHBA⁺Mac and iCAF on Treg.

TGFβ/SMAD activation plays a crucial role in Treg [56] and TGFβ1 is the most representative subtype among three isoforms of TGFβ in tumor immune suppression [57]. INHBA, a member of the TGFβ superfamily, shares similar structure and canonical SMAD2/3 pathway to TGFβ and when TGFβ signaling is compromised, INHBA can compensate for the deficiency in SMAD2/3 phosphorylation [55]. GSEA enrichment analysis of CD4⁺Tn and Treg revealed that SMAD-related pathways were enriched in Treg (Fig. 6G). This underscores the crucial role of the TGFβ/SMAD pathway in the activation of Treg in OSCC. To analyze whether the effect of INHBA on Treg is related to ActivinRI/II-SMAD signaling pathway, the activated TGFβ superfamily members acting through SMAD2/3 pathway were compared between ODSCC and NODSCC. Results discovered that INHBA was most obviously upregulated among TGFβ superfamily members in ODSCC versus NODSCC (Fig. 6H, Suppl. Table 5). Correspondingly, the moderate expression level of ActivinRI/II was found in Treg in both ODSCC and NODSCC (Fig. 6I). However, TGFβ1 was downregulated (Fig. 6F), whereas GSEA enrichment analysis showed no significant difference in SMAD2/3-related signaling activation (Suppl. Figure 8 A), indicating that INHBA compensates for the insufficient TGFβ1 expression to activate SMAD2/3 pathway.

In summary, these results highlights the critical role of INHBA in activating the SMAD pathway to promote Treg formation in ODSCC.

Fig. 6

INHBA is involved in the modulation effect of INHBA⁺Mac and iCAF on Treg through the SMAD pathway in ODSCC. (A) Mountain plot of the expression of INHBA in different cell subsets in ODSCC and NODSCC.(B) mIFC demonstrating a higher expression of INHBA in ODSCC than NODSCC. (C) Bar chart of the relationship between different concentrations of arecoline stimulation and TGF(Left)/ INHBA (Right) expression in macrophages. (D) Spatial transcriptomics analysis of the overlapping area of INHBA-ACVR1/ACVR2A/ACVR2B and iCAF/Tregs/INHBA⁺Mac. (E) mIFC confirmed that the distribution of Tregs (CD4 + FOXP3+) is correlated with INHBA (left). (F) Statistical graph showing the spatial distribution of Treg around INHBA. (G) GSEA analysis of SMAD2/3-related signaling pathways in Tregs. (H) Mountain plot of the expression distribution of TGFβ superfamily members in two groups. (I) Heatmap of the expression distribution of the receptor of INHBA in different cell types of CD4 + T cells

Discussion

Previous studies have identified clinical and pathological differences between ODSCC and NODSCC [5, 6]. However, the differential features of TME, especially TIME between the both are not fully understood. Here, we conducted an in-depth analysis of the data, focusing on differences in transcriptomic profiles between ODSCC and NODSCC. Our analysis revealed a higher proportion of tumor epithelial cells, a reduced presence of stromal components such as B cells and endothelial cells, and a more pronounced TISME in ODSCC. Especially, we propose for the first time that increased proportion and immune suppression activity of INHBA⁺Mac and iCAF are characteristics of ODSCC. INHBA is involved in the modulation effect of INHBA⁺Mac and iCAF on Treg (see Fig. 7).

Fig. 7

Graphical illustration of the working model. In ODSCC, upregulated INHBA mediates crosstalk between INHBA⁺Mac and iCAF via the INHBA-ACVR1/ACVR2A/ACVR2B ligand-receptor axis, activating the SMAD signaling pathway to induce Treg differentiation and functionally exert immunosuppressive activity

In the TIME, cytotoxic T cells are commonly manifested by dysfunction, presenting an exhausting status and the number and immune suppression function of Treg are generally enhanced [57]. In present study, compared to NODSCC, ODSCC showed a significant increase in CD8⁺Tn, CD8⁺Tex and Treg cells, while CTLs were reduced. These findings suggest a potentially more severe immune-suppressive landscape in ODSCC than in NODSCC.

Immune checkpoint blockade (ICB) aims to boost CD8⁺T cell responses against cancer [24] and different exhausting status of CD8⁺T cells exerts varied influence on ICB sensitivity [27]. Exhausted CD8⁺T cells contain a subset of progenitor exhausted and terminally differentiated T cells (Tterm) that differentiated from the former. Progenitor exhausted T cell can well respond to anti-PD-1 therapy, but Tterm cannot [27]. Our analysis reveals a marked increase in Tterm in ODSCC (Ro/e = 1.22). Moreover, Tterm (PD1^hiHAVCR2⁺TOX⁺) in ODSCC exhibited higher expression of exhaustion markers (p < 0.001) such as INFG, CXCL13, CCL3, PDCD1, and LAG3 than that in NODSCC. Additionally, exhaustion markers and immune dysfunction-related transcripts such as SOX4, FOXP3 and PRDM1 were significantly upregulated in ODSCC compared to NODSCC. All these results suggest that ODSCC may be less amenable to reversal by anti-PD1 immunotherapy than NODSCC.

Macrophages are central regulators in the TIME, while M2 is one of the key immune suppressive cells [30]. With the advancement of single-cell sequencing technology, more precise classification methods based on gene expression profiles have gradually replaced the traditional M1/M2 macrophage classification [34] and more macrophage subsets have been identified, including TREM2⁺TAMs related to immunosuppressive TME [58], the SPP1⁺TAMs associated with angiogenesis in colon cancer [59], and PD-L1⁺macrophage mediating immune evasion in melanoma [60]. Our previous study also found that Macro-IDO1 was main macrophage subset in oral leukoplakia-derived OSCC and had a strong immunosuppressive role and contributed to oral carcinogenesis [33]. However, pro-cancer activity of INHBA + monocytes/macrophages has been seldom evaluated until recently to our limited knowledge [61, 62]. The present study novelly found that INHBA⁺Mac was more prevalent in ODSCC (Ro/e = 1.30) compared to NODSCC (Ro/e = 0.88), making it the most significantly different subcluster between the two groups.

In addition to INHBA high expression, INHBA⁺Mac highly expressed common immunosuppressive molecule SPP1 (p < 0.001), and possessed strong immune-suppressive functions, evidenced by higher immune checkpoint scores, increased expression of immunosuppressive molecules such as CD274, ADORA2A, and PVR, and having the lowest MHC score among macrophage subclusters. This implies that the increase in INHBA⁺macrophages is potentially responsible for the stronger immunosuppression in ODSCC.

CAFs, another key component of the TME that exert a crucial effect on TIME [44, 45], present complex heterogeneity [33, 48]. Based on highly variable genes and functional enrichment, CAFs were classified into four subclusters (iCAF, mCAF, myCAF, apCAF) in this study. We unveiled that iCAF, characterized by elevated expression of cytokines and chemokines IL6, IL11, CXCL1 and CXCL8, represented a high proportion among the identified CAF subsets and was obviously enriched (Ro/e = 1.25), while apCAF depleted in ODSCC (Ro/e = 0.50) compared to NODSCC (Ro/e = 1.23). GSVA enrichment analysis indicated that iCAF mainly enriched in immune-related pathways and collagen/ECM pathways. More importantly, further analysis demonstrated that iCAF in ODSCC possessed stronger immune-suppressive functions than those in NODSCC, as shown by differential immune-suppression gene expression, e.g. TDO2 and IDO1, between the two groups, the upregulated T-cell activation pathway inhibitor DUSP4 in the ODSCC group. apCAF, defined by HLA-DRB1, HLA-DRA, and CD74, can activate T cells and induce tumor suppression [63]. The present analysis showed that apCAF mainly enriched in pro-angiogenic and antigen-presenting functions. Thus, the reduction in apCAF may be partial causes of the immunosuppressive microenvironment exacerbation and the reduced number of blood and lymphatic vessels in ODSCC relatively to NODSCC. Anyway, this study suggests that the increase in iCAF and the decrease in apCAF may be another distinctive TIME landscape of ODSCC from NODSCC.

INHBA is a member of the TGFβ superfamily and has been reported to promote the formation of a TISME [64], promoting lapatinib resistance [65]. In HNSCC, INHBA is expressed at higher levels in tumors compared to normal tissues, and its overexpression is associated with a poor prognosis [66, 67]. The classical TGFβ signaling pathway regulates tumor immunity through the activation of the SMAD pathway [57]. Recent clinical trials have validated the role of TGFβ-targeted drugs in inhibiting Treg production and enhancing the cytotoxicity of CD8⁺T cells [68]. Targeting TGFβ signal (anti-TGFβ antibodies, TβR inhibitors) can synergistically enhance the effects of other immunotherapy approaches [69]. Similar to TGFβ, INHBA can activate the downstream SMAD pathways through its receptor, exhibiting similar functions to TGFβ and compensating for the functional defects caused by TGFβ molecule deficiency [55]. Interestingly, significant overexpression of INHBA in iCAF was observed. Nagaraja et al. show that inhibin β A is an important regulator of the CAF phenotype in ovarian cancer [70]. Hu et al. identified an INHBA(+) subset of immunomodulatory pro-tumoral CAFs as a potential therapeutic target in advanced ovarian cancers which typically show a poor response to immunotherapy [55]. The bioinformatics by Zheng et al. demonstrated that CAFs producing INHBA promotes colorectal cancer development and correlates with poor prognosis [14]. Yu et al. by bioinformatics reveal that INHBA expression strongly correlated with various markers of monocytes/macrophages and cancer-associated fibroblasts in breast cancer [71]. Our study revealed that both INHBA⁺Mac and iCAF are main origins of INHBA in ODSCC, both of which are related to the formation of a TISME, in keeping with the previous studies [55, 61, 62]. Moreover, in ODSCC, although TGFβ1 expression is relatively low (p < 0.001), INHBA expression is elevated (p < 0.01) and the degree of SMAD pathway activation remains similar between ODSCC and NODSCC (p > 0.05). This further emphasizes the role of INHBA in OSCC, particularly in ODSCC, where it may partially substitute for TGFβ, thereby activating SMAD and downstream pathways, affecting the tumor microenvironment and patient survival. Analysis in the TCGA database showed that, similar to TGFβ, INHBA expression has a significant negative correlation with prognosis (p = 0.0012). All these results suggest that INHBA is a distinctive immunosuppressive molecule, highlighting the potential of INHBA as a therapeutic target, especially in ODSCC.

In summary, we utilized scRNA-seq and ST data from the GEO database and experimental validation to reveal the distinctive TIME landscape of ODSCC, a subtype of OSCC with relatively poor prognosis [5]. Our results suggest that compared to NODSCC, ODSCC shows a more severe TISME and the poorer sensitivity to immunotherapy. The increased INHBA⁺Mac and iCAF seem to be responsible for these immune characteristics in ODSCC. The upregulated INHBA in ODSCC and INHBA-ACVR1/ACVR2A/ACVR2B interaction may mediate the modulation effect of INHBA⁺Mac and iCAF on Treg differentiation and functionality. This underscores the therapeutic potential of INHBA and provides a theoretical basis for developing personalized treatment plans for OSCC. However, there are some limitations in our study. Firstly, the analysis of ODSCC data is based on only 3 cases, which may lead to internal errors and less generalizability of sequencing results. Secondly, single-cell RNA sequencing, while powerful, has inherent limitations in detecting low-abundance transcripts and resolving rare cell subsets. Manual cell annotation based on marker genes may also introduce subjectivity. Thirdly, our proposed INHBA-SMAD-Treg axis is primarily supported by spatial co-localization and pathway enrichment analyses. Direct experimental validation (e.g., SMAD2/3 phosphorylation assays or Treg differentiation assays with INHBA blockade) is needed to establish causality. Fourthly, another important issue is the absence of a stable and reliable ODSCC tumor model, which prevented us from validating our findings in vivo. The recent emergence and maturation of organoid models may partially replace animal models and further validate our discoveries in future.