Identification of a novel immune gene panel in tongue squamous cell carcinoma

Abstract

Background: Tongue squamous cell carcinoma (TSCC) is one of the most common oral cancers. Immune activity is significantly related to the initiation and progression of TSCC. Systemic analysis of the immunogenomic landscape and identification of crucial immune-related genes (IRGs) would help understanding of TSCC. Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) provide multiple TSCC cases for use in an integrated immunogenomic study.

Methods: Immune landscape of TSCC was depicted by expression microarray data from GSE13601 and GSE34105. Univariate Cox analysis, in combination with survival analysis, was applied to select candidate IRGs with significant survival value. Survival predicting models were constructed by multivariate Cox regression and logistic regression analysis. Unsupervised clustering analysis was used to construct an immune gene panel based on prognostic IRGs to distinguish TSCC subgroups with different prognostic outcomes. Finally, IHC staining was performed to validate the clinical value of this immune-gene panel.

Results: Differentially expressed IRGs were identified in two TSCC microarray datasets. Functional enrichment analysis revealed that ontology terms associated with variations in T cell function, were highly enriched. Infiltration status of activated CD8+ T cells, central memory CD4+ T cells and type 17 T helper cells, had great prognostic value for TSCC progression. Unsupervised clustering analysis was further performed to classify TSCC patients into three subgroups. CTSG, CXCL13, and VEGFA were finally combined together to form an immune-gene panel, todistinguish different TSCC subgroups. IHC staining of TSCC sections further validated the clinical efficiency of the immune-gene panel consisting of prognostic IRGs to distinguish TSCC patients.

Conclusion: VEGFA, CXCL13, and CTSG, correlated with T cell infiltration and prognostic outcome. They were screened to form an immune-gene panel to identify TSCC subgroups with different prognostic outcomes. Clinical IHC further validated the efficacy of this immune-gene panel to evaluate aggressiveness of TSCC development.

Keywords: CTSG; CXCL13; Oral squamous cell carcinoma; VEGFA; immune-related genes; tongue squamous cell carcinoma.

Figure 1

Immunogenomic landscape of tongue squamous cell carcinoma. A. The immune cell infiltration status in GSE13601 depicted by heatmap. B. Correlation between infiltration of immune cells undertaking anti-tumor functions and those executing pro-tumor functions in GSE13601. C. Immune cell infiltration status in GSE34105 depicted by heatmap. D. Correlation between infiltration of immune cells undertaking anti-tumor functions and those executing pro-tumor functions in GSE34105. E, F. Representative images of activated CD4+ T cell and activated CD8+ T cell infiltration status between normal and tumor tissues in GSE13601 and GSE34105 respectively, depicted by boxplots.

Figure 2

Differentially expressed immune-related genes. A. Differentially upregulated and downregulated genes selected from GSE13601. P-value <0.05 and |log2 Fc| >1 were set as cutoff values for the selection of differentially expressed genes. B. Differentially expressed genes in GSE13601 depicted by heatmap. C. Differentially upregulated and downregulated genes selected from GSE34105. P value <0.05 and |log2 Fc| >1 were set as cutoff values for the selection of differentially expressed genes. D. Differentially expressed genes in GSE34105 depicted by heatmap. E. Venn diagram showing the relationship between differentially expressed genes and immune-related genes. The common part indicates IRGs differentially expressed by tongue squamous cell carcinoma.

Figure 3

Gene functional enrichment analysis of differentially expressed immune-related genes. A. Terms of significantly enriched gene ontology and pathways in GSE13601. B. A network plot depicting interactions among those enriched ontological functions in GSE13601. C. Terms of significantly enriched gene ontology and pathways in GSE34105. D. A network plot depicting interactions among those enriched ontological functions in GSE34105. E. Enriched GO terms by shared immune-related genes between GSE13601 and GSE34105. F. Enriched KEGG terms by shared immune-related genes between GSE13601 and GSE34105.

Figure 4

Protein-protein interaction (PPI) analysis. A. PPI network constructed by differentially expressed immune-related genes in GSE13601. B. PPI network constructed by differentially expressed immune-related genes in GSE34105. C. Core PPI modules selected by MCODE function in GSE13601. D. Core PPI modules selected by MCODE function in GSE34105.

