Single-cell RNA sequencing reveals tumor immune microenvironment in human hypopharygeal squamous cell carcinoma and lymphatic metastasis

Background Human hypopharygeal squamous cell carcinoma (HSCC) is a common head and neck cancer with a poor prognosis in advanced stages. The occurrence and development of tumor is the result of mutual influence and co-evolution between tumor cells and tumor microenvironment (TME). Tumor immune microenvironment (TIME) refers to the immune microenvironment surrounding tumor cells. Studying TIME in HSCC could provide new targets and therapeutic strategies for HSCC. Methods We performed single-cell RNA sequencing (scRNA-seq) and analysis of hypopharyngeal carcinoma, paracancerous, and lymphoid tissues from five HSCC patients. Subdivide of B cells, T cells, macrophages cells, and monocytes and their distribution in three kinds of tissues as well as marker genes were analyzed. Different genes of IGHG1 plasma cells and SPP1+ macrophages between HSCC tissues, adjacent normal tissues and lymphatic tissues were analyzed. Additionally, we studied proliferating lymphocytes, T cells exhaustion, and T cell receptor (TCR) repertoire in three kinds of tissues. Results Transcriptome profiles of 132,869 single cells were obtained and grouped into seven cell clusters, including epithelial cells, lymphocytes, mononuclear phagocytics system (MPs), fibroblasts, endothelial cells (ECs), plasmacytoid dendritic cells (pDCs), and mast cells. Tumor metastasis occurred in three lymphoid tissues. Four distinct populations were identified from lymphocytes, including B cells, plasma cells, T cells and proliferating lymphocytes. We found IGHA1 and IGHG1 specific plasma cells significantly overexpressed in HSCC tissues compared with normal hypopharygeal tissues and lymphatic tissues. Five distinct populations from MPs were identified, including macrophages, monocytes, mature dendritic cells (DCs), Type 1 conventional dendritic cells (cDC1) and Type 2 conventional dendritic cells (cDC2). SPP1+ macrophages were significantly overexpressed in HSCC tissues and lymphatic tissues compared with normal hypopharygeal tissues, which are thought to be M2-type macrophages. Exhaustion of CD8+ Teff cells occurred in HSCC tissues. At last, we verified that IgA and IgG1 protein expression levels were significantly up-regulated in HSCC tissues compared to adjacent normal tissues. Conclusion Overall, this study revealed TIME in HSCC and lymphatic metastasis, and provided potential therapeutic targets for HSCC.


Introduction
Hypopharyngeal cancer accounts for 3% to 5% of head and neck squamous cell carcinomas (1). The vast majority of hypopharyngeal cancers are hypopharyngeal squamous cell carcinomas (HSCC) (2). Patients with HSCC are predominantly male and usually have a history of smoking or heavy alcohol consumption (3). Hypopharyngeal has the worst prognosis of all head and neck cancers, with a reported 5-year overall survival rate of 30% to 35% (4). Hypopharyngeal cancer is usually detected late, with 70% to 85% of cases diagnosed at stage III or IV, approximately 60% to 80% of patients have ipsilateral cervical lymph node metastases (5), and 40% have contralateral cervical occult lymph node metastases (6). Relapse is fairly common, with nearly 50% of patients relapse within one year of diagnosis, often with distant metastases diagnosed (7). For patients with locally advanced hypopharyngeal carcinoma, surgery combined with chemoradiotherapy can significantly improve local control rates and systemic efficacy (8). Although the application of emerging therapeutic approaches such as immunotherapy and targeted therapy have greatly improved the prognosis of cancer patients, the survival of HSCC patients has not been significantly prolonged (9).
Single-cell transcriptome sequencing (scRNA-seq) technology is a technology which isolates single cells from cell populations in tissue or body fluid samples, obtains relevant data through unbiased, high-throughput and high-resolution transcriptome sequencing, and finally performs informative analysis (10).
ScRNA-seq can not only reveal the unique changes of each cell, but also discover completely new cell types (11). Single-cell genome sequencing can be applied to detect genomic stability and genomic variation, thus providing new insights into human understanding of the physiological and pathological functions of cells (12). ScRNAseq technology provides a powerful new way to characterize the clonal diversity of tumor cells and explore the role of rare cells in tumor development (13). ScRNA-seq has been used in a variety of human tumors, including glioma (14), lung cancer (15), hepatocellular carcinoma (16), nasopharyngeal carcinoma (17), colorectal cancer (18), ovarian cancer (19), gallbladder cancer (20), oral cancer (21), laryngeal carcinoma (22), etc. Chen et al. constructed the first single cell transcriptome map of hypopharyngeal cancer, revealed the complex crosstalk in HSCC and found BMPR promotes HSCC cells proliferation and migration (23).
The tumor microenvironment (TME) refers to the special environment in which tumor cells grow by interacting with the extracellular matrix during the growth process (24). Tumor immune environment (TIME) is an important part of TME, including tumor-associated macrophages (TAM), mast cells, T lymphocytes, B lymphocytes, and natural killer (NK) cells, myeloid-derived suppressor cells (MDSC), and other subgroups (25). Single-cell RNA sequencing can accurately identify different immune cell populations in the microenvironment based on the single-cell level and biomarkers that can be used to characterize such cells, providing the cellular composition and distribution characteristics of the tumor immune microenvironment from a holistic perspective, thereby revealing their functional states and discovering potential immunotherapy targets (26). Studying the tumor immune microenvironment by using single-cell sequencing has been applied to melanoma (27), liver cancer (28), lung cancer (29), breast cancer (30), colorectal cancer (31) and so on. However, the immune microenvironment in HSCC has not been explored.
In the present study, TIME in HSCC and lymphatic metastasis was revealed for the first time by performing scRNA-seq in hypopharyngeal carcinoma, paracancerous, and lymphoid tissues from HSCC patients, which may improve our current understanding of the mechanisms of HSCC development and progression.

