The global landscape of single-cell primary and recurrent OTSCC
We performed single-cell and TCR sequencing of fresh human tissue samples of PT, paired with non-tumor tongue tissue and RT, in order to elucidate the dynamic landscape and tumor heterogeneity in the development of OTSCC chemoresistance and immune inhibition subsequent to chemotherapy. Consequently, single-cell 5’ gene expression and V(D)J libraries were constructed on a 10X platform (10X genomics, California, United States). A patient with four years of carboplatin and DTX combination chemotherapy had recurrent OTSCC removed. In addition, 12 additional single-cell datasets from primary OTSCC were obtained from the GEO database (GSE103322 and GSE172577) in order to better comprehend the heterogeneity of OTSCC tumors. They were contrasted to investigate in depth the transformation of the tumor microenvironment in response to chemoresistance. Through stringent quality filtration and doublet elimination, 62,375 cells and 39,156 genes per cell were retained for further analysis (Fig. 1a). Table 1 lists the patient information and quality of cells for each sample.
To discover distinct cell composition based on the gene expression pattern, we performed dimensional reduction, unsupervised cell clustering, and batch effect removal among multiple samples using Harmony and SCTransform software from the Seurat suite [28, 29]. The selected predominant ten genes for each cluster were list in the additional file 1: Table S1. Some canonical marker genes expression was used to manually annotate the cell clusters, and three major cell types, including epithelial cells (EPCAM), stromal cells (COL1A1, DCN), and immune cells (PTPRC), were identified and visualized using t-SNE (Fig. 1b). As anticipated, nearly all immune and stromal cells congregated together in the tissues obtained from various patients. In contrast, there was greater heterogeneity among the epithelial cells of the patients, specifically when comparing RT, PT, and NT (Figs. 1c, e, Additional file 2: Figures S1a-b). Consistent with previous studies [6, 19], a comparison of cell numbers among the three groups revealed a significantly higher proportion of immune cells in RT than in PT and NT (Fig. 1c), which was supported by the number of immune cells per 3000 cells, which was 2522 in RT compared to 747 in PT and 691 in NT (Addition file 2: Figure S1c). Furthermore, immune cells were clustered into T lymphocyte (T_cell, CD3D, CD3E), B lymphocyte/naive B lymphocyte (B_cell/Naive B_cell, CD79A, MS4A1), plasma cells (Plasma_cell, MZB1), monocyte (CD14, LYZ), mast cell (Mast_cell, CPA3), dendritic cells (DC, CD1C), macrophage (LYZ, CSF1R, C1QC, SPP1) (Fig. 1f). The proportion of immune cell subtypes varied among three groups, particularly macrophage (3,899 in RT vs. 1,044 in PT per 6,000 cells) and T cells (1,212 in RT vs. 2,890 in PT per 6,000 cells). T cells represented the largest proportion of immune cells in normal tissue (76.8%), while macrophage represented only 7%. (Fig. 1g, Additional file 2: Figure S1e). In addition, as previously reported, epithelial, endothelial, and fibroblast cells (42,743 cells) comprised non-immune cells. The selected predominant ten differently expressed genes for 16 clusters were list in additional file 3: Table S2. Furthermore, the majority of non-immune cells in PT and RT were inferred to be cancer-related using NT cells as a reference (Fig. 1h-j). However, for PT and RT, the cell composition was highly differentiation between tumor samples, indicating the chemotherapy-induced changes in tumor ecosystems, particularly the tumor immune microenvironment (TIM), and indicating the need to further decipher the underlying differences between the tumor environments of primary and recurrent OTSCC.
Combined chemotherapy expanded tissue-resident macrophage, reprogrammed TAMs, and promoted immunosuppression in recurrent OTSCC
From 19,632 immune cells, myeloid cells were grouped into monocyte, macrophage, DC, and mast cells. As depicted by t-SNE, monocyte is closely associated with macrophage, which comprises the majority of immune cells in RT (Figs. 1f-g). Therefore, monocytes and macrophages were isolated to investigate their correlation with chemotherapy resistance and their biological role in tumor chemoresistance. Approximately 60% of extracted cells originated from RT (4,822 cells), 38.87% from PT (3,112 cells), and only 97 cells from NT (1.2%), indicating that macrophages are significantly expanded during tumor initiation and progression. In addition, 8,031 cells were aggregated into 18 clusters and manually classified into seven subtypes based on canonical marker genes and the ten most significant DEGs (Figs. 2a-b, Additional file 4: Table S3). We identified four SPP1 + TAM clusters (C1, C2, C3, C17) marked by SPP1 and FN1, five C1Qs + TAM clusters (C4, C8, C11, C13, C16) with high expression of C1QA, C1QB, C1QC, APOE, Etc.), four HLAs + Macro clusters with high expression of HLAs related genes and CCL3L1 (C5, C6, C7, C9), and two non-classic macrophages, without any expression of canonical marker genes (C12, C18). C10, C14, and C15 clusters were labeled as CD14 + Macro, CD68 + Macro, and SELENOP + Macro, indicating that they were highly expressed with CD14, CD68, and SELENOP, respectively (Fig. 2b). All clusters existed in PT and RT, whereas in NT, eight clusters (C4, C5, C7, C8, C9, C13, C14, C17, C18) disappeared. Particularly, SPP1 + TAM accumulated in RT, whereas PT and NT demonstrated a significant decline (3793:336:5) (Fig. 2c). In contrast, there were significantly fewer C1Qs + TAM cells in NT (52 cells) than in PT (1,060 cells) and RT (355 cells) (Fig. 2c). The distribution of SPP1 + TAM and CAQs + TAM in RT and PT suggested that these two macrophage subpopulations may have a significantly relationship with tumor origins, progression, and chemoresistance. Intriguingly, the matrix remodeling-related genes SPP1 and FN1 were substantially up-regulated in RT. More than 93.2% of RT cells expressed SPP1 and 66.3% of RT cells expressed FN1, whereas only 39% and 23% of PT cells expressed SPP1 and FN1, respectively (Fisher exact test, p < 2.22e-16) (Fig. 2d, Additional file 2: Figure S1g). To explore the transformation of the biological function of variable genes in macrophages from normal tongue tissue, primary OTSCC, and recurrent OTSCC, GSVA was performed based on the 50 hallmark pathways, suggesting that gene sets in RT were mostly enriched in EPITHELIAL MESENCHYMAL TRANSITION, NOTCH SIGNALING, MYC TARGETs and WNT BETA CATENIN SINGNALING (Fig. 2e), while for PT, genes were most enriched in regulation of DNA REPAIR, HYPOXIA, P53 pathways (Fig. 2e). The gene ontology enrichment (GO) analysis revealed that the biological process of macrophages in RT was associated with cytokine-mediated signaling pathway and T cell differentiation, in contrast to macrophages in PT and NT, which were associated with positive regulation of response to external stimulus and regulation of cell-cell adhesion, respectively. (Additional file 5: Figure S2 a-c)
In a previous study, it was demonstrated that SPP1 + TAM are associated with tumor angiogenesis and facilitate immune evasion by up-regulating PD-L1 [18, 22]. Consistently, SPP1 + TAM comprised the majority of macrophages in recurrent tumor tissue (Fig. 2f) and also played a crucial role in tumor angiogenesis and epithelial-mesenchymal transition (Fig. 2g). When scoring the pathway activation by calculating the AUC of macrophage, the results indicated that EPITHELIAL MESENCHYMAL TRANSFORMATION and ANGIOGENESIS were powerfully activated in OTSCC after chemotherapy (Figs. 2h-i).
M0 (non-polarized or neutral), M1 (pro-inflammatory and anti-tumor), and M2 (anti-inflammatory or pro-tumor) are the three traditional subtypes of macrophages. We calculated M0, M1, and M2 signature scores with CIBERSORTx (http://cibersort.stanford.edu/) in order to better characterize the phenotype of macrophage subsets. Clusters of SPP1 + TAM and CD68 + Macro displayed an M0 phenotype, SELENOP + Macro and one unique macrophage displayed a more robust M2 signature, and HLAs + Macro displayed a more robust M1 signature (Fig. 2j). Intriguingly, C1Qs + TAM expressed the M1 and M2 signatures, indicating that macrophages in primary OTSCC were polarized, but that polarization vanished following chemotherapy. To better understand the transformation of macrophages from primary OTSCC to recurrent OTSCC, pseudotime trajectory analyses revealed that macrophages in RT and PT were completely separated into two stages, with M0 macrophages and M2 macrophages predominant in the early stages of the RT branch. However, the other branch of PT was enriched in polarized M1 macrophages, supporting the findings that macrophage in primary OTSCC support anti-tumor activity, whereas the biological function of macrophages became pro-tumor and pro-angiogenesis after chemotherapy (Fig. 2k-m). In addition, based on DEGs and TF analyses, despite the C1Qs-related genes (C1QA, C1QB, C1QC), the C1Qs + TAM subtype was confirmed to be highly expressed cathepsin family genes, such as CTSL, CTSZ, etc., which have been linked to the immune response [37]. These genes in C1Qs + TAM were probably affected by the DNA repair regulation genes RAD21 and BCL11A. Additionally, additional MMP-related genes are expressed in the SPP1 + TAM cluster and are likely regulated by POLE3. (Additional file 5: Figure S2 d,e).
PT and RT T cells have a distinct transcriptome and cellular composition
As there were more T cells in the primary OTSCC than in the recurrent OTSCC, a total of 7,957 T cells were extracted for unsupervised clustering, and 16 distinct clusters representing different cell types were identified, including 3 clusters of CD4 + T cells, 8 clusters for CD8 + T cells, and 3 clusters for double-negative T cells with no expression of CD4 and CD8 (Fig. 3a). By evaluating the expression of conventional marker genes, cell categories were manually annotated (Fig. 3b). Regarding CD4 + T cells, two clusters (C4, C9) were categorized as regulatory CD4 + T cells due to the presence of regulatory T cell markers FOXP3, CTLA4, TIGIT, and IL2RA; cluster C16 was identified as exhausted CD4 + T cells due to the high expression of exhausted T cells markers ENTPD1 and HAVCR2. For CD8 + T cells, 5 clusters (C3, C5, C6, C8, C10) were annotated as cytotoxic CD8 + T cells with high expression of cytotoxic markers GZMA, GZMB, PRF1, NKG7, CCL5; C2, C12, and C14 clusters were annotated as transitional CD8 + T cells with high expression of transitional marker genes GZMK; C15 cluster was annotated as naive CD8 + T cell with high expression of naive T cell marker genes SELL, TCF7. The remaining clusters (C1, C7, C11, C13) could not be distinguished by the double-negative T cells annotated as CD8 and CD4. The composition of each cluster of T cells suggested heterogeneity in the immune microenvironment between groups (Fig. 3c). To avoid the bias error generated by the difference in sample size and increase accuracy, we selected 1,326 T cells at random (Counts of T cells in RT) and summarized the subtype composition in RT and PT (Fig. 3d). RT contained a greater number of regulatory T cells than PT.
In contrast, cytotoxic T cells were reduced by more than 50%, indicating that immune responses were significantly different after chemotherapy and that immune inhibition occurred in recurrent tumors. To further investigate the underlying molecular mechanism, we dissected the various expression genes and their enrichment function for each subtype of PT and RT. ECM_RECEPTOR_INTERACTION and UBIQUITI_MEDIATED_PROTEOLYSIS were significantly up-regulated in RT (p < 0.05) according to the GSVA enrichment analysis (Fig. 3e).
TCR sequencing of T cells in recurrent OTSCC, primary OTSCC, and normal tongue tissue
To obtain insight into T lymphocyte alteration and clonal expansion of T cells in RT, we analyzed the single-cell transcriptome and TCR profiling of T cells from recurrent OTSCC (1,604 cells), primary OTSCC (881 cells) and normal tongue tissue (795 cells). All T cells (3,280 cells) from three samples were reintegrated and unsupervised clustering was conducted on them. 4 clusters (C1, C2, C10, C13) for CD8 + T cells, 2 clusters (C4, C7) for CD4 + T cells, 2 clusters (C9, C11) for NK T cells, and 5 clusters (C3, C5, C6, C8, C12) for double-negative T cells were were manually annotated using the same gene expression profiling as discribed above (Fig. 4a,b). Figures 4C-D highlight four main cell types, including CD8 + T cells (CD8T), CD4 + T cells (CD4T), NK T cells (NKT), and double-negative T cells (DNT). To better comprehend the subtype of the T lymphocyte variety, the proportion of all subtypes in three samples. In RT, the proportion of CD8T cells was 29.9%, which was greater than the proportion of CD4T cells (21.4%), whereas in PT, the proportion of cells was 40.6%, which was greater than the proportion of CD4T cells (17.5%).
To gain a clearer understanding of the clonal expansion of T cells in RT, we compared TCR-seq data from RT, PT, and NT. Here, cells with identical CDR3 sequences for the TCR ɑ- and β-chain were categorized as belonging to the same clonotype. From single-cell TCR sequencing data, we detected 3,017 cells, forming 2,177 unique clonotypes, of which 296 contained at least two cells, indicating the clonal proliferation of T cells (Additional file 6: Table S4). Figure 4e demonstrates that the T cell diversity of RT was substantially greater than that of PT and NT. Figure 4f demonstrates that the number of clonotypes with the same clone size increased significantly in RT relative to PT and NT. These results indicate that T-cell clonal expansion is present in RT.
In each cluster, clonally expanded T cells were extensively dispersed, particularly among CD8 + and CD4 + T cells (Fig. 4g). Figure 4H displays the T cell distribution by clone size (= 1, ≥ 2, ≥10, ≥ 20) for each sample. The proportion of T cells with clone sizes ≥ 20 increased considerably in PT and RT relative to NT. The composition of cell types in each sample varied based on clone size (Fig. 4j). Larger clonotype displayed a nonuniform distribution of cell types with an enrichment of CD8 + and CD4 + T cells, whereas only a small number of NK T cells (6%) were detected with ɑβTCR and the clone size was less than three (Fig. 4j). Nonetheless, cytotoxic CD8 + T cells with large clone sizes (≥ 10 and ≥ 20) predominated in RT and PT, and the proportion of cells increased in RT along with clone expansion, indicating that cytotoxic CD8 + T cells promoted tumor chemoresistance in the OTSCC immune microenvironment after chemotherapy.
Clonal association between CD4 + and CD8 + T cells generates a gradient of transcriptional states in RT
We conducted a comprehensive analysis of all subtypes of T cells in three samples. Intriguingly, CD8 + and CD4 + T cells displayed a nonuniform distribution of function stages, with a significant enrichment of cytotoxic CD8 + T cells (C2), regulatory CD8 + T cell (C10) and regulatory CD4 + T cells (C4) in an OTSCC patient (Fig. 5a). In this investigation, two CD4 + T cells were identified from CD3 + T cells and annotated as two distinct subtypes due to the significant variation in gene expression. The high expression of the regulatory marker genes FOXP3, CTLA4, and TIGIT identified the C4 cluster as regulatory T cells, as shown in Fig. 4b, whereas the C7 cluster expressed no canonical marker genes and was classified as CD4 + T cells. To elucidate the relationship between CD4 + T cell clusters, we visualized these cells on a pseudotime trajectory using diffusion maps (Fig. 5b).
Interestingly, CD4 + T cell clusters were divided into two sections by the first two diffusion components, and we observed two clusters that were tightly connected. Using Go enrichment, we found that C7 had a similar biological function as C4, including regulation of T cell differentiation and proliferation, as well as highly expressed regulatory-related genes, such as ANXA1, IL17F, CXCL13 (Additional file 7: Figure S3a,b). Figure 5C demonstrates that the distribution of T cell clonotypes with identical TCRs facilitated their interaction. Fisher’s exact test identified clusters of clonally expanded RT-specific CD4 + T cells, revealing an increase in regulatory CD4 + T cells in RT (Fisher’s exact test, FDR = 0.0244*(RT vs. PT) (Fig. 5g). Specifically, regulatory CD4 + T cells in RT comprised 62.8% of regulatory CD4 + T cells (C4 cluster), and this proportion increased to 88.7% when the background consisted of clonally expanded regulatory T cells (Fig. 5j).
C1 (transitional CD8 + T cells), C2 (Cytotoxic CD8 + T cells), C10 (regulatory CD8 + T cells), and C13 (Cytotoxic CD8 + T cells) were further investigated for CD8T cells. Signature genes expression was highly variable in these four CD8 + T cell clusters. The expression of regulatory-related genes such as ANXA1 and CXCR4, was more extensive and elevated in transitional CD8 + T cells and two cytotoxic genes. Fisher’s exact test revealed that the number of transitional CD8 + T cells (C1) and cytotoxic CD8 + T cells (C13) was significantly greater in RT (Fisher’s exact test, FDR = 0.000107 and 0.0204, respectively).
To better understand the relationship between these four clusters, we first performed pseudotime trajectory analyses, which revealed that these four subtypes were well-distributed in different states (Fig. 5d), and that regulatory CD8 + T cells and transitional CD8 + T cells were predominantly present in the beginning stage, whereas cytotoxic CD8 + T cells accumulated in the end stage (Fig. 5e). Similarly, clonally expanded T cells were predominantly enriched for cytotoxic CD8 + T cells (Fig. 5f). Then, we identified RT-specific clonally expanded CD8 + T cell clusters using Fisher’s exact test. Compared to normal tongue tissue and primary OTSCC, clonally expanded CD8 + T cells increased significantly in transitional CD8 + T cells (FDR = 0.0000504 (NT vs. RT) and 0.00000435 (PT vs. RT)) and cytotoxic CD8 + T cells (C1, C2 and C13 cluster) in recurrent tumor tissue (FDR = 3.72e-8 (NT vs. RT) and 5.16e-15 (PT vs. RT); FDR = 0.00918 (NT vs. RT) and 0.0143(PT vs. RT)) (Fig. 5g, Table 2). Specifically, transitional CD8 + T cells in RT comprised 48.1% of transitional cells (C1 cluster), and this proportion increased to 75.1% when the background consisted of clonally expanded transitional T cells (Fig. 5h). In addition, we visualized the relationships between C1, C2, C10, and C3 using diffusion maps (Fig. 5k). Both transitional CD8 + T cells and cytotoxic CD8 + T cells were derived from regulatory CD8 + T cells and correlated with each other. Figure 5L depicted the distribution of larger colonocytes at the branching point of regulatory CD8 + T cells. The expression of regulatory T cell-related marker genes FOXP3, IL2RA accumulated at the commencement of the branch and regulated T cell activation (Additional file 7: Figure S3c, d). The transitional-related marker gene GZMK was highly expressed in the center of the branch (Additional file 7: Figure S3e). The cytotoxic-related genes GZMB and PRF1 were highly expressed and tightly correlated with cytotoxic CD8 + T cells (C2 cluster) at extremity of the branch (Additional file 7: Figure S3f-h). These findings indicate that CD8 + cytotoxic T cells derived from TCR-activated regulatory T cells clonally expanded substantially in RT.
Immunosuppressive macrophage cells and T cells function inhibition
We investigated further the communication between macrophages and T cells in recurrent OTSCC. The CellPhoneDB analysis of communication between subsets of macrophage and T cells revealed that CD8 + T cells had a stronger interaction with macrophage than CD4 + T cells. Macrophage inhibition factor (MIF) was found to be widely distributed in all cell subpopulations and to play an essential role in the interaction between T cells and macrophage (Fig. 6a), as well as between tumor-related cells and immune cells. In addition, the immune inhibitory pathway (SPP1) and cytotoxic pathway (GALECTIN) were prominent in the interplay between macrophages and T cells in tumor tissue. The analysis also revealed that nearly all subpopulations of macrophages and T cells expressed the immune-inhibitory ligand MIF and CD74, whereas ligand encoding galectin 9 (LGALS9) was primarily expressed in C1Qs + TAM cells (Fig. 6b) and was significantly elevated in RT, indicating that macrophages play a crucial role in T cell suppression. The expression levels and interaction intensities of these immune-inhibitory ligands were lowest in normal tissue, highest in primary OTSCC, and lowest in recurrent OTSCC (Fig. 6c, d). These findings suggest that macrophage cells progressively formed an immune-inhibitory microenvironment as recurrent malignant tumors advanced.