Germinal center B-cell subgroups in the tumor microenvironment cannot be overlooked: Their involvement in prognosis, immunotherapy response, and treatment resistance in head and neck squamous carcinoma

Abstract

Background: More than 60 % of patients with head and neck squamous carcinoma (HNSCC) are diagnosed at advanced stages and miss radical treatment. This has prompted the need to find new biomarkers to achieve early diagnosis and predict early recurrence and metastasis of tumors.

Methods: Single-cell RNA sequencing (scRNA-seq) data from HNSCC tissues and peripheral blood samples were obtained through the Gene Expression Omnibus (GEO) database (GSE164690) to characterize the B-cell subgroups, differentiation trajectories, and intercellular communication networks in HNSCC and to construct a prognostic model of the associated risks. In addition, this study analyzed the differences in clinical features, immune cell infiltration, functional enrichment, tumor mutational burden (TMB), and drug sensitivity between the high- and low-risk groups.

Results: Using scRNA-seq of HNSCC, we classified B and plasma cells into a total of four subgroups: naive B cells (NBs), germinal center B cells (GCBs), memory B cells (MBs), and plasma cells (PCs). Pseudotemporal trajectory analysis revealed that NBs and GCBs were at the early stage of B cell differentiation, while MBs and PCs were at the end. Cellular communication revealed that GCBs acted on tumor cells through the CD99 and SEMA4 signaling pathways. The independent prognostic value, immune cell infiltration, TMB and drug sensitivity assays were validated for the MEF2B⁺ GCB score groups.

Conclusions: We identified GCBs as B cell-specific prognostic biomarkers for the first time. The MEF2B⁺ GCB score fills the research gap in the genetic prognostic prediction model of HNSCC and is expected to provide a theoretical basis for finding new therapeutic targets for HNSCC.

Keywords: Drug sensitivity; Germinal center B cells; Head and neck squamous carcinoma; Immune cell infiltration; Single-cell RNA-Seq.

Fig. 1

| Integration and analysis of scRNAseq data of B and plasma cells in HNSCC. (A–B) The Uniform Manifold Approximation and Projection (UMAP) of the expression profiles of the 7082 B and plasma cells derived from tumor tissues and peripheral blood (PBL) of HNSCC. B cells are classified into 16 distinct transcriptional clusters. HPV infection status (HPV⁺ vs HPV⁻), Phase (G1, S, G2M). (C) The distribution of CNVscore, nCount_RNA, S.Score, G2M.Scorce for each subgroup was visualized using UMAP. (D) Pie chart showing the proportions of each subgroup in the HNSCC samples. (E) The proportion of subgroups derived from PBL and tumor tissues. (F) Visualization of sample characteristics associated with different B and plasma cell subgroups and different HPV infection status (HPV⁺ vs HPV⁻), Phase (G1, S, G2M), G2M.Score, S.Score, nCount_RNA, CNVscore, differential gene expression. (G) Volcano plots showing DEGs in subgroups.

Fig. 2

| Trajectory of B and plasma cells. (A) The differentiation trajectory of B and plasma cells shows the derivation of HNSCC cells, based on the pseudotime of cells (from left to right). (B) Pseudotime differentiation trajectory plots of B and plasma cell subgroups. (C) The UMAP visualization delineates the pseudotime scores of the entire spectrum of cells. The violin plot illustrates the pseudotime scores attributed to each distinct subgroup. The ridge map provides an insight into the density distribution of pseudotime values. (D–E) Bar charts of the changes in the proportion of B and plasma cell subgroups in the pseudotime state, respectively.

Fig. 3

| Landscape of NBs subgroups in the HNSCC. (A) The distribution of UMAP for the 6 subgroups in NBs. (B) Visualization of sample characteristics associated with different NBs subgroups and different HPV infection status (HPV⁺ vs HPV⁻), Phase (G1, S, G2M), G2M.Scorce, S.Score, nCount_RNA, CNVscore, differential gene expression. (C) Expression of top5 marker genes in different subgroups. (D) The distribution of CNVscore, nCount_RNA, S.Score, G2M.Scorce for each subgroup was visualized using UMAP. (E) Results of GO-BP enrichment analysis of differential genes and their corresponding transcription factors and surface proteins in 6 NBs subgroups. (F) Histogram of HPV infection status of NBs subgroups. (G) The cell cycle of NBs subgroups. (H–I) Violin plots of GO-BP Negative regulation of B cell activation and CNV score. (J–L) GSEA enrichment analysis among NBs subgroups.

Fig. 4