Materials and methods
Sample collection and processing HSCC tissues, adjacent normal tissues and lymph tissues from five HSCC patients were collected and stored at -20°C. All patient data used in the study were approved by the Ethics Committee of Qilu Hospital of Shandong University. The patients participating in the program were informed.

Tissue dissociation and preparation
Fresh tumor tissue was stored in GEXSCOPE ® tissue preservation solution (Singleron) and shipped on ice to the Singleron laboratory as soon as possible. Samples were washed 3 times with Hanks' Balanced Salt Solution (HBSS) and cut into 1-2 mm pieces. The tissue sections were then digested with 2 ml of GEXSCOPE ® tissue dissociation solution (Singleron) and placed in 15 ml centrifuge tubes at 37°C with continuous agitation for 15 minutes. After digestion, the samples were filtered through a 40 mm sterile filter, centrifuged at 1000 RPM for 5 min, the supernatant was discarded, and the pellet was resuspended in 1 ml of PBS (HyClone). To remove erythrocytes, 2 ml of GEXSCOPE ® erythrocyte lysis buffer (Singleron) was added for 10 minutes at 25°C. Centrifuge at 500 × g for 5 min and suspend in PBS. Samples were stained with trypan blue (Sigma) and evaluated under a microscope. Cell activity and cell concentration were measured using a fluorescence cell analyzer (Countstar Rigel S2). If there were a lot of dead cells or debris, Biotec Dead Cell Removol Kit (Miltenyi Biotec) were used for removal.

Single cell RNA sequencing
A single-cell suspension was prepared at a concentration of 1 × 10 5 cells/mL in PBS (HyClone). The scRNA-Seq library was constructed by the GEXSCOPE ® Single Cell RNA Library Kit (Singleron Biotechnologies) according to the Singleron GEXSCOPE ® protocol (32). Individual libraries were diluted to 4 nM and pooled for sequencing. 150 bp paired-end sequencing was performed using Illumina HiSeq X.