Figure 5

Survival analysis of differentially expressed immune-related genes. A. Hazard ratio univariate Cox regression model constructed by immune-related genes with potential survival value. B-I. Kaplan-Meier curves depicting survival analysis of immune-related genes with significant prognostic value (BPIL2, CTSG, CXCL13, VEGFA, ADIPOQ, DKK1, LEFTY1, NPY).

Figure 6

Survival analysis of cancer-related infiltrating immune cells. A-C. Kaplan-Meier curves depicting survival analysis of three types of infiltrating immune cells (activated CD8+ T cells, central memory CD4+ T cells and type 17 T helper cells) with high prognostic value.

Figure 7

Description of relationship between IRG expression and immune cell infiltration. A-C. Negative correlations between VEGFA expression level and enrichment scores for activated CD8+ T cells, central memory CD4+ T cells and type 17 T helper cells. D-F. Positive correlations between CXCL13 expression level and enrichment scores for activated CD8+ T cells, central memory CD4+ T cells and type 17 T helper cells. G-I. Positive correlations between VEGFA expression level and enrichment scores for activated CD8+ T cells, central memory CD4+ T cells and type 17 T helper cells. The significance of all these correlation analysis depicted above were set by the cutoff P-value at 0.05.

Figure 8

Construction of prognosis-predicting models using selected IRGs with significant prognostic value. A-C. ROC curves depicting the test efficiency of expression levels of VEGFA, CXCL13 and CTSG respectively. The measures of area under curve were defined as AUC value, which indicated their clinical outcome-predicting efficiency precisely. D, E. Information of prognostic-predicting models constructed by logistic regression analysis and multivariate Cox regression analysis. F, G. ROC curves depicting the test efficiency of logistic regression model and multivariate Cox regression model respectively. The measures of area under curve were defined as AUC value, which indicated their clinical outcome-predicting efficiency precisely.

Figure 9

Exploration of an immune-gene panel to distinguish different TSCC clusters. A. Three types of TSCC subgroups, named cluster A, cluster B and cluster C, were constructed by unsupervised cluster analysis based on prognostic IRGs selected above. B-D. Prognostic analysis between different TSCC clusters. E-G. Comparison of activated CD8+ T cells, central memory CD4+ T cells and type 17 T helper cells infiltration status among different TSCC clusters.

Figure 10

Identification of TSCC immune subgroups by immune-gene panel. A-C. Relative expression levels of VEGFA, CXCL13 and CTSG in tongue squamous cell carcinoma samples among different TSCC subgroups in TCGA dataset. D-F. Positive and negative rate of VEGFA, CXCL13 and CTSG among different TSCC subgroups in TCGA dataset. Expression data of TSCC samples from TCGA dataset were homogenized for calculation of positive rate corresponding to each IRG.

Figure 11

Clinical validation for the predicting efficiency of immune-gene panel. A. Relative expression levels of CTSG, CXCL13 and VEGFA in tongue squamous cell carcinoma samples among different T stages in TCGA dataset. B. Relative expression levels of CTSG, CXCL13 and VEGFA in tongue squamous cell carcinoma samples among different N stages in TCGA dataset. C. Clinical validation of VEGFA, CXCL13 and CTSG expression levels in tumor samples at pathologic T1, T3 and T4 stage with the help of IHC staining images. D. Clinical validation of VEGFA, CXCL13, and CTSG expression levels in tumor samples at pathological N0 and N2 stages with the help of IHC staining images. E-G. Clinical validation for accumulation of Th17 cells, CD4+ T cells and CD8+ T cells in TSCC by immunofluorescent staining.

Figure 12

Validation for the role of BPIL2, CXCL13 and CTSG in TSCC development and metastasis. A-C. Role of BPIL2, CXCL13 and CTSG in TSCC malignant development validated by Cell Counting Kit (CCK8) assay. D-F. Role of BPIL2, CXCL13 and CTSG in TSCC metastasis confirmed by transwell assay. Si-RNA technique was implemented for efficient knockdown of BPIL2, CXCR5 and CTSG, and two target sites for each gene were applied in our experiment.