| Landscape of GCBs subgroups in the tumor microenvironment of HNSCC. (A) 332 GCBs were subjected to UMAP dimension reduction, divided into two cell clusters including C11 MEF2B⁺ Germinal center B cells (C11 GCBs, 184) and C14 STMN1⁺ Germinal center B cells (C14 GCBs, 148). The contribution of pie chart to patients (13 patients), HPV infection status (HPV⁺ vs HPV-) and Phase (G1, S, G2M) were visualized. (B–C) The bubble chart and UMAP showed the expression of GCBs top marker genes. (D) Visualization of relevant sample features for 332 GCBs using UMAP dimensionality reduction clustering: nCount_RNA, CNVscore, G2M.Scorce, S.Score. (E) The violin plot shows the nCount_RNA, CNVscore, G2M.Scorce, S.Score of C11 GCBs and C14 GCBs and their differences. ∗P < 0.05, ∗∗P < 0.01, ∗∗∗P < 0.001, means significant difference, ns means no significant difference. The boxplots showed the median (middle line), upper and lower quantiles (boxplot), and range of the data (violin plot). (F) nCount_RNA, CNVscore, G2M.Scorce, S.Score of GCBs with different HPV infection status (HPV⁺ vs HPV-) and their differences. ∗∗P < 0.05, ∗∗P < 0.01, ∗∗∗P < 0.001, means significant difference, ns means no significant difference. The boxplots showed the median (middle line), upper and lower quantiles (boxplot), and range of the data (violin plot). (G) The proportion of C11 GCBs and C14 GCBs in PBL or tumor group among patients (above). The ratio of C11 GCBs and C14 GCBs in different groups of HPV infection status (HPV⁺ vs HPV-) among patients (below). (H) The C11 GCBs and C14 GCBs ratios between the HPV⁺ and HPV-groups. (I–J) Bar graphs depicted the difference in cell proportion between subgroups and phases, as well as the difference in cell proportion between different HPV infection statuses and phases. (K) GSEA enrichment analysis of GCBs subgroups. The DEGs of C11 GCBs and C14 GCBs were screened, and the enrichment items of GOBP were scored by GSEA. NES >0 was up-regulated and <0 was down-regulated. NES, N stands for standardization and ES for enrichment scores. (L) Volcano plots showing expression of DEGs in GCBs subgroups. (M) The heat map showed the active biological processes of each subgroup of GCBs. (N) Word cloud for GO-BP enrichment analysis of GCBs DEG.

Fig. 5

| Landscape of MBs subgroups in the HNSCC microenvironment. (A) 3253 MBs were subjected to UMAP dimension reduction, divided into 6 cell clusters including C1 CRIP1⁺ Memory B cells(C1 MBs, 868), C2 LINC01781⁺ Memory B cells(C2 MBs, 738), C3 PPP1R15A⁺ Memory B cells(C3 MBs, 620), C4 IFI30⁺ Memory B cells(C4 MBs, 614), C9 CD1C⁺ Memory B cells(C9 MBs, 258), C13 TLE5⁺ Memory B cells(C13 MBs, 155). The contribution of pie chart to HPV infection status (HPV⁺ vs HPV⁻) and Phase (G1, S, G2M) were visualized. (B–C) The bubble chart and UMAP showed the expression of MBs top5 marker genes. (D) Visualization of relevant sample features for 3253 MBs using UMAP dimensionality reduction clustering: nCount_RNA, CNVscore, G2M.Scorce, S.Score. (E–F) The histogram showed the distribution of HPV infection status (HPV⁺ vs HPV⁻) and Phase (G1, G2M, S) in each subgroup of MBs intuitively. (G–H) The violin plot shows the nCount_RNA, CNVscore, G2M.Scorce, S.Score of MBs subgroups and their differences. ∗∗P < 0.05, ∗∗P < 0.01, ∗∗∗P < 0.001, means significant difference, ns means no significant difference. The boxplots showed the median (middle line), upper and lower quantiles (boxplot), and range of the data (violin plot). (I) The activity of different pathways in MBs subgroups was demonstrated by GO-BP enrichment analysis and word cloud. (J) The correlation of DEGs among MBs subgroups in positive regulation of cell communication, positive regulation of cell death, positive regulation of cellular metabolic process was analyzed by GSEA enrichment analysis.

Fig. 6

| Landscape of PCs subgroups in the HNSCC microenvironment. (A) The distribution of UMAP for the 2 subgroups in PCs. (B) UMAP was used to illustrate the distribution of CNVscore, nCount_RNA, S.Score, and G2M.Score for each subgroup. (C) Visualization of sample characteristics associated with different PCs and different HPV infection status (HPV⁺ vs HPV⁻): Phase (G1, S, G2M), G2M.Scorce, S.Score, nCount_RNA, CNVscore, differential gene expression. (D) Expression of marker genes in different PCs. (E–G) The histogram displayed the differences in cell percentage across different HPV groups and phase groups. (H–I) Results of GO-BP enrichment analysis of differential genes and their corresponding transcription factors and surface proteins in 2 PCs.

Fig. 7