ScRNA-seq quantifications and statistical analysis
The batch-effect was assessed and corrected using the Harmony algorithm. Raw reads are processed through an internal pipeline to generate gene expression profiles. Briefly, cell barcodes and UMIs were extracted after filtering reads in the absence of multiple t-tails. Read 2 will be trimmed with splicing and poly A tails (FASTP V1) and then aligned to GRCh38 using integrated version 92 gene annotation (FASTP 2.5.3A and featurecount 1.6.2) (33). Reads for the same cell barcode, UMI, and gene were combined to calculate the number of UMIs per gene per cell. The UMI count table for each cell barcode was used for further analysis. Cell type identification and cluster analysis were performed using the Seurat program (34,35). The RNA sequencing data were analyzed using the Seurat program (http://satijalab.org/seurat/, R package, v. 3.0.1). Loaded the UMI count table into R using read table function, then set the parameter resolution of the FindClusters function to 0.6 for cluster analysis. Used the findmarker function to identify differentially expressed genes (DEGs) between different samples or consecutive clusters. GO function enrichment analysis was performed on gene sets using ClusterProfiler software to find biological functions or pathways significantly associated with specific expressed genes (36).

T cell receptor library preparation and scRNA-seq
The single cell suspension (1×10 5 cells/mL) was loaded into microfluidic devices. Subsequently, the scTCR-seq libraries were constructed according to the protocol of GEXSCOPE Single Cell Immuno-TCR Kit (Singleron Biotechnologies). In brief, the magnetic beads with molecular labels captured the poly (A) tail and the T-cell receptor (TCR) region on the mRNA to label the cells and mRNA after the cells were lysed. Afterwards, the magnetic beads in the chip were collected and then mRNAs captured by the magnetic beads were reverse transcribed into complementary DNA (cDNA) and amplified. Sequencing libraries suitable for the Illumina sequencing platform were constructed after partial cDNA fragmentation and splicing. The remaining cDNA was enriched for the immune receptor (TCR), and the enriched products were amplified by PCR to construct a sequencing library suitable for the Illumina sequencing platform. Finally, each library was sequenced on Illumina HiSeq X with 150 bp paired-end reads.

TCR library analysis
TCR clonotypes assignment was performed using Cell Ranger (v4.0.0) vdj pipeline with GRCh38 as reference. In brief, a TCR diversity metric was obtained, which contains the frequency of clonotype and barcode information. Only cells with one productive TCR a-chain (TRA) and one productive TCR b-chain (TRB) were kept for further analysis. Each unique TRA(s)-TRB(s) pair was defined as a clonotype. If one clonotype was present in at least two cells, cells harboring this clonotype were considered to be clonal and the number of cells with such pairs indicated the degree of clonality of the clonotype. The TCR diversity index was calculated using the vegan package in R, with the shanno and invsimpson indices being computed through the diversity() function, while Chao and ACE were calculated using the estimateR() function. After computing the diversity for each sample based on the frequency of different clonotypes, plotted the boxplot by R package ggplot2.

Quality control, dimension-reduction and clustering
Cells with gene counts less than 200 or in the top 2% were excluded, and cells with UMI counts in the top 2% were excluded. Removed cells containing more than 20% mitochondria. Used functions in Seurat V3.1.2 for reduction and clustering (42). All gene expressions were normalized and scaled using NormalizeData () and ScaleData(). FindVariableFautres() selects the top 2000 variable genes for PCA analysis. FindClusters() divides the cells into 31 clusters using the top 20 principal components and a resolution parameter of 1.2. For subclustering of seven cell types, set the resolution to 0.8. Applied the Uniform Manifold Approximationand Projection (UMAP) algorithm to visualize cells in two-dimensional space.

Differentially expressed genes analysis
Seurat FindMarkers() selected genes expressed by more than 10% of the cells in the cluster with an average log (fold change) greater than 0.25 as DEGs based on the Wilcox likelihood ratio test with default parameters.

Cell type annotation
The cell type identity of each cluster was determined by the expression of canonical markers found in DEGs combined with the knowledge of the published literature (43). Seurat DoHeatmap ()/DotPlot()/Vlnplot() generated a heatmap/dotplot/violin plot showing the markers used to identify each cell type.

Pathway enrichment analysis
To investigate the potential functions of subdivided cells, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were used with the "clusterProfiler" R package version (36). Pathways with P_ADJ values less than 0.05 were considered significantly enriched. For GSVA pathway enrichment analysis, the average gene expression per cell type was used as input data using the GSVA software package (44). Reference Gene Ontology genes include Molecular Function (MF), Biological Process (BP), and Cellular Component (CC) categories. The protein-protein interactions (PPIs) of DEGs in different clusters were predicted based on the interactions of known genes with associated GO terms in StringDB (1.22.0) (45).

Trajectory analysis
To map the differentiation and transformation of cell subtypes in HSCC tissues, adjacent normal tissues, and lymphatic metastatic tissues, a spurious time trajectory analysis was performed using Monocle2 (46). To construct trajectories, using differentially expressed genes to rank cells in order of spatiotemporal differentiation, we performed FindVairable Features and dimensionality reduction by using DDRTree. Trajectories were visualized by plot_cell_trajectory().

RNA velocity
For RNA velocity, BAM files containing tumor or epithelial cells and the reference genome GRCh38 were used for analysis in Python with velocyto (v 0.2.3) and scVelo (v 0.17.17), default parameters. The results were projected onto the UMAP map from the Seurat cluster analysis for visual consistency.

Expression programs analysis
Transcription programs were extracted using the cNMF algorithm, the top 50 genes were used as meta-signatures, and the scores for each program in each cell were calculated based on the meta-signatures. The metaprogram performed computations and hierarchical clustering based on artificial correlations between each program.

UCell gene set scoring
Gene set scoring Gene set scoring was performed using the R package UCell v 1.1.0 (47). UCell scores are based on the Mann-Whitney U statistic by ranking query genes' in order of their expression levels in individual cells. Because UCell is a rank-based scoring method, it is suitable to be used in large datasets containing multiple samples and batches.

Immunohistochemistry staining and statistical analysis
Tissue sections were fixed, dehydrated, and antigenically repaired. At first, the first antibody, a specific antibody that recognizes the antigen was added, and then, the second antibody, a biotin-labeled antibody that recognizes the FC segment of the first antibody was added. After that, the lecitin, biotin, and horseradish peroxidase complex was added, which was generally configured in the first 30 minutes and finally displayed with the substrate of the enzyme. At last, Observed and photographed by using a microscope. IgA antibody (ImmunoWay, YT2281, 1:200) and IgG1 antibody (ImmunoWay, YT2293, 1:200) were used for IHC staining. SPSS 22.0 statistical software package was used for statistical analysis. The two-tailed Student's t test was used to assess the statistical differences between the groups. The data were consistent with normal distribution, and P<0.05 was considered statistically significant.

Single-cell RNA expression profiling in HSCC
To explore the cellular diversity and microenvironment composition of HSCC, we performed scRNA-seq and T cell receptor (TCR) analysis of primary HSCC, adjacent normal, and lymphoid tissues from five HSCC patients ( Figure 1A). After quality control assessment, we obtained transcriptomes of 132,869 single cells using the Singleron ™ Single-Cell mRNA Whole Transcriptome Analysis System, of which 52,145 cells were derived from primary tumor tissue, 39,757 cells were derived from lymph tissue cells, and 40,967 cells were derived from adjacent normal tissue. Tumor metastasis occurred in three lymphoid tissues while it did not occur in two lymphoid tissues. Seven distinct cell populations were identified from whole singlecell analysis based on t-distributed stochastic neighbor embedding (t-SNE) analysis and canonical marker expression, including lymphocytes, MPs, fibroblast cells, ECs, epithelial cells, pDCs and mast cells ( Figure 1B). These cell populations were unevenly distributed in different kinds of tissues ( Figures 1C; S1A). The proportion of lymphocytes was highest in lymphoid tissues, followed by HSCC tissue, and lowest in normal hypopharyngeal tissue. The proportion of MPs was highest in HSCC tissues. A dot plot of the top 5 differential genes in each cell subset were shown ( Figure 1D). These results reveal the different cell types distribution in normal hypopharyngeal tissues, HSCC tissue and tumor metastasis tissues.

Subpopulations and transcriptome landscape of lymphocytes and MPs in HSCC and lymphoid tissues
We then discovered the role of immune cells in HSCC, which mainly lymphocytes and MPs. Four distinct populations were identified from lymphocytes, including B cells, plasma cells, T cells and proliferation lymphocytes (Figure 2A). These cell populations distributed differently in different tissues (Figures 2B,  S1B). The stacked vin plot showed marker genes in each cell subset ( Figure 2C). UMAP plot showed six subtypes colored in HSCC tissues, adjacent normal tissues and lymphatic tissues ( Figure 2D). The heatmap of the top10 differential genes in each cell subset was shown ( Figure 2E). Five distinct populations were identified from MPs, including macrophages, monocytes, mature dendritic cells (DCs), cDC1 and cDC2 ( Figure 2F). These cell populations distributed differently in different tissues (Figures 2G, S1C). Stacked vin plot showed marker genes in each cell subset ( Figure 2H). UMAP plot showed five subtypes colored in HSCC tissues, adjacent normal tissues and lymphatic tissues ( Figure 2I). The heatmap of top10 differential genes in each cell subset was shown ( Figure 2J).  Figures 3B, S1D). IGHA1 and IGHG1 specific plasma cells were significantly overexpressed in HSCC tissues compared with normal hypopharygeal tissues. Stacked vin plot showed marker genes in each cell subset ( Figure 3C). UMAP plot showed seven subtypes colored in HSCC tissues, adjacent normal tissues and lymphatic tissues ( Figure 3D). The heatmap of the top10 differential genes in each cell subset was shown ( Figure 3E). Seven distinct populations were identified from T cells, including NK, Th2, Tfh, CD8Teff, Proliferating T cells, Naive T, and Treg cells ( Figures 3F). These cell populations distributed differently in different tissues ( Figures 3G, S1E). Treg cells were significantly overexpressed in HSCC tissues and lymphatic tissues compared with normal hypopharygeal tissues. Naive T cells were significantly overexpressed in lymphatic tissues compared with normal hypopharygeal tissues and HSCC tissues. Stacked vin plot showed marker genes in each cell subset ( Figure 3H). UMAP plot showed seven subtypes colored in HSCC tissues, adjacent normal tissues and lymphatic tissues ( Figure 3I). The heatmap of the top10 differential genes in each cell subset was shown ( Figure 3J). The volcano plots and the heatmaps showed different genes of IGHG1 plasma cells between HSCC tissues, adjacent normal tissues, and lymphatic tissues. IGLV1-40, SPP1, CXCL8, IGLV3-10, and COL1A1 were the most up-regulated genes in HSCC tissues compared with adjacent normal hypopharynx tissues (Figures 4A, B). IGKV2D-29, IGHA2, IGHV3-48, IGHV4-34, and IGKV3-11 were the most upregulated genes in lymphatic tissues compared with adjacent normal hypopharynx tissues ( Figures 4C, D). IGKV2D-29, IGHV4-34, IGHV3-48, IGKV3-11, and CD52 were the most up-regulated genes in lymphatic tissues compared with HSCC tissues (Figures 4E, F).

