HYPOTHESIS & THEORY

Pathol. Oncol. Res., 05 July 2023
https://doi.org/10.3389/pore.2023.1611210

Single-cell transcriptome analysis reveals the clinical implications of myeloid-derived suppressor cells in head and neck squamous cell carcinoma

www.frontiersin.orgWenru Jiang1, www.frontiersin.orgKangyao Hu1, www.frontiersin.orgXiaofei Liu1, www.frontiersin.orgJili Gao1* and www.frontiersin.orgLiping Zhu2*
  • 1Department of Implant and Prosthodontics, The First Affiliated Hospital of Harbin Medical University, Harbin, China
  • 2Department of Implant and Prosthodontics, Harbin First Hospital, Harbin, China

Head and neck squamous cell carcinoma (HNSC) is the most common malignant tumor that arises in the epithelium of the head and neck regions. Myeloid-derived suppressor cells (MDSCs) are one of the tumor-infiltrating immune cell populations, which play a powerful role in inhibiting anti-tumor immune response. Herein, we employed a single-cell RNA sequencing (scRNA-seq) dataset to dissect the heterogeneity of myeloid cells. We found that SPP1+ tumor-associated macrophages (TAMs) and MDSCs were the most abundant myeloid cells in the microenvironment. By cell cluster deconvolution from bulk RNA-seq datasets of larger patient groups, we observed that highly-infiltrated MDSC was a poor prognostic marker for patients’ overall survival (OS) probabilities. To better apply the MDSC OS prediction values, we identified a set of six MDSC-related genes (ALDOA, CD52, FTH1, RTN4, SLC2A3, and TNFAIP6) as the prognostic signature. In both training and test cohorts, MDSC-related prognostic signature showed a promising value for predicting patients’ prognosis outcomes. Further parsing the ligand-receptor pairs of intercellular communications by CellChat, we found that MDSCs could frequently interact with cytotoxic CD8+ T cells, SPP1+ TAMs, and endothelial cells. These interactions likely contributed to the establishment of an immunosuppressive microenvironment and the promotion of tumor angiogenesis. Our findings suggest that targeting MDSCs may serve as an alternative and promising target for the immunotherapy of HNSC.

Introduction

Head and neck squamous cell carcinoma (HNSC) is the sixth most common form of malignant tumor occurring in the epithelial tissues of the head and neck regions [13]. Despite the advancement in treatment approaches, such as surgical resection and multidisciplinary treatments involving radiotherapy and chemotherapy, the overall 5-year survival rate remains below 70%, particularly for patients in advanced stages of the disease [24]. Hence, comprehending the tumor microenvironment of HNSC, identifying therapeutic targets, and formulating novel treatment approaches are imperative.

Myeloid-derived suppressor cells (MDSCs) are a heterogeneous population of pathologically activated neutrophils and monocytes that accumulate in the tumor microenvironment and play a critical role in promoting tumor growth and immune evasion [5]. MDSCs have been identified in various types of cancer, including lung cancer, breast cancer, melanoma, and others [68]. MDSCs promote tumor growth by suppressing the immune response, particularly T cell activation, and function, which are critical for controlling cancer growth [9, 10]. Immunotherapies that target MDSCs include monoclonal antibodies and small molecule inhibitors. These therapies work by blocking the recruitment of MDSCs to the tumor microenvironment, inhibiting their function, and promoting T cell activation [1113]. A recent clinical study revealed that the administration of tadalafil, a phosphodiesterase-5 (PDE5) inhibitor, in HNSC patients resulted in a decrease in circulating MDSCs [14, 15]. Moreover, the treated patients exhibited reduced expression of iNOS and arginase in these cells, along with an increased presence of spontaneously generated tumor-specific T cells [14, 15]. Notably, in a co-culture system, HNSC cells had the potential to induce MDSCs differentiation from peripheral blood mononuclear cells and upregulate the expression of iNOS and arginase [16], which further indicated that immunotherapy strategies targeting MDSCs hold great promise.

In this study, we utilized a public single-cell RNA sequencing (scRNA-seq) dataset to investigate the heterogeneity of myeloid cells in HNSC. Our analysis revealed that SPP1+ tumor-associated macrophages (TAMs) and MDSCs were the most abundant myeloid cells in the tumor microenvironment. We also discovered that highly infiltrated MDSCs were associated with poorer overall survival rates in HNSC patients. To further explore the potential clinical application of MDSCs as a prognostic marker, we identified a set of six MDSC-related genes that could be used as a prognostic signature. Using this signature, we were able to predict patients’ prognosis outcomes with promising accuracy in both training and test cohorts. By examining intercellular communications, we found that MDSCs were able to suppress the activity of cytotoxic CD8+ T cells and recruit SPP1+ TAMs to shape an immunosuppressive microenvironment that promoted tumor angiogenesis. Overall, our findings suggest that targeting MDSCs may provide a promising therapeutic strategy for the immunotherapy of HNSC.

Materials and methods

Data collection

Single-cell dataset from Peng et al. was employed in our study, containing six tumor tissues of head and neck cancer (HNSC) patients (GSE172577) [17]. Bulk gene expression profiles of HNSC were downloaded from The Cancer Genome Atlas (TCGA) program. We also obtained the clinical information of corresponding samples of TCGA. In addition, two independent datasets GSE65858 [18] and GSE41613 [19] were employed as the test cohorts, which included 270 and 97 samples, respectively. The patient’s clinicopathological information is listed in Supplementary Table S1.

scRNA-seq data processing, batch correction, and clustering