| Cell communication network. (A) Circle visualization of interaction strength (left) and number (right) among the B and plasma cell subgroups and tumor cells. (B) Cell and communication patterns. (C) Illustration of the incoming and outgoing interaction strengths for each of the B and plasma cell subgroups and tumor cells. (D) Incoming communication patterns of target cells (left). Outgoing communication patterns of secreting cells (right). (E) Heatmap shows an overview of cell-cell interactions in both outgoing and incoming signaling patterns. (F) Interaction bubble diagram among GCBs, MBs and tumor cells. (G) The circle plots show the interaction strength (left) and number (right) of B and plasma cell subgroups on tumor cells. (H) Circle visualization of interaction strength (left) and number (right) about tumor cells on B cell subgroups.

Fig. 8

| GCBs based regulatory network detect tumor-specifc signaling pathways in HNSCC. (A, F) Illustration of the incoming and outgoing interaction strengths for CD99 and SEMA4 signaling pathways. (B, G) Heatmaps show dominant senders, receivers, mediators and infuencers in CD99 and SEMA4 signaling pathways of C11 GCBs inferred by network centrality score. (C, H) The violin diagram depicts the gene expression levels in the CD99 and SEMA4 signaling pathway networks including B and plasma cells as well as malignant cells. (D) The circle diagram shows the interaction of CD99 ligand on MBs and tumor cells. (I) The circle plot shows the interaction of SEMA4 ligand on GCBs and tumor cells. (E, J) The hierarchy plots of signaling pathway networks indicate the regulation of CD99 and SEMA4 on B cells and tumor cells in HNSCC, as predicted by CellChat anlysis.

Fig. 9

| Screening for prognosis-related memory B cells in HNSCC patients. (A, F) Univariate Cox regression analysis was performed to identify the prognosis-related genes. (B, G) Lasso regression analysis of prognosis-related genes. (C, H) Multivariate COX regression analysis on OS in the TCGA cohort. (HR > 1 protective factor, <1 risk factor.) (D, I) Kaplan–Meier curves were created to estimate OS for high- and low-score groups from TCGA cohort (P < 0.0001). (E, J) Kaplan-Meier curve showed that the prognosis of HNSCC with high expression of ATP6V0E1, PGK1, TPM3, ADK, LGALS1, TKT was poor in OS, while that of HNSCC with high expression of MS4A1, TRAC, BLK, ITGB7 was better in OS.

Fig. 10

| Screening for prognosis-related germinal center B cells in HNSCC patients. (A, F) Univariate Cox regression analysis was performed to identify the prognosis-related genes. (B, G) Lasso regression analysis of prognosis-related genes. (C, H) Multivariate COX regression analysis on OS in the TCGA cohort. (HR > 1 protective factor, <1 risk factor.). (D, I) Kaplan–Meier curves were created to estimate OS for high- and low-score groups from TCGA cohort (P < 0.0001). (E, J, K) Kaplan-Meier curve showed that the prognosis of HNSCC with high expression of HMMR, HSP90AA1, PRDX6, ACTB, BASP1, EZR, RRAS2 was poor in OS, while that of HNSCC with high expression of LRMP, RGS13, CD38, LCK, SMIM14, TCL1A was better in OS. (L) The bar diagrams show the distribution of the expression of CD38, SMIM14 and TCL1A in N0, N1, N2, N3, NX.

Fig. 11

MEF2B+GCB score. (A) Visualization of genes constituted the MEF2B⁺ GCB score in different GCB subgroups and different HPV infection statuses (HPV⁺ vs. HPV-), Phase (G1, S, G2M), and G2M.Scorce, S.Score, nCount_RNA, CNVscore, differential gene expression. (B–C) Expression of 18 prognosis-related genes in C11 and C14 subgroups. (D) Distribution and median value of risk scores of HNSCC patients. Distribution of OS status and risk scores of HNSCC patients. The heatmap showed the expression of 18 prognosis-related genes in high and low-risk groups. (E) Pie charts depict the differences in clinical characteristics (status, stage, TNM, phase) between high and low-risk groups in HNSCC. (F) The bar chart shows the risk coefficients that make up the risk score genes. (G) Time-dependent ROC curves and AUC values evaluate the prognostic performance of the risk score model at 1, 3, and 5 years. (H) An overview of the correlation between 18 prognosis-related genes. Red indicates a positive correlation, and blue indicates a negative correlation. (I) ACTB, BASP1, EZR, PRDX6, and RRAS2 were positively correlated with risk scores. CD38, LCK, RGS13, SMIM14, and TCL1A were negatively correlated with risk score (P < 0.001). (J) The differential expression of ACTB, BASP1, EZR, PRDX6, and RRAS2 in high and low-risk groups. (K) ACTB, BASP1, EZR, PRDX6, RRAS2, and risk score were negatively correlated with OS. (L–M) The correlation analysis among ACTB, BASP1, EZR, PRDX6, and RRAS2. (N) The differences of EZR and RRAS2 in high and low risk groups, different age groups, different races and different TNM stages.

