The GSE150321 dataset was downloaded from gene expression omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) including LSCC01 and LSCC02 tumor samples obtained from 66-year-old patient and a 64-year-old male patient. Neither patient received clinical treatment.

The Seurat package11 (version 4.1.0) in R language (version 4.1.2) was performed to read the single cell gene expression matrix, based on which we conducted the following analysis. First, we delete the cells with <200 expressed genes and genes with <3 cells present. The PercentageFeatureSet function was performed to caLSCCulate the propertion of mitochondrial genes in cells. Next, the cells meet the threshold of UMI ranged from 200 to 6000, mitochondrial gene ratio <50% were included in this study. Then, the data was normalized using SCTransform function, analyzed for principal component analysis using RunPCA function, and removed batch effects using harmony package.12 After that, according to the top 30 principal components, the FindNeighbors function was carried out to construct the public nearest neighbor graph based on euclidean distance metric. The clustering cells were visualized using RunUMAP function, which was divided using FindClusters function. The resolution set to 0.26 for all the cells and epithelial cells, and 0.1 for T-cell. Finally, we performed DoubletFinder package13 to remove double cells. Subsequently, a total of 711 double cells were removed and the rest 11,990 cells were included in our study for further analysis.

To identify the Differentially Expressed Genes (DEGs) in cell subsets, the FindAllMarker function in Seurat package was performed to caLSCCulate the significance of the DEGs with the parameter of min.pct = 0.25, logfc.threshold = 0.5, and only.pos = TRUE. Then, to further understand the function of DEGs, we searched the biological process of DEGs from DAVID database14 (https://david.ncifcrf.gov/) using the threshold of p < 0.05.

The monocle2 in R package15 (version 2.22.0) was used to performed the pseudo-time trajectory analysis of single cells. The data processed using Seurat was regarded as the gene-cell matrix of raw UMI counts. Then, we used the newCellDaraset function to produce an object for pseudo-time trajectory analysis (expressionFamily = negbinomial.size). Subsequently, we included genes with the mean expression value was above 0.1 and their expression detected in at least 10 cells for pseudo-time trajectory analysis. Then, the differential GeneTest function was used to screen the DEGs, and the genes with q-value <1e-30 was used to dimensionality reduction. Next, the reduce Dimension function was used for dimensionality reduction with the parameter of reduction_method = “DDRTree” and max_components = 2. We used the orderCells function to construct pseudotime trajectories and arrange cells, and plot_cell_trajectory to sort and visualize cells. After that, differential GeneTest with the parameter of fullModelFormulaStr = “∼sm.ns(Pseudotime)” was performed to caLSCCulate the genes changed with pseudotime. The genes were visualized using “plot_pseudotime_heatmap”, which was grouped based on expression pattern. For branch analysis of single-cell differentiation trajectories, we used the BEAM function to count genes with significantly different expression patterns between the two branches. Also, the genes were visualized using “plot_genes_branched_heatmap”, which were grouped according to expression patterns.

CelLSCChat16 is an analytical tool for exploring ligand–receptor interactions of specific signaling pathways between cells. To explore the ligand–receptor interactions between different types of epithelial cells and T-cells, the CelLSCChat in R package (Version 1.1.3) was used to perform the analysis of cell–cell communication.

All analysis procedures and image generation were performed by R language (version 4.1.2), and p < 0.05 was considered to be statistically significant.

Through standardization by SCTransform, UMAP dimensionality reduction and clustering and removal of batch effects by harmony, we found that the cells in LSCC001 and LSCC002 were evenly distributed together (Fig. 1A). Next, all the cells were divided into 16 subpopulation that included 7 epithelial cell subsets, 3 T-cell subsets (Fig. 1B). Besides, the 7 epithelial cell subsets were found related to their high expression of epithelial cell markers, such as KRT5, KRT14 and KRT19. T-cell subsets have high expression of the marker gene CD3D. Meanwhile, regulatory T-cells were correlated to high expression of CTLA4, TIGIT and other genes that inhibit the activity of immune cells, proliferating T-cells highly expressed TOP2A, MKI67 and other cell proliferation-related genes. These results indicated that the T-cell subset was in an active state of cell division. We have listed the high expression genes in each cell subset as Table S1.

Figure 1.

Grouping and identification of cells in laryngeal cancer tissues. (A) Distribution of cells from different samples. (B) Location distribution of UMAP in different types of cells. (C) Bubble plot of marker genes in different cell types.

(0.42MB).

To explore heterogeneity and differentiation trajectories of epithelial cells in LSCC, we merged all the epithelial cells and re-clustered them into a total of 10 cell subsets (Fig. 2). We listed the genes with high expression in epithelial cells as Table S2. Then, we studied the biological process of the highly expressed genes in each epithelial cell. Interestingly, we found that epithelial cell subsets 1, 3, 6 were mainly involved in cell cycle and cell division related functions. Epithelial cell subsets 2, 4, 5, 7 and 10 were mainly involved in immune-related functions such as immune response, T-cell activation, and inflammatory response. Epithelial cell subsets 8 and 9 were mainly involved in cell migration-related functions (Fig. 2B).

Figure 2.

Heterogeneity and differentiation trajectories of epithelial cells in laryngeal carcinoma tissue. (A) UMAP plot of epithelial cells divided into 10 subpopulations. (B) Bubble plot of the enrichment of highly expressed genes in 10 epithelial cell subpopulations to representative biological processes. (C) Differentiation trajectories of 10 epithelial cell subsets and the proportion of each type of epithelial cell in each state. (D) Mountain plot of density distribution of 10 epithelial cell subsets in pseudotime. (E) Variation with pseudotime the biological process to which the genes are mainly enriched. (F) The biological process to which the branch-related genes between State 2 and State 3 are mainly enriched.

(0.77MB).

Moreover, we further explored the differentiation trajectory of the epithelial cells, the results indicated that the subgroup 3 of epithelial cells in the branch of State 2 accounted for the largest proportion, while the subgroups 1, 2, 4 and 5 of epithelial cells in the branch of State 3 accounted for the largest proportion. This result implies that epithelial cells are mainly differentiated into two types in the process of differentiation and evolution, one is mainly involved in cell proliferation, and the other is mainly involved in immune activation (Fig. 2C‒D). Besides, with the progress of epithelial cell differentiation, the ability of cell division and migration, and the activation of the immune system gradually increased, while the ability of apoptosis, cytoplasmic translation, and glucose metabolic process gradually decreased (Fig. 2E). Finally, we explored the differences in gene expression between the two branches of State 2 and State 3, the results indicated that with the pseudotime, the expression of genes related to cell division and mitotic cell cycle in State 2 gradually increased, while the expression of genes associated with immune response and inflammatory response in State 3 also gradually increased (Fig. 2F).

To explore heterogeneity and differentiation trajectories of NK/T-cells in LSCC, we merged all the epithelial cells and re-clustered them into 4 cell subsets (Fig. 3A). We found that the NK/T-cell subsets 1 and 3 had high expression of GZMA, GZMB, NKG7 and other cytotoxic T-cell marker genes, while NK/T-cell subset 2 were highly expressed CTLA4, TIGIT and other exhausted T-cell marker genes (Fig. 3B). We have listed the gene markers of NK/T-cell at Table S3.

Figure 3.

Heterogeneity and differentiation trajectories of T-cells in laryngeal cancer tissues. (A) UMAP plot of T-cells divided into 4 subpopulations. (B) Bubble plot of known different T-cell marker genes in 4 subpopulations. (C) Differentiation trajectories of the 4 T-cell subsets and the proportion of each T-cell in each state. (D) Mountain plot of the density distribution of the 4 T-cell subsets in pseudotime. (E) Genes as a function of pseudotime the biological process that is mainly enriched. (F) The biological process that the branch-related genes between State 2 and State 3 are mainly enriched.

(0.75MB).

Next, we explored the differentiation trajectories of NK/T-cells, the results indicated that NK/T-cell subsets 1, 2, 3 in State 2 had the largest proportion, and the NK/T-cell subset 2 in State 3 had the biggest proportion (Fig. 3C and D), indicating that during the differentiation and evolution of NK/T-cells, one type tends to activate immune cells and kill tumor cells, and the other type tends to be exhausted. Then we found that with the progress of NK/T-cell differentiation, the ability of cell division and aerobic respiration gradually increased, and the ability of T-cell activation gradually decreased (Fig. 3E). Furthermore, we analyzed the differential expression of genes between State 2 and State 3, the results found that with the change of pseudotime, the expression of genes related to cell cycle, oxidative phosphorylation in State 2 were gradually increased, while the expression of genes related to immune response, inflammatory response, T-cell activation were also improved (Fig. 3F).

Finally, we explored the ligand–receptor interactions of different types of epithelial cells and NK/T-cells in LSCC tissues. In particular, we pooled epithelial cell subsets 1, 3, 6 together and named proliferative epithelial cells; epithelial cell subsets 2, 4, 5, 7, 10 pooled together and named immune-related epithelial cells; epithelial cells Cell subsets 8, 9 were merged together and named migration-related epithelial cells. For NK/T-cells, we pooled NK/T cell subsets 1 and 3 together and named cytotoxic T-cells; NK/T-cell subset 2 named exhausted T-cells.

Interestingly, we found that immune related and proliferative epithelial cells activate cytotoxic T-cells by acting on CD8A through HLA-A, B, C. Migration-related epithelial cells act on all T cells through ligand–receptor pairs such as CNTN2-CNTN2, CDH1-KLRG1, CDH1-CDH1. All epithelial cells had stronger CD86-CTLA4 and CD80-CTLA4 interactions with exhausted T-cells than with cytotoxic T-cells (Fig. 4A).

Figure 4.

Communication analysis of cell contact types between epithelial cells and immune cells in laryngeal cancer tissues. (A) Ligand–receptor bubble diagram of different types of epithelial cells acting on T-cells. (B) Different types of T-cells acting on epithelial cells ligand–receptor bubble diagram.

(0.65MB).

Conversely, all T-cells act on immune-related epithelial cells through ligand–receptor pairs such as HLA-C-KIR2DL3, HLA-B-CD8A, and HLA-E-CD8A. All T-cells act on migration-related epithelial cells through ligand–receptor pairs such as CNTN2-L1CAM, CNTN2-CNTN2, CDH1-KLRG1, CDH1-CDH1, CD99-CD99. Exhausted T-cells act on all epithelial cells through ligand–receptor pairs such as CD86-CTLA4, CD86-CD28, CD80-CTLA4, CD80-CD28 (Fig. 4B).