We imported the unique molecular identifier (UMI) count data generated by 10x genomics into Seurat (V4.1.0) [20]. To remove the low-quality cells, we filtered (1) the cells with more than 20% mitochondrial counts; (2) cells expressing lower than 300 genes or more than 4,000 genes (Supplementary Figure S1). We also employed Scrublet with the default parameters to identify putative doublets [21]. The remaining 45,876 cells from six patients were normalized, scaled, and then used for batch correction. We took the Seurat-v3 batch correction strategy, anchors across patients were identified using the function FindIntegrationAnchors, and the data were finally integrated using the “IntegrateData” function. The assay “integrated” was used for downstream analysis. We next used the “FindVariableFeatures” function to choose the top 2000 highly variable genes based on the “vst” selection method. Principal component analysis (PCA) was performed and the top 30 PCA components were used for Uniform Manifold Approximation and Projection (UMAP) [17, 2224]. Subsequently, the cells were clustered on UMAP space using the Lovain algorithm on the k-nearest neighbors graph constructed using gene expression data as implemented in FindNeighbors and FindClusters. We further annotated major cell types according to the gene expression of well-known markers: T/NK cells (PTPRC, CD3D), mast cells (TPSAB1, TPSB2), B cells (CD79A, MS4A1), myeloid cells (C1QA, C1QB), fibroblasts (LUM, DCN), endothelial cells (PECAM1), keratinocytes (KRT15, KRT19), epithelial cells (KRT13, KRT14, EPCAM), and proliferating cells (TOP2A, MKI67).

Identification of myeloid subpopulations

We next subgrouped the myeloid cells. Based on the UMAP space, we applied the Lovain algorithm on the k-nearest neighbors graph using the function of FindNeighbors and FindClusters in the “Seurat” package. Markers were identified using the Wilcoxon Rank Sum test in each myeloid subgroup using the function FindAllMarkers.

Cell-type deconvolution