Subpopulations, pseudotime trajectory and transcriptome landscape of macrophages and monocytes in HSCC and lymphoid tissues
Six distinct populations were identified from macrophages cells ( Figure 5A). These cell populations distributed differently in different tissues ( Figures 5B, S1F). Macrophage1 was significantly overexpressed in HSCC tissues and lymphatic tissues compared with normal hypopharygeal tissues, with marker genes including CXCL5, SPP1, INHBA, MMP1, FN1, MMP12, TNFAIP6, CHI3L1 and CCL20. UMAP plot showed six subtypes colored in HSCC tissues, adjacent normal tissues and lymphatic tissues ( Figure 5C). Macrophage1 can also be called SPP1+ macrophages according to its marker genes. Heatmap of the top10 differential genes in each cell subset was shown ( Figure 5D). Pseudo-time of macrophages expression profiles was reconstituted ( Figure 5E). Molecular functions, cytological components, biological and KEGG pathways of macrophages1 were shown ( Figures 5E, F). We also detected the activation of multiple key regulators and TFs in six macrophages populations ( Figure 5G). The heatmaps and the volcano plots showed different genes of macrophage1 cells between HSCC tissues, adjacent normal tissues, and lymphatic tissues. CXCL5, MMP12, IGHG4, IGHG1, and CCL18 were the most up-regulated genes in HSCC tissues compared with adjacent normal hypopharynx tissues ( Figures 6A, B). SPP1, APOC1, FN1, IGHG3, and MMP12 were the most up-regulated genes in lymphatic tissues compared with adjacent normal hypopharynx tissues (Figures 6C, D). FN1, CD36, APOC1, CD52, and SPP1 were the most up-regulated genes in lymphatic tissues compared with HSCC tissues (Figures 6E, F). Four distinct monocytes populations were identified from monocytes cells ( Figure  S2A). These cell populations distributed differently in different tissues ( Figures S1G, S2B). UMAP plot showed four subtypes colored in HSCC tissues, adjacent normal tissues and lymphatic tissues ( Figure  S2C). Dotplot of the top5 differential genes and heatmap of the top10 differential genes in each cell subset were shown (Figures S2D, E).