Fig. 12

| Construction and verification of Nomogram and analysis of immune cell infiltration. (A) Multivariate analyses of the prognosis model. (B) Nomogram plot constructed by combining the risk score and clinical features. (C–D) AUC values evaluate the prognostic performance of the risk score model at 1, 3 and 5 years. (E) Deviations of actual 1-year, 3-year, and 5-year survival probabilities. (F)A heat map of principal component analysis for the immune cell infiltration between high and low risk groups. (G–H) Differential analysis of estimate score, immune score, stromal score, and tumor purity between high and low risk groups. (I)A boxplot of immune cells in the two risk subgroups. (J)Differential analysis of immune cell infiltration distribution between high and low risk groups. (K)The correlation between risk score and immune cells. The horizontal axis represents the correlation with the risk score, the vertical axis represents the proportion of immune cells, and the size of the circle represents the absolute value of the correlation coefficient. The larger the circle, the higher the correlation. (L)The correlation between prognostic genes and immune cells. Red indicates positive correlation, blue indicates negative correlation, ∗P < 0.05, ∗∗P < 0.01. (M)An overview of the correlation of immune cells. (N)The distribution of prognostic genes, estimate score, immune score, stromal score, tumor purity and immune cells in high and low risk groups. (O) RRAS2 was positively correlated with NBs, MBs, PB. (P) RRAS2 was negatively correlated with stromal score, immune score and tumor purity. (Q) Analysis of the correlation between RRAS2 and immune cells.

Fig. 13

| (A–B) DEGs analysis between high and low risk groups. (C–E) Functional enrichment analysis between the high and low risk groups of HNSCC patients. (C) Results of the GO enrichment analysis. (D) Results of the KEGG enrichment analysis. (E) Enrichment plots from the gene set enrichment analysis. Biological processes enriched in the high-risk group of HNSCC patients. Signaling pathways enriched in the low-risk HNSCC samples. (F) The waterfall map showed the distribution of the top 20 genes of TMB in the high and low risk groups. The horizontal coordinates represent the different samples and the vertical coordinates represent the genes. The different colors in the note below represent different types of mutations, the upper bars represent the probability of mutation in each sample, and the right bars represent the proportion of different types of mutations in each gene. (G) The CNV amplifcation, deletion or absence of data of risk genes. (H) Analysis of the correlation between gene and gene mutation. Mutually exclusive: if one gene is mutated, there is a high probability that no mutation will occur in the corresponding gene. Co−occurrence: two genes mutate at the same time. (I) The sample numbers, type, and location of mutations in each gene. The figure notes below represent the different mutation types. The dots in the figure represent the samples in which the mutation occurred, and the connected vertical lines indicate the location of the mutation. (J) TMB in high and low risk groups. (K) Correlation analysis between risk scores and TMB (P < 0.001). (L) Survival analysis of high risk-high TMB, high risk-low TMB, low risk-high TMB, low risk-low TMB. (M) AUC values evaluate the prognostic performance of the TMB model at 1, 3 and 5 years. (N–P) Chemotherapy and targeted drugs sensitivity of high- and low-risk subgroups.

Fig. 14

| PLXNB2 significantly affects the proliferation and migration of head and neck cancer cell lines. (A, B) CCK-8 experiment. After PLXNB2 knockdown, the proliferation ability of HN-5 and UMSCC-47 cell lines decreased significantly. (C, D) Plate cloning experiment. After PLXNB2 knockdown, the colony-forming ability of HN-5 and UMSCC-47 cell lines was significantly decreased. (E, F) wound healing test. After PLXNB2 knockdown, the migration ability of HN-5 and UMSCC-47 cell lines decreased significantly. (G, H) Transwell experiment. After PLXNB2 knockdown, the migration and invasion ability of HN-5 and UMSCC-47 cell lines were significantly reduced. (∗P < 0.05; ∗∗P < 0.01; ∗∗∗P < 0.001).

Integration and analysis of scRNAseq data of HNSCC cells. (A-D) The UMAP of the expression profiles of the 116212 HNSCC cells derived from tumor tissues and PBL. Cells are classified into 10 distinct transcriptional clusters. (E) The expression of marker genes in these cell types was visualized by UMAP plots.

The distribution of CNV. (A) The distribution of CNV in epithelial/tumor cells. (B) The distribution of CNV in the B and plasma cell subgroups.

The tumor microenvironment feature was captured by mIHC in three HNSCC samples. (A-C) B cells were stained with anti-CD79a, and CD8⁺ T cells were stained with anti-CD8, and the tumor epithelium was stained with anti-EPCAM. Scale bar, 200um.