We employed the BayesPrism algorithm to deconvolute the infiltration levels of myeloid cells [25]. BayesPrism is a statistical framework for cell-type deconvolution, which is the process of inferring the proportion of cell types present in a heterogeneous mixture of cells based on gene expression data. The analysis was performed using the default parameters of the R package “BayesPrism” (https://github.com/Danko-Lab/BayesPrism).

CopyKAT analysis

CopyKAT [26] uses integrative Bayesian approaches to find genome-wide aneuploidy at 5MB resolution, and cells with large genome-wide aneuploidy were identified as tumor cells (Supplementary Figure S2). UMI count matrix was used as input, and others were default parameters (id.type = “S”, ngene.chr = 5, win.size = 25, and KS.cut = 0.1).

CytoTRACE analysis

We utilized the CytoTRACE [27] to compare differentiation states among HNSC tumor cells (https://cytotrace.stanford.edu/). CytoTRACE analyzes the number of uniquely expressed genes per cell, as well as other factors like distribution of mRNA content and the number of RNA copies per gene, to calculate a score assessing the differentiation and developmental potential of each cell (lowest differentiation and highest developmental potential at 1; highest differentiation and lowest developmental potential at 0). CytoTRACE analysis was conducted using default parameters. In addition, we evaluated the activities of cancer hallmark pathways from MSigdb (https://www.gsea-msigdb.org/gsea/) using the R package AUCell [28] with default parameters. Subsequently, Spearman’s correlation was calculated between the hallmark activity and CytoTRACE score.

Construction of the MDSC-related prognostic signature

We first divided the TCGA-HNSC cancer samples into three parts (two parts as the “training” set and one part as the “test” set) to apply 3-fold cross-validation. Then, we applied the univariate Cox regression model to screen the MDSC-related genes (log2FC > 0.25 and expression percentage >25%) that were associated with patients’ overall survival (OS) in the training set of the TCGA-HNSC cohort. Genes with p-value <0.05 were identified as the candidate prognosis-related genes (Supplementary Table S2). Afterward, we used a stepwise multivariate Cox regression model based on the Akaike information criterion (AIC) value to analyze the candidate MDSC-related genes and selected the ones that minimized AIC to achieve the best model fit [29]. We subsequently calculated the risk score [30] for each patient by the linear combination of expression values weighted by the coefficient from the multivariate Cox regression model, riskscore=in=6coefi*expi, where i represented the ith MDSC-related gene, exp denotes the expression levels of MDSC-related genes. We used the median value of patients’ risk scores to determine the high-risk and low-risk groups. Kaplan-Meier (KM) analysis with log-rank test was applied to compare the survival difference between patients’ risk groups using the R package “survival.”

CellChat analysis

We employed the CellChat computational tool to analyze communication among microenvironment cells [31]. CellChat uses a scoring system to identify the most likely interactions between cells based on the expression of genes encoding ligands and receptors.

In brief, we followed the official workflow and imported gene expression data using the function of createCellChat. We then applied the functions of identifyOverExpressedGenes, identifyOverExpressedInteractions, projectData to detect significant cell-cell interactions among the investigated cells. The analysis was conducted by the R package CellChat (https://github.com/sqjin/CellChat).

Results

The single-cell landscape of HNSC tumor microenvironment

To dissect the landscape of the tumor microenvironment of head and neck squamous cell carcinoma (HNSC), we obtained scRNA-seq data from 6 untreated HNSC primary patients (GSE172577) [17]. A total of 45,876 passed the initial quality control and were retained for downstream analysis (Supplementary Figures S1A,B). We merged all scRNA-seq data and performed gene expression normalization, scaling, dimension reduction, batch correction, and cell clustering to identify coarse cell types (Materials and methods). Nine major cell types were detected based on the gene expression of canonical cell markers, including epithelial cells, keratinocytes, fibroblasts, endothelial cells, myeloid cells, T/NK cells, mast cells, B cells, and proliferating cells (Figures 1A, B). The proportions of these major cell types exhibited significant variation among the different patients (Figure 1C). Notably, compared to epithelial cells, tumor-infiltrating immune cells (myeloid cells, T/NK cells, mast cells, and B cells) showed lower sample heterogeneity (Supplementary Figures S1C).

FIGURE 1
www.frontiersin.org

FIGURE 1. A global overview of TME in HNSC at single-cell resolution. (A) UMAP plot of major cell populations, dots colored by different cell types. (B) The expression level of cell type specific gene markers. (C) Bar plot showing the sample type fractions relative to the total cell count per cell type.

Dissecting the heterogeneity of myeloid cells

Myeloid cells play critical roles in regulating the immune response in the HNSC tumor microenvironment [32, 33]. Overall, 3,026 myeloid cells were further reclassified into nine populations (Figure 2A). Known representative genes were used to recognize cell identities (Figure 2B), including three subtypes of dendritic cells (DCs): cDC1 (BATE3, XCR1), cDC2 (CLEC10A, FCER1A), DC3 (LAMP3, CCR7); two types of tumor-associated macrophages (TAMs): SPP1+ TAM (SPP1, MRC1), CXCL9+ TAM (CXCL9); monocytes (FCN1); myeloid-derived suppressor cells (MDSCs) (S100A8, S100A9, IL1B); Langerhans cells (CD1A, CD207); proliferating cells (MKI67, TOP2A). We next investigated the proportions of myeloid subpopulations among HNSC patients (Figure 2C). The result showed that SPP1+ TAM and MDSC had higher fractions compared to other myeloid cells. Notably, both SPP1+ TAM and MDSC played a critical role in shaping the immunosuppressive tumor microenvironment [3437].

FIGURE 2
www.frontiersin.org

FIGURE 2. Subpopulations of myeloid cells in HNSC microenvironment. (A) UMAP plot of myeloid subpopulations, dots colored by different cell types. (B) Dot plot showing the expression level of cell-type-specific gene markers. (C) Box plot showing the fractions of cells in each sample.

Highly-infiltrated MDSCs prompt poor prognostic risks

To understand the role of myeloid cells in HNSC patients, we performed the deconvolution analysis using the BayesPrism algorithm [25]. We evaluated the infiltration levels of myeloid cells in three independent bulk RNA-seq datasets (TCGA, GSE65858 [18], and GSE41613 [19]) (Figure 3A). By constructing the univariate Cox regression model, we found only the infiltration levels of MDSC showed a significant association with the patient’s overall survival (OS) probability in all three datasets (Figure 3B). Patients with highly-infiltrated MDSCs usually suffer poor prognostic outcomes, suggesting MDSC infiltration was a prognostic risk factor (HR > 1 and p-value <0.05). In addition, we constructed the multivariate Cox regression model to explore the correlation between MDSC infiltration and patient OS probability adjusting the effects of tumor stage, age, and gender. The result showed that the MDSC infiltration was an independent factor for predicting the prognostic risk (Figure 3C).

FIGURE 3
www.frontiersin.org

FIGURE 3. Tumor-infiltrated MDSC was associated with patients’ prognostic outcomes. (A) Bar plot showing the infiltration levels of myeloid subpopulations evaluated by the BayesPrism algorithm. (B) Forest plot showing the result of univariate Cox-regression analysis for correlation between the myeloid cell infiltration levels and the overall survival. (C) Forest plot showing the result of multivariate Cox-regression analysis of the MDSC infiltration levels and the overall survival adjusting the effects of tumor stage, age, and gender. Three independent datasets were used: TCGA cohort (n = 491) (top), GSE65858 (n = 270) (middle), and GSE41613 (n = 97) (bottom), respectively. (D) Boxplot depicting the comparison of CytoTRACE scores between tumor cells in HIM [high-infiltrated MDSC tumor microenvironment (TME)] group and those in LIM [low-infiltrated MDSC TME] group. p-value was calculated by the Kruskal-Wallis test. (E) Volcano plot showing the correlation between CytoTRACE score and the activity of cancer hallmark pathway. (F) Violin plot showing the expression levels of EMT-related genes in tumor cells in high/low-infiltrated MDSC TME. *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001, as calculated by Mann Whitney U test.

We next explored the impact of MDSC infiltration on tumor cells in the tumor microenvironment (TME). Based on the median levels of MDSC infiltration, we divided the high-infiltrated-MDSC TME (HIM) (patients SYSMH2, SYSMH3, SYSMH5) and low-infiltrated-MDSC TME (LIM) (patients SYSMH1, SYSMH4, and SYSMH6). We then used CopyKAT [26] to infer copy number alterations at 5 Mb resolution by averaging large chromosomal regions (1 Mbp) and identified the HNSC tumor cells in each sample (Supplementary Figure S2). By applying the CytoTRACE algorithm to evaluate tumor cell differentiation states [27], we found that tumor cells in HIM TME had higher developmental potential than those in LIM TME (Figure 3D). Further correlation analysis indicated that apical junctions, mTORC1 signaling, and epithelial-mesenchymal transition (EMT) were the most associated processes with tumor cell differentiation (Figure 3E) (See the ‘Methods’ section). Specifically, we also found some EMT-related genes were more highly expressed in the HIM tumor cells than in LIM tumor cells (Figure 3F). These results suggested that MDSC infiltration might play an important role in tumor cell development.

An MDSC-related gene signature shows robust prognostic predictive values

To better apply MDSC in survival prediction, we explored the prognosis values of MDSC-related genes (log2FC > 0.25 and expression percentage >25%). We divided the TCGA-HNSC cancer samples into three parts, where the TCGA training set had twice as many patients as the TCGA test set, and constructed the univariate Cox regression model based on the expression levels of MDSC-related genes in the TCGA-training set (Supplementary Table S2). Subsequently, multivariate Cox regression and stepwise regression models were employed for identifying the prognostic signature. A set of six MDSC-related genes (ALDOA, CD52, FTH1, RTN4, SLC2A3, and TNFAIP6) were finally trained as the prognostic signature (Figure 4A) (Materials and methods). Further analyzing the expression levels of the signature genes, we found that TNFAIP6 (TNF alpha induced protein 6) was an MDSC-specific marker among myeloid subpopulations (Figure 4B; Supplementary Table S3). In addition, we employed a risk-scoring model based on the MDSC prognostic signature. The risk score of each patient could be expressed as Risk score = (0.0026681*ALDOA) + (−0.0126845*CD52) + (0.0015119*FTH1) + (0.0125243*RTN4) + (0.0154432*SLC2A3) + (0.0052426*TNFAIP6). We next subgrouped patients into high- and low-risk groups based on the median value of risk scores. Compared to the low-risk group, patients in the high-risk group showed worse prognosis outcomes in the TCGA training set, TCGA test set, and two external test datasets (GSE65858 and GSE41613) (Figure 4C). Moreover, we also evaluated the capacity of the MDSC risk group as the independent prognostic predictor. By adjusting the effects of tumor stage, age, and gender, the risk group is still a robust poor prognostic signature (Figure 4D). These results showed that the MDSC-related genes exhibited outstanding performance in HNSC prognostic prediction.

FIGURE 4
www.frontiersin.org

FIGURE 4. The MDSC-related prognostic signature of HNSC. (A) Forest plot showing the result of multivariate Cox-regression analysis for correlation between the MDSC-related genes and the overall survival. (B) Violin plot showing the expression levels of the MDSC-related prognostic signature across myeloid subpopulations. (C) Kaplan-Meier curve of HNSC samples stratified by the risk groups with log-rank test p-value provided. (D) Forest plot showing the result of multivariate Cox-regression analysis for correlation between the risk group and the overall survival by adjusting the effects of gender, age, and tumor stages.

Intercellular communication associated with MDSCs in HNSC

We next characterized the role of MDSC in the tumor microenvironment by dissecting intercellular communications. We first subgrouped the T/NK cells (Figure 5A) and fibroblasts (Supplementary Figure S4). For T/NK cells, we identified five subpopulations and annotated them by the known gene markers: naïve/memory T-cells (Tn/Tm) (IL7R, SELL), cytotoxic CD8+ T-cells (CD8A, CD8B, GZMA, GZMB), regulatory T-cells (Treg) (FOXP3), exhausted T-cells (Tex) (PDCD1, CTLA4), and natural killer cells (NK) (KLRF1, KLRD1, TRDC) (Figure 5B). For fibroblasts, we recognized five subgroups, including APOE+ cancer-associated fibroblasts (CAFs), APOD+ CAFs, myofibroblastic CAFs (myoCAFs), inflammatory CAFs (iCAFs), and proliferating cells (Supplementary Figure S4). Subsequently, we employed CellChat to explore the crosstalk between MDSCs and other microenvironment cells (Figure 5C). The result showed that MDSC mainly communicated with cytotoxic CD8+ T-cells, SPP1+ TAMs, and endothelial cells. Further parsing of the ligand-receptor (LR) pairs between MDSC and its three partners, we found that several non-classical MHC class I molecules, including HLA-E and HLA-F, mediated the communications from MDSC to cytotoxic CD8+ T-cells (Figure 5D), which was important for shaping immunosuppressive microenvironment [38, 39]. In addition, MDSC could interplay with SPP1+ TAMs via the LR pair CCL3/CCL3L3-CCR1 (Figure 5D), which might facilitate the recruitment of SPP1+ TAMs in tumor tissue [40, 41]. We also observed that MDSC-secreted VEGFA and multiple chemokines CCL2, CXCL1, CXCL2, CXCL3, and CXCL8 could act on endothelial cells via VEGFA-KDR/FLT1 and CCL2/CXCLs-ACKR1 interactions, which were crucial for tumor angiogenesis [42, 43]. We also evaluated the clinical implications of MDSC-related LR pairs by analyzing their associations with the probabilities of patients’ OS and progression-free survival (PFS) (Supplementary Figure S5). We found some immunosuppressive (CD86-CTLA4, HLA-E-CD8A/CD94_NKG2C/CD94_NKG2A, HLA-F-CD8A) and angiogenesis-related (CXCL1-ACKR1) LR pairs showed a significant association with patients’ OS and PFS. These results highlight the potential of MDSC in shaping the immunosuppressive microenvironment and promoting tumor angiogenesis, which is also associated with tumor development and progression.

FIGURE 5
www.frontiersin.org

FIGURE 5. Intercellular communications between MDSCs and other microenvironment cells. (A) UMAP plot of T/NK subpopulations, dots colored by different cell types. (B) Dot plot showing the expression level of cell-type-specific gene markers. (C) Circle plots displaying the inferred network of MDSC-related cell communications. Edge width is proportional to the inferred communication counts. (D) Dot plot showing the ligand-receptor pairs from MDSC to endothelial cell, cytotoxic CD8+ T cell, and SPP1+ TAM. Dots are sized by the inferred communication probabilities.

Discussion

Head and neck squamous cell carcinoma (HNSC) is one of the most common malignancies in the head and neck region. In this study, scRNA-seq data was used to dissect the landscape of the tumor microenvironment of HNSC. We found that there were nine major cell types, including epithelial cells, keratinocytes, fibroblasts, endothelial cells, myeloid cells, T/NK cells, mast cells, B cells, and proliferating cells. Among these, myeloid cells play a critical role in regulating the immune response in the HNSC tumor microenvironment. We further analyzed the subpopulations of myeloid cells, including dendritic cells (DCs), tumor-associated macrophages (TAMs), monocytes, myeloid-derived suppressor cells (MDSCs), Langerhans cells, and proliferating cells. We investigated the proportions of myeloid subpopulations among HNSC patients and found that SPP1+ TAM and MDSC had higher fractions compared to other myeloid cells. We also found that highly-infiltrated MDSCs prompted poor prognostic risks, suggesting that MDSC infiltration was a prognostic risk factor. Recent studies also demonstrated that the infiltration levels of MDSC were increased in the HNSC tumor tissues, and their presence was positively associated with advanced T stage, higher pathological grade, lymph node metastasis, and poor prognosis [16]. Furthermore, in a co-culture system, tumor-related MDSCs were found to promote the progression of HNSC by enhancing cell proliferation, migration, epithelial-mesenchymal transition (EMT), and vasculogenic mimicry (VM). These findings suggest a reciprocal interaction between MDSCs and tumor cells, facilitating the malignant progression of HNSC and enhancing the immunosuppressive properties of MDSCs.

To better apply MDSCs in survival prediction, we explored the prognosis values of MDSC-related genes. We identified a set of six MDSC-related genes that were trained as the prognostic signature, including ALDOA, CD52, FTH1, RTN4, SLC2A3, and TNFAIP6. We also employed a risk-scoring model based on the MDSC prognostic signature, and patients were subgrouped into high- and low-risk groups based on the median value of risk scores. The high-risk group showed worse prognosis outcomes in the TCGA training set, TCGA test set and two external test datasets, suggesting that the MDSC risk score could be a useful independent prognostic predictor.

Further analysis of MDSC-related cell communications, we found MDSC showed the potential to recruit SPP1+ TAM, inhibit the activity of cytotoxic T cells, and promote endothelial outgrowth, which was associated with immunosuppressive TME. In addition to MDSC, SPP1+ TAM also played an important role in shaping the immunosuppressive TME [44]. Notably, the prognostic signature also showed an association with hypoxic immunosuppressive TME. In addition to TNFAIP6, another five exhibited a general expression across myeloid subtypes (Figure 4B). These genes may also be regulated by the immunosuppressive microenvironment [4547], and their association with MDSCs will be investigated in our future study. For instance, ALDOA and SLC2A3 were found to be linked to glycolysis, indicating their potential role in energy metabolism [4850]; RTN4 was identified as a possible contributor to tumor angiogenesis, suggesting its involvement in the formation of new blood vessels to support tumor growth [51]; TNFAIP6 emerged as a significant regulator of extracellular matrix organization, implying its influence on the structural integrity and composition of the tumor microenvironment [5258].

Overall, this study provides valuable insights into the landscape of the tumor microenvironment of HNSC and highlights the critical role of myeloid cells in regulating the immune response. The findings also suggest that MDSC infiltration is a prognostic risk factor in HNSC patients and that the MDSC prognostic signature could be a useful tool for predicting patient outcomes. These results may have important implications for the development of novel immunotherapeutic strategies for HNSC.

Data availability statement

All the data used in this study could be accessed freely in GEO and TCGA.

Author contributions

WJ, KH, and XL designed the study and performed the analysis. JG and LZ supervised the research. All authors contributed to the article and approved the submitted version.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.por-journal.com/articles/10.3389/pore.2023.1611210/full#supplementary-material

References

1. Sung, H, Ferlay, J, Siegel, RL, Laversanne, M, Soerjomataram, I, Jemal, A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin (2021) 71(3):209–49. doi:10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Johnson, DE, Burtness, B, Leemans, CR, Lui, VWY, Bauman, JE, and Grandis, JR. Head and neck squamous cell carcinoma. Nat Rev Dis Primers (2020) 6(1):92. doi:10.1038/s41572-020-00224-3

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Pulte, D, and Brenner, H. Changes in survival in head and neck cancers in the late 20th and early 21st century: A period analysis. Oncologist (2010) 15(9):994–1001. doi:10.1634/theoncologist.2009-0289

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Chen, W, Zheng, R, Baade, PD, Zhang, S, Zeng, H, Bray, F, et al. Cancer statistics in China, 2015. CA: A Cancer J Clinicians (2016) 66(2):115–32. doi:10.3322/caac.21338

CrossRef Full Text | Google Scholar

5. Veglia, F, Sanseviero, E, and Gabrilovich, DI. Myeloid-derived suppressor cells in the era of increasing myeloid cell diversity. Nat Rev Immunol (2021) 21(8):485–98. doi:10.1038/s41577-020-00490-y

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Al-Khami, AA, Rodriguez, PC, and Ochoa, AC. Metabolic reprogramming of myeloid-derived suppressor cells (MDSC) in cancer. Oncoimmunology (2016) 5(8):e1200771. doi:10.1080/2162402X.2016.1200771

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Ma, P, Beatty, PL, McKolanis, J, Brand, R, Schoen, RE, and Finn, OJ. Circulating myeloid derived suppressor cells (MDSC) that accumulate in premalignancy share phenotypic and functional characteristics with MDSC in cancer. Front Immunol (2019) 10:1401. doi:10.3389/fimmu.2019.01401

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Albeituni, SH, Ding, C, and Yan, J. Hampering immune suppressors: Therapeutic targeting of myeloid-derived suppressor cells in cancer. Cancer J (Sudbury, Mass) (2013) 19(6):490–501. doi:10.1097/PPO.0000000000000006

CrossRef Full Text | Google Scholar

9. Noman, MZ, Desantis, G, Janji, B, Hasmim, M, Karray, S, Dessen, P, et al. PD-L1 is a novel direct target of HIF-1α, and its blockade under hypoxia enhanced MDSC-mediated T cell activation. J Exp Med (2014) 211(5):781–90. doi:10.1084/jem.20131916

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Schouppe, E, Mommer, C, Movahedi, K, Laoui, D, Morias, Y, Gysemans, C, et al. Tumor-induced myeloid-derived suppressor cell subsets exert either inhibitory or stimulatory effects on distinct CD 8+ T-cell activation events. Eur J Immunol (2013) 43(11):2930–42. doi:10.1002/eji.201343349

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Fultang, L, Panetti, S, Ng, M, Collins, P, Graef, S, Rizkalla, N, et al. MDSC targeting with Gemtuzumab ozogamicin restores T cell immunity and immunotherapy against cancers. EBioMedicine (2019) 47:235–46. doi:10.1016/j.ebiom.2019.08.025

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Li, T, Liu, T, Zhu, W, Xie, S, Zhao, Z, Feng, B, et al. Targeting MDSC for immune-checkpoint blockade in cancer immunotherapy: Current progress and new prospects. Clin Med Insights: Oncol (2021) 15:11795549211035540. doi:10.1177/11795549211035540

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Ghansah, T. A novel strategy for modulation of MDSC to enhance cancer immunotherapy. Oncoimmunology (2012) 1(6):984–5. doi:10.4161/onci.20201

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Califano, JA, Khan, Z, Noonan, KA, Rudraraju, L, Zhang, Z, Wang, H, et al. Tadalafil augments tumor specific immunity in patients with head and neck squamous cell carcinoma. Clin Cancer Res (2015) 21(1):30–8. doi:10.1158/1078-0432.CCR-14-1716

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Weed, DT, Vella, JL, Reis, IM, De la Fuente, AC, Gomez, C, Sargi, Z, et al. Tadalafil reduces myeloid-derived suppressor cells and regulatory T cells and promotes tumor immunity in patients with head and neck squamous cell carcinoma. Clin Cancer Res (2015) 21(1):39–48. doi:10.1158/1078-0432.CCR-14-1711

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Pang, X, Fan, HY, Tang, YL, Wang, SS, Cao, MX, Wang, HF, et al. Myeloid derived suppressor cells contribute to the malignant progression of oral squamous cell carcinoma. PLoS One (2020) 15(2):e0229089. doi:10.1371/journal.pone.0229089

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Peng, Y, Xiao, L, Rong, H, Ou, Z, Cai, T, Liu, N, et al. Single-cell profiling of tumor-infiltrating TCF1/TCF7+ T cells reveals a T lymphocyte subset associated with tertiary lymphoid structures/organs and a superior prognosis in oral cancer. Oral Oncol (2021) 119:105348. doi:10.1016/j.oraloncology.2021.105348

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Wichmann, G, Rosolowski, M, Krohn, K, Kreuz, M, Boehm, A, Reiche, A, et al. The role of HPV RNA transcription, immune response-related gene expression and disruptive TP53 mutations in diagnostic and prognostic profiling of head and neck cancer. Int J Cancer (2015) 137(12):2846–57. doi:10.1002/ijc.29649

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Lohavanichbutr, P, Méndez, E, Holsinger, FC, Rue, TC, Zhang, Y, Houck, J, et al. A 13-gene signature prognostic of HPV-negative OSCC: Discovery and external validation. Clin Cancer Res (2013) 19(5):1197–203. doi:10.1158/1078-0432.CCR-12-2647

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Butler, A, Hoffman, P, Smibert, P, Papalexi, E, and Satija, R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol (2018) 36(5):411–20. doi:10.1038/nbt.4096

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Wolock, SL, Lopez, R, and Klein, AM. Scrublet: Computational identification of cell doublets in single-cell transcriptomic data. Cell Syst (2019) 8(4):281–91.e9. doi:10.1016/j.cels.2018.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Kamath, T, Abdulraouf, A, Burris, SJ, Langlieb, J, Gazestani, V, Nadaf, NM, et al. Single-cell genomic profiling of human dopamine neurons identifies a population that selectively degenerates in Parkinson's disease. Nat Neurosci (2022) 25(5):588–95. doi:10.1038/s41593-022-01061-1

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Huang, L, Wang, X, Pei, S, Li, X, Dong, L, Bian, X, et al. Single-cell profiling reveals sustained immune infiltration, surveillance, and tumor heterogeneity in infiltrative BCC. J Invest Dermatol (2023). doi:10.1016/j.jid.2023.04.020

CrossRef Full Text | Google Scholar

24. Li, X, Zhao, S, Bian, X, Zhang, L, Lu, L, Pei, S, et al. Signatures of EMT, immunosuppression, and inflammation in primary and recurrent human cutaneous squamous cell carcinoma at single-cell resolution. Theranostics (2022) 12(17):7532–49. doi:10.7150/thno.77528

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Chu, T, Wang, Z, Pe’er, D, and Danko, CG. Cell type and gene expression deconvolution with BayesPrism enables Bayesian integrative analysis across bulk and single-cell RNA sequencing in oncology. Nat Cancer (2022) 3(4):505–17. doi:10.1038/s43018-022-00356-3

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Gao, R, Bai, S, Henderson, YC, Lin, Y, Schalck, A, Yan, Y, et al. Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes. Nat Biotechnol (2021) 39(5):599–608. doi:10.1038/s41587-020-00795-2

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Gulati, GS, Sikandar, SS, Wesche, DJ, Manjunath, A, Bharadwaj, A, Berger, MJ, et al. Single-cell transcriptional diversity is a hallmark of developmental potential. Science (2020) 367(6476):405–11. doi:10.1126/science.aax0249

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Aibar, S, González-Blas, CB, Moerman, T, Huynh-Thu, VA, Imrichova, H, Hulselmans, G, et al. Scenic: Single-cell regulatory network inference and clustering. Nat Methods (2017) 14(11):1083–6. doi:10.1038/nmeth.4463

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Wang, J, Yin, X, Zhang, YQ, and Ji, X. Identification and validation of a novel immune-related four-lncRNA signature for lung adenocarcinoma. Front Genet (2021) 12:639254. doi:10.3389/fgene.2021.639254

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Pan, X, Zhang, C, Wang, J, Wang, P, Gao, Y, Shang, S, et al. Epigenome signature as an immunophenotype indicator prompts durable clinical immunotherapy benefits in lung adenocarcinoma. Brief Bioinform (2022) 23(1):bbab481. doi:10.1093/bib/bbab481

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Jin, S, Guerrero-Juarez, CF, Zhang, L, Chang, I, Ramos, R, Kuan, C-H, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun (2021) 12(1):1088. doi:10.1038/s41467-021-21246-9

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Chu, M, Su, YX, Wang, L, Zhang, TH, Liang, YJ, Liang, LZ, et al. Myeloid-derived suppressor cells contribute to oral cancer progression in 4NQO-treated mice. Oral Dis (2012) 18(1):67–73. doi:10.1111/j.1601-0825.2011.01846.x

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Jiang, Y, Wang, C, Wang, Y, Zhang, W, Liu, L, and Cheng, J. Prognostic role of CD11b(+) myeloid-derived suppressor cells in oral squamous cell carcinoma. Arch Med Sci (2023) 19(1):171–9. doi:10.5114/aoms/116683

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Qi, J, Sun, H, Zhang, Y, Wang, Z, Xun, Z, Li, Z, et al. Single-cell and spatial analysis reveal interaction of FAP+ fibroblasts and SPP1+ macrophages in colorectal cancer. Nat Commun (2022) 13(1):1742. doi:10.1038/s41467-022-29366-6

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Georgoudaki, AM, Prokopec, KE, Boura, VF, Hellqvist, E, Sohn, S, Östling, J, et al. Reprogramming tumor-associated macrophages by antibody targeting inhibits cancer progression and metastasis. Cell Rep (2016) 15(9):2000–11. doi:10.1016/j.celrep.2016.04.084

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Groth, C, Hu, X, Weber, R, Fleming, V, Altevogt, P, Utikal, J, et al. Immunosuppression mediated by myeloid-derived suppressor cells (MDSCs) during tumour progression. Br J Cancer (2019) 120(1):16–25. doi:10.1038/s41416-018-0333-1

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Li, K, Shi, H, Zhang, B, Ou, X, Ma, Q, Chen, Y, et al. Myeloid-derived suppressor cells as immunosuppressive regulators and therapeutic targets in cancer. Signal Transduction Targeted Ther (2021) 6(1):362. doi:10.1038/s41392-021-00670-9

CrossRef Full Text | Google Scholar

38. Kochan, G, Escors, D, Breckpot, K, and Guerrero-Setas, D. Role of non-classical MHC class I molecules in cancer immunosuppression. Oncoimmunology (2013) 2(11):e26491. doi:10.4161/onci.26491

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Liu, X, Song, J, Zhang, H, Liu, X, Zuo, F, Zhao, Y, et al. Immune checkpoint HLA-E:CD94-NKG2A mediates evasion of circulating tumor cells from NK cell surveillance. Cancer Cell (2023) 41(2):272–87.e9. doi:10.1016/j.ccell.2023.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Ma, R-Y, Black, A, and Qian, B-Z. Macrophage diversity in cancer revisited in the era of single-cell omics. Trends Immunol (2022) 43(7):546–63. doi:10.1016/j.it.2022.04.008

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Nagarsheth, N, Wicha, MS, and Zou, W. Chemokines in the cancer microenvironment and their relevance in cancer immunotherapy. Nat Rev Immunol (2017) 17(9):559–72. doi:10.1038/nri.2017.49

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Potente, M, and Mäkinen, T. Vascular heterogeneity and specialization in development and disease. Nat Rev Mol Cel Biol (2017) 18(8):477–94. doi:10.1038/nrm.2017.36

CrossRef Full Text | Google Scholar

43. Salazar, N, and Zabel, BA. Support of tumor endothelial cells by chemokine receptors. Front Immunol (2019) 10:147. doi:10.3389/fimmu.2019.00147

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Zhang, L, Li, Z, Skrzypczynska, KM, Fang, Q, Zhang, W, O’Brien, SA, et al. Single-cell analyses inform mechanisms of myeloid-targeted therapies in colon cancer. Cell (2020) 181(2):442–59.e29. doi:10.1016/j.cell.2020.03.048

PubMed Abstract | CrossRef Full Text | Google Scholar

45. He, Z, Chen, Q, He, W, Cao, J, Yao, S, Huang, Q, et al. Hepatocellular carcinoma subtypes based on metabolic pathways reveals potential therapeutic targets. Front Oncol (2023) 13:1086604. doi:10.3389/fonc.2023.1086604

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Chiu, C-F. Abstract 378: Ferritin regulates oxaliplatin sensitivity and immunosuppression in pancreatic ductal adenocarcinoma. Cancer Res (2023) 83:378. 7_Supplement. doi:10.1158/1538-7445.am2023-378

CrossRef Full Text | Google Scholar

47. Li, L, Yang, L, Chen, X, Chen, X, Diao, L, Zeng, Y, et al. TNFAIP6 defines the MSC subpopulation with enhanced immune suppression activities. Stem Cel Res Ther (2022) 13(1):479. doi:10.1186/s13287-022-03176-5

CrossRef Full Text | Google Scholar

48. Kuang, Q, Liang, Y, Zhuo, Y, Cai, Z, Jiang, F, Xie, J, et al. The ALDOA metabolism pathway as a potential target for regulation of prostate cancer proliferation. Onco Targets Ther (2021) 14:3353–66. doi:10.2147/OTT.S290284

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Zhang, F, Lin, J-D, Zuo, X-Y, Zhuang, Y-X, Hong, C-Q, Zhang, G-J, et al. Elevated transcriptional levels of aldolase A (ALDOA) associates with cell cycle-related genes in patients with NSCLC and several solid tumors. BioData Mining (2017) 10(1):6. doi:10.1186/s13040-016-0122-4

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Yao, X, He, Z, Qin, C, Deng, X, Bai, L, Li, G, et al. SLC2A3 promotes macrophage infiltration by glycolysis reprogramming in gastric cancer. Cancer Cel Int (2020) 20:503. doi:10.1186/s12935-020-01599-9

CrossRef Full Text | Google Scholar

51. Xiao, P, Gu, J, Xu, W, Niu, X, Zhang, J, Li, J, et al. RTN4/Nogo-A-S1PR2 negatively regulates angiogenesis and secondary neural repair through enhancing vascular autophagy in the thalamus after cerebral cortical infarction. Autophagy (2022) 18(11):2711–30. doi:10.1080/15548627.2022.2047344

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Sai, B, Dai, Y, Fan, S, Wang, F, Wang, L, Li, Z, et al. Cancer-educated mesenchymal stem cells promote the survival of cancer cells at primary and distant metastatic sites via the expansion of bone marrow-derived-PMN-MDSCs. Cell Death Dis (2019) 10(12):941. doi:10.1038/s41419-019-2149-1

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Desbois, M, Udyavar, AR, Ryner, L, Kozlowski, C, Guan, Y, Dürrbaum, M, et al. Integrated digital pathology and transcriptome analysis identifies molecular mediators of T-cell exclusion in ovarian cancer. Nat Commun (2020) 11(1):5583. doi:10.1038/s41467-020-19408-2

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Choi, H, Lee, RH, Bazhanov, N, Oh, JY, and Prockop, DJ. Anti-inflammatory protein TSG-6 secreted by activated MSCs attenuates zymosan-induced mouse peritonitis by decreasing TLR2/NF-κB signaling in resident macrophages. Blood (2011) 118(2):330–8. doi:10.1182/blood-2010-12-327353

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Mittal, M, Tiruppathi, C, Nepal, S, Zhao, YY, Grzych, D, Soni, D, et al. TNFα-stimulated gene-6 (TSG6) activates macrophage phenotype transition to prevent inflammatory lung injury. Proc Natl Acad Sci U S A (2016) 113(50):E8151–e8158. doi:10.1073/pnas.1614935113

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Kuznetsova, SA, Mahoney, DJ, Martin-Manso, G, Ali, T, Nentwich, HA, Sipes, JM, et al. TSG-6 binds via its CUB_C domain to the cell-binding domain of fibronectin and increases fibronectin matrix assembly. Matrix Biol (2008) 27(3):201–10. doi:10.1016/j.matbio.2007.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Lawrance, W, Banerji, S, Day, AJ, Bhattacharjee, S, and Jackson, DG. Binding of hyaluronan to the native lymphatic vessel endothelial receptor LYVE-1 is critically dependent on receptor clustering and hyaluronan organization. J Biol Chem (2016) 291(15):8014–30. doi:10.1074/jbc.M115.708305

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Azzaoui, I, Uhel, F, Rossille, D, Pangault, C, Dulong, J, Le Priol, J, et al. T-cell defect in diffuse large B-cell lymphomas involves expansion of myeloid-derived suppressor cells. Blood (2016) 128(8):1081–92. doi:10.1182/blood-2015-08-662783

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: head and neck squamous cell carcinoma, myeloid-derived suppressor cell, single cell RNA sequencing (scRNA-seq), prognosis, immunotherapeutic

Citation: Jiang W, Hu K, Liu X, Gao J and Zhu L (2023) Single-cell transcriptome analysis reveals the clinical implications of myeloid-derived suppressor cells in head and neck squamous cell carcinoma. Pathol. Oncol. Res. 29:1611210. doi: 10.3389/pore.2023.1611210

Received: 27 March 2023; Accepted: 28 June 2023;
Published: 05 July 2023.

Edited by:

Andrea Ladányi, National Institute of Oncology (NIO), Hungary

Copyright © 2023 Jiang, Hu, Liu, Gao and Zhu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Jili Gao, gjl8977@163.com; Liping Zhu, 1920636476@qq.com

These authors have contributed equally to this work