Identification of a cytokine-dominated immunosuppressive class in squamous cell lung carcinoma with implications for immunotherapy resistance

Abstract

Background: Immune checkpoint blockade (ICB) therapy has revolutionized the treatment of lung squamous cell carcinoma (LUSC). However, a significant proportion of patients with high tumour PD-L1 expression remain resistant to immune checkpoint inhibitors. To understand the underlying resistance mechanisms, characterization of the immunosuppressive tumour microenvironment and identification of biomarkers to predict resistance in patients are urgently needed.

Methods: Our study retrospectively analysed RNA sequencing data of 624 LUSC samples. We analysed gene expression patterns from tumour microenvironment by unsupervised clustering. We correlated the expression patterns with a set of T cell exhaustion signatures, immunosuppressive cells, clinical characteristics, and immunotherapeutic responses. Internal and external testing datasets were used to validate the presence of exhausted immune status.

Results: Approximately 28 to 36% of LUSC patients were found to exhibit significant enrichments of T cell exhaustion signatures, high fraction of immunosuppressive cells (M2 macrophage and CD4 Treg), co-upregulation of 9 inhibitory checkpoints (CTLA4, PDCD1, LAG3, BTLA, TIGIT, HAVCR2, IDO1, SIGLEC7, and VISTA), and enhanced expression of anti-inflammatory cytokines (e.g. TGFβ and CCL18). We defined this immunosuppressive group of patients as exhausted immune class (EIC). Although EIC showed a high density of tumour-infiltrating lymphocytes, these were associated with poor prognosis. EIC had relatively elevated PD-L1 expression, but showed potential resistance to ICB therapy. The signature of 167 genes for EIC prediction was significantly enriched in melanoma patients with ICB therapy resistance. EIC was characterized by a lower chromosomal alteration burden and a unique methylation pattern. We developed a web application ( http://lilab2.sysu.edu.cn/tex & http://liwzlab.cn/tex ) for researchers to further investigate potential association of ICB resistance based on our multi-omics analysis data.

Conclusions: We introduced a novel LUSC immunosuppressive class which expressed high PD-L1 but showed potential resistance to ICB therapy. This comprehensive characterization of immunosuppressive tumour microenvironment in LUSC provided new insights for further exploration of resistance mechanisms and optimization of immunotherapy strategies.

Keywords: Immune checkpoint blockade resistance; Immunogenomics; Immunosuppressive cytokine; LUSC; T cell exhaustion; Tumour microenvironment.

Fig. 1

The identification and molecular characterization of EIC. A The heatmap of gene expression clusters for 250 late-stage (IIA-IV) LUSC samples by unsupervised NMF illustrates 4 distinct expression patterns. B Stromal and immune enrichment analysis defined the cluster 2 of four expression patterns as an immune-stromal cluster. High and low gene enrichment scores are delineated in red and grey, respectively. C The enrichment scores of gene signatures identified the immune cells for the immune-stromal and other clusters. D The comparison of the absolute fractions of TME cells inferred by CIBERSORT between two classes. E,F Box plots show the differences of leukocyte fraction and TIL percentage between two classes. G Box plots show different expression levels of multiple inhibitory receptors in the immune-stromal cluster compared to the other clusters. H The consensus-clustered heatmap of 250 LUSC samples defined the immune-stromal cluster as exhausted immune class (EIC). High and low gene enrichment scores are represented in red and grey, respectively. I GSEA analysis indicated the EIC showed significant enrichments of hallmark gene sets and KEGG pathways related to cytokine, T cell receptor, epithelial mesenchymal transition, and apoptosis. J The functionally grouped network of KEGG pathways by ClueGO/CluePedia for the interpretation of metagene-specific genes’ biological roles. Colourless and colour nodes represent metagene-specific genes and KEGG pathway terms, respectively. Node colours represent distinct functional groups. Node size represents the significance of KEGG pathways. The more significant KEGG pathways are, the larger highlighted nodes. All statistical differences of two classes were compared by Wilcoxon rank-sum test; *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001

Fig. 2

Internal validation of EIC on 247 early-stage (I–II) LUSC samples. A The heatmap of gene expression clusters for 247 early-stage (IIA–IV) LUSC samples by unsupervised NMF illustrates 4 distinct expression patterns. B Heatmap shows the cluster 2 (defined as EIC) exhibited high enrichment scores of gene signatures of T cell exhaustion, immunosuppressive cells, and immunosuppressive cytokine. C The comparison of the absolute fraction of TME cells between the EIC and the rest class. D,E Box plots show the differences of leukocyte fraction and TIL percentage between two classes. F Box plots shows the different expression levels of multiple inhibitory receptors between two classes. G Cytokine-, T cell receptor-, and epithelial mesenchymal transition-related hallmark gene sets and KEGG pathways enriched in the EIC. H ROC curve evaluated the predictive capacity of 167 exhausted immune classifier genes in early-stage LUSC samples. All statistical differences of two groups were computed by Wilcoxon rank-sum test; *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001

Fig. 3

Prognosis analysis for the EIC and the rest class across different stages of LUSC. A–C Kaplan–Meier estimates of overall survival for the EIC and the rest class across late-stage, early-stage, and all-stage LUSC. Survival data was limited to maximum 5 years (60 months). D–F Kaplan–Meier estimates of progression-free survival for the EIC and the rest class across late-stage, early-stage, and all-stage LUSC. P-values were calculated by log-rank test. Survival data was limited to maximum 5 years (60 months). G,H Multivariate Cox regression analysis on four variables (class, gender, tumour stage, and age) for late-stage and all-stage LUSC

Fig. 4

Prediction of ICB therapy resistance. A,B Patients in the EIC showed a relatively higher expression level of PD-L1 and higher TIDE prediction score for ICB therapy. C Metastatic melanoma patients with no response to anti-PD-1 therapy had higher enrichment scores of 167 exhausted immune classifier genes compared to patients with response. D The EIC showed higher enrichment score than the rest class in all tumour stage TCGA LUSC. E Box plot shows higher expression of TGFB1 in the EIC than rest class

Fig. 5

Distinctive methylation signatures characterized the EIC of LUSC. A Hierarchical clustering heatmap of 216 CpG methylation values located within 162 immunosuppression-related gene promoters show significant difference between the EIC and the rest class (FDR< 0.05, diff > 0.2). B Boxplot displays the mean methylation levels of 216 CpG sites within 162 immune exhaustion-related gene promoters for 2 classes. Wilcoxon rank-sum test (p<0.0001). Exhausted vs Rest: p=5.1E−15. C–E Correlations between expression and promoter methylation levels for deregulated genes in the EIC to the rest class. CARTN exhibited significantly lower expression in the EIC, while SMAD7, IRF7, CCR4, and MYO1G were synergistically overexpressed (P < 0.001, Wilcoxon rank-sum test). D The expressions of ARTN, SMAD7, IRF7, CCR4, and MYO1G negatively correlated with their promoter methylation level for the whole training cohorts. ESMAD7, IRF7, CCR4, and MYO1G had lower promoter methylation levels mirroring higher expression levels in the EIC, whereas ARTN had an opposite status. Red dots in the plotting represent the members of EIC, and blue dots represent the members of the rest class

Fig. 6

Association of EIC with somatic mutations, neoantigens, and copy number alteration. A, D The landscape of most frequently mutated genes between the EIC and the rest class in late-stage and early-stage LUSC, respectively. B, E Box plots show the number of mutations between two classes in late-stage and early-stage LUSC, respectively. C, F Box plots show the number of neoantigens between two classes in late-stage and early-stage LUSC, respectively. G Pearson correlation analysis between scaled the EIC score and the number of mutations across LUSC cohorts of late-stage, early-stage, and all-stage, respectively. H, L Box plots show significant difference of amplification burden of cytoband between two classes in late-stage and early-stage LUSC, respectively. I, M Box plots show significant difference of deletion burden of cytoband between two classes in late-stage and early-stage LUSC, respectively. For late-stage (J) and early-stage (N) LUSC, the frequency of patients with amplification of driver genes in two classes. For late-stage (K) and early-stage (O) LUSC, the frequency of patients with deletion of driver genes in two classes. All statistical significances of two classes were computed by Wilcoxon rank-sum test; *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001

Fig. 7

IDO protein expression analysis in EIC. A, B Boxplots showing protein expression difference between the EIC and the rest class in late-stage and early-stage LUSC, respectively. Wilcoxon rank-sum test was utilized to perform comparision analysis. C The immunohistochemically stained tissue images show different expression levels of IDO protein across three LUSC patients from HPA