Subpopulations and transcriptome landscape of proliferating lymphocytes and T cell exhaustion in HSCC and lymphoid tissues
Four distinct populations were identified from proliferating lymphocytes, including proliferating NK, proliferating B cells, proliferating plasma cells and proliferating T cells ( Figure 7A). These cell populations distributed differently in different tissues ( Figures 7B; S1H). Proliferating T cells and proliferating NK cells were significantly overexpressed in HSCC tissues and lymphatic tissues compared with normal hypopharygeal tissues. Proliferating B cells were significantly overexpressed in lymphatic tissues compared with HSCC tissues and normal hypopharygeal tissues. Bubble chart showing 5 typical genes expressed in each subtype ( Figure 7C). UMAP plot showed four subtypes colored in HSCC tissues, adjacent normal tissues and lymphatic tissues ( Figure 7D). The heatmap of the top10 differential genes in each cell subset was shown ( Figure 7E). Exhausted CD8+ effector T (Teff) cells levels and exhausted NK T cells levels in HSCC tissues, adjacent normal tissues and lymphatic tissues were measured (Figures 7F, G). Exhaustion of CD8+ Teff cells occurred in HSCC tissues.

Comprehensive analysis of the TCR repertoire in HSCC and lymphoid tissues
We then studied the TCR repertoires in HSCC, adjacent normal and lymphatic tissues. The clonotype ratio showed clonal diversity of T-cell repertoires in HSCC tissues, adjacent normal tissues, and lymphatic tissues ( Figure 8A). Clonal diversity in HSCC tissues, adjacent normal tissues, and lymphatic tissues was analyzed by using Shannon score, Inv. Simpson score, Chao score and ACE score (Figures S3A-D). The diversity of TCR clonal was lower in HSCC and lymphoid tissues than in normal hypopharynx tissues. VDJ TCR features and TOP10 colontypes were shown (Figures 8B, C). Most of TOP10 colontypes related to proliferating T cells and CD8 Teff cells. Clonotypes distribution of grouped clonotypes were analyzed ( Figure 8D). Treg cells showed significant higher clonal expansion in HSCC tissues compared with lymphoid metastasis and normal hypopharynx tissues. VJ pairs heatmap showed different VJ rearrangement in HSCC tissues, adjacent normal tissues and lymphatic tissues (Figures 8E-G).

IHC staining of IgA and IgG1 in HSCC tissues and adjacent normal tissues
At last, since IGHA2 was also upregulated in plasma cells of HSCC tissues, we verified IgA protein expression levels which encoded by IGHA1 and IGHA2 in HSCC tissues and adjacent normal tissues by using IHC staining. IgA positive stain in HSCC tissues was significantly higher than that in adjacent normal tissues ( Figures 9A, B). We also verified IgG1 protein expression levels which encoded by IGHG1 in HSCC tissues and adjacent normal tissues by using IHC staining. IgG1 positive stain in HSCC tissues was significantly higher than that in adjacent normal tissues ( Figures 9C, D). These results showed IgA and IgG1 may be potential diagnostic markers of HSCC.

Discussion
HSCC is a malignant disease with a poor prognosis. Over the past few decades, TIME targeting strategies have provided new therapeutic options for cancer therapy. However, these strategies have not yet been applied clinically in HSCC because the cellular characteristics and immune microenvironment of HSCC are largely unknown. In this study, we performed single-cell transcriptomic profiling of HSCC tissues, adjacent normal tissues, and lymphoid tissues to reveal TIME in HSCC and lymphatic metastasis for the first time. These results may improve our current understanding of HSCC development and progression, and provide new therapeutic targets for HSCC. TIME is an important component of TME, including T lymphocytes, B lymphocytes, NK cells, macrophages and other cell subsets; these immune cells not only play a role in killing tumor cells, for example, CD8+ T cells, natural killer cells, and M1 macrophages; but also promote tumor development, for example, Th2 cells, Treg cells, and M2 macrophages (48). Therefore, the tumor microenvironment has become a potential anticancer target, and its research has become a hot spot in tumor biology (49). Single-cell sequencing technology is a powerful tool for analyzing the heterogeneity of cellular components in the TME, and scRNAseq analysis can detect more diverse TME immune cells in tumor tissues, thus serving as a highly feasible platform for analyzing the tumor microenvironment (50).
In this study, we studied the immune microenvironment in HSCC and lymphatic tissues. Seven distinct cell populations were identified from the whole single-cell analysis, including lymphocytes, MPs, fibroblasts cells, ECs, epithelial cells, pDCs, and mast cells. We focused on lymphocytes and MPs in HSCC tissues, adjacent normal tissues, and lymphatic tissues. Four distinct populations were identified from lymphocytes, including B cells, plasma cells, T cells and proliferation lymphocytes. Subdivide of B cells and T cells and their distribution in three kinds of tissues as well as marker genes were analyzed. Four distinct populations from MPs were identified, including macrophages, monocytes, cDC2, MatureDCs, and cDC1. Subdivide of macrophages cells and monocytes cells and their distribution in three kinds of tissues as well as marker genes were analyzed. We also performed TCR repertoire analysis, which including clonal diversity, clonotype distribution and V-J pairing in HSCC tissues, adjacent normal tissues, and lymphatic tissues. We firstly verified that IGHA1 and IGHA2 up-regulated at mRNA level, and IgA up-regulated at protein level in HSCC by using single-cell sequencing and IHC staining. The basic structure of immunoglobulin is a monomer composed of four symmetrical polypeptide chains, including two identical heavy chains with larger molecular weight and two identical light chains with smaller molecular weight, and there are two light chains between the light and heavy chains. Immunoglobulin A is divided into two subclasses, IGHA1 and IGHA2, which are composed of the heavy chain a1 or a2 and the light chain respectively, and its heavy chain constant region is located on chromosome 14q32.33. IGHA1 was upregulated and reported as an unfavorable biomarker in renal cell carcinoma (51) and prostate cancer (52), while reported as a favorable biomarker in breast cancer (53). We firstly verified that IGHG1 up-regulated at mRNA level and IgG1 up-regulated at protein level in HSCC by using single-cell sequencing and IHC staining. IGHG1 is immunoglobulin gamma-1 heavy chain constant region, belongs to immunoglobulin G (IgG), which accounts for approximately 80% of total immunoglobulins (54). Although only B cells and plasma cells are known to produce IgG, there are increasing reports suggesting that IgG can be produced by several malignant cells, such as tumor cells from the esophagus, breast, liver, or prostate (55). IGHG1 promotes tumor development in gastric cancer, breast cancer and prostate cancer via AKT and MEK pathway (55)(56)(57). Overall, IGHA1 and IGHG1 may be potential diagnostic markers and serum markers of HSCC. TAMs play an important role in immunosuppressive macrophages. TAMs occupy a large proportion of tumor mesenchymal cells, most of which are migrated and differentiated from peripheral blood mononuclear cells and developed under the influence of tumor cells and their microenvironment. TAMs can secrete a variety of growth factors, cytokines, immunosuppressive mediators and proteolytic enzymes to promote tumor progression and metastasis (58). We found Macrophage1 was significantly overexpressed in HSCC tissues and lymphatic tissues compared with normal hypopharygeal tissues. Macrogphage1 may be potential tumor-associated macrophages (TAMs) in HSCC. Macrophage1 can also be called SPP1+ macrophages according to its marker genes. According to the previous literatures, the overexpression of SPP1 in macrophages is thought to correlate with the phenotype of M2-type macrophages (23). The number of M2-type macrophages was negatively correlated with progression-free survival, distant metastasis-free survival and overall survival of HNSCC patients (59). M2-type macrophages may play an important role in assisting HNSCC tumors to evade immune surveillance, and promoting tumor invasion and metastasis in HNSCC (60).
We measured exhausted CD8+ Teff cells levels and exhausted NK T cells levels in HSCC tissues, adjacent normal tissues and lymphatic tissues. We found exhaustion of CD8+ Teff occurred in HSCC tissues. CD8+ T cells exhaustion is an important factor that affecting tumor progression, its heterogeneity is closely related to the strength of immunotherapy response and survival prognosis of tumor patients (61).
In conclusion, we investigated TIME in HSCC and lymphatic metastasis by performing single-cell RNA sequencing in hypopharyngeal carcinoma, paracancerous tissues and lymphatic tissues of five HSCC patients. Different cell populations and subpopulations and their marker genes, and differentially expressed genes were identified. These results may not only elucidate the mechanism of occurrence and development of HSCC, but also provide potential therapeutic targets for HSCC.

Data availability statement
The data presented in the study are deposited in the GEO repository, accession number GSE206038, GSE227156. The analysis of raw data will be shared on reasonable request to the corresponding author.

Ethics statement
The studies involving human participants were reviewed and approved by Ethics Committee of Qilu Hospital of Shandong University. The patients/participants provided their written informed consent to participate in this study.

Author contributions
DL contributed to the conception of the study; CL performed the study and wrote the manuscript; FC and PW helped in the early design. LC acquired HSCC tissues, RG, WL, DW, SC and CX helped perform the analysis with constructive discussions. All authors contributed to the article and approved the submitted version.