Granzymes expression patterns predict immunotherapy response and identify the heterogeneity of CD8 + T cell subsets

Abstract

BACKGROUND:

Recent studies illustrated the effects of granzymes (GZMs) gene alterations on immunotherapy response of cancer patients. Thus, we aimed to systematically analyze the expression and prognostic value of GZMs for immunotherapy in different cancers, and identified heterogeneity of the GZMs expression-based CD8 ${}^{+}$ T cell subsets.

METHODS:

First, we analyzed GZMs expression and prognostic value at pan-cancer level. Meanwhile, we established a GZMs score by using the single-sample gene set enrichment analysis (ssGSEA) algorithm to calculate the enrichment scores (ES) based on a gene set of five GZMs. The potential value of GZMs score for predicting survival and immunotherapy response was evaluated using the tumor immune dysfunction and exclusion (TIDE) and immunophenoscore (IPS) algorithm, and we validated it in immunotherapy cohorts. CellChat, scMetabolism, and SCENIC R packages were used for intercellular communication networks, quantifying metabolism activity, and regulatory network reconstruction, respectively.

RESULTS:

The GZMs score was significantly associated with IPS, TIDE score. Patients with high GZMs score tended to have higher objective response rates of immunotherapy in melanoma and urothelial carcinoma. GZMs expression-based CD8 ${}^{+}$ T cell subsets presented heterogeneity in functions, metabolism, intercellular communications, and the tissue-resident memory programs in lung adenocarcinoma (LUAD). The transcription factors RUNX3 and ETS1, which may regulate the expression of GZMs, was found to be positively correlated with the tissue-resident memory T cells-related marker genes.

CONCLUSIONS:

The higher GZMs score may indicate better response and overall survival (OS) outcome for immunotherapy in melanoma and urothelial carcinoma but worse OS in renal cell carcinoma (RCC). The GZMs score is a potential prognostic biomarker of diverse cancers. RUNX3 and ETS1 may be the potential targets to regulate the infiltration of GZMs expression-based CD8 ${}^{+}$ T cell subsets and affect the tissue-resident memory programs in LUAD, which may affect the prognosis of LUAD patients and the response to immunotherapy.

1. Introduction

Cancer has ranked as a leading cause of death worldwide, and immunotherapy rapidly emerged with effectiveness in treating various cancers by influencing the interaction between the human immune system and cancer [1, 2]. Despite immunotherapy being widely applied in several cancers, only a minority of patients experience durable survival [3]. There is a general recognition of the urgent need for excavating predictive biomarkers to assess the response to immunotherapy.

Human Granzymes (GZMs) comprise a family of five (granzyme A/GZMA, granzyme B/GZMB, granzyme H/GZMH, granzyme K/GZMK, granzyme M/GZMM) secreted serine proteases, which were mainly found in the cytotoxic granules of cytotoxic T lymphocytes (CTLs) and natural killer (NK) cells. CTLs and NK cells used the granule exocytosis pathway to target cell death by delivering the contents of cytotoxic granules to the surface of target cells [4]. The perforin, one of the granule proteins, was thought to facilitate GZMs entry into target cells by physically forming holes in the cell membrane [5, 6]. Previous studies of GZMs mainly focused on the two most abundant members (GZMA and GZMB) and their role in killing cancer cells [7]. GZMB was believed to induce target cell apoptosis by cleaving caspase-3 or caspase-3 substrates [7, 8, 9]. Recently, Zhou et al. [10] reported that GZMA from cytotoxic lymphocytes cleaves gasdermin B (GSDMB) to trigger pyroptosis in GSDMB-positive cells. Kontani et al. [11] reported that GZMB might play a role in suppressing nodal metastasis of cancer cells in breast and lung cancers.

The high expression of the GZMs in cytotoxic cells and the well-established role of these cells in cancer immune surveillance lead to a general belief that GZMs are mainly protective cytotoxic proteases that block cancer development. However, GZMs may also be crucial regulators of biological processes involved in cancer development and progression, such as immune homeostasis, inflammation, angiogenesis, and invasiveness [12]. Findings accumulated over the past several years on functions of GZMs other than cell death. Metkar et al. [13] have shown that extracellular granzyme A can promote the release of pro-inflammatory cytokines from human monocytes and murine peritoneal macrophages. The identification of GZMs expression in cells other than CTLs and NK cells, for example, certain GZMs are expressed and secreted by B cells, mast cells, keratinocytes, basophils, macrophages, dendritic cells, and regulatory T-cells, in the absence of detectable perforin, indicating that GZMs may have hitherto unsuspected roles in immunity [14, 15, 16, 17, 18, 19, 20, 21].

Furthermore, a comprehensive understanding of GZMs in pan-cancer is still lacking. The effects of GZMs have not been illustrated clearly. The heterogeneity of different cancer types may cause the different cellular components of tumor tissue and make GZMs perform different functions. We integrated the expression of GZMs as a score and explored its associations with immune infiltration, biological processes, immunotherapy response, and prognostic value in pan-cancer. Furthermore, we identified the subsets of CD8 ${}^{+}$ T cells based on the GZMs expression and explored their heterogeneity.

2. Materials and methods

2.1 Data sources

This study was performed mainly based on the Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. The TCGA data, including Transcripts Per Kilobase Million (TPM) expression and clinical information (survival status, survival time, grades, stages, etc.) and Copy number variation (CNV) data, were downloaded from UCSC XENA (http://xenabrowser.net/). Furthermore, eight immunotherapy-related datasets were screened to conduct the survival analysis: Three (Gide et al. 2019 [23], Nathanson et al. 2017 [24], Van Allen et al. 2015 [25]) were downloaded from Tumor Immune Dysfunction and Exclusion website (http://tide.dfci.harvard.edu/), three (GSE100797 [26], GSE91061 [27], GSE78220 [28]) were downloaded from the GEO database, two (IMvigor210CoreBiologies [29], Braun et al. 2020 [30]) were obtained from two research articles.

2.2 GZMs score

The ssGSEA algorithm [22] was used to calculate the ES. The ES calculated based on a gene set of five granzymes by ssGSEA was defined as GZMs score. The GZMs scores represented the integrated level of the expression of five GZM genes in each sample.

2.3 Gene expression and clinical analysis

We used packages edgeR [23] for differential gene expression analysis and calculated the log ${}_2$ (Fold Change) to compare the tumor to adjacent normal tissue groups. The relations between the clinical features and GZMs gene expression/GZMs scores were conducted by the Kruskal-Wallis test. The age was divided into two groups (age $>$ 65 and $⩽$ 65 years old). All of the other clinical features were categorical variables. We also explored the expression of GZMs based on the classification developed by Thorsson et al. [24], which was a novel global immune classification based on the transcriptomic profiles from all 33 TCGA cancer types.

2.4 Survival analysis

Before Kaplan-Meier curve analysis, the continuous variables were divided into binary variables by the optimal values calculated by package survminer. Log-rank tests were used to test the different survival rates between the two groups. For immunotherapy datasets from the same cancer type, the meta-analysis assay was used to calculate the hazard ratio (HR) and describe the pooled effects of the GZMs score as a continuous variable.

2.5 Gene set enrichment analysis

We calculated the correlation coefficients between the GZMs scores and every gene of the whole transcript, respectively, in each cancer type. The gene list was then re-ranked by the correlation coefficients and inputted into gene set enrichment analysis (GSEA) based on the package ClusterProfiler [25]. The signaling pathways information in the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/) were used to carry out the GSEA analysis, including the gene sets from the Kyoto Encyclopedia of Genes and Genomes (KEGG), gene ontology (GO) terms, and Hallmarks. The universally up-regulated (NES $>$ 2, $q<$ 0.05, and more than 17 cancer types met this standard) and down-regulated (NES $<-$ 1, $q<$ 0.05, and more than 17 cancer types met this standard) pathways associated with GZMs scores were partly visualized.

2.6 Estimation of immune-cell type fractions

Cell-type identification is a deconvolution algorithm that calculates the putative proportion of immune cell fraction based on gene expression profiles by estimating relative subsets of RNA transcripts (CIBERSORT algorithm) [26]. We used the reference set with 22 sorted immune cell subtypes to analyze the abundance of tumor-infiltrating immune cells in TCGA data. Normalized gene expression data were analyzed using the CIBERSORT algorithm. In the scRNA-seq analysis, we used Bayesian cell proportion reconstruction inferred using statistical marginalization (BayesPrism) [27], to predict the cellular composition and gene expression in individual cell types from bulk RNA-seq, using patient-derived, scRNA-seq as prior information.

2.7 Methylation profiles and mutation profiles

GSCALite platform (http://bioinfo.life.hust.edu.cn/web/GSCALite/) was used to analyze differential methylation of GZMs in tumors and normal tissues in pan-cancer [28]. Subsequently, we explored the correlation between GZMs gene methylation and its expression in different cancers.

The copy number alteration (CNA) and mutation landscape of GZMs in pan-cancer were gained from cBioPortal (http://www.cbioportal.org) [29]. Furthermore, GSCALite was used to investigate the SNP frequencies and heterozygous/homozygous CNA of GZM genes, as well as the correlation between CNA and mRNA expression of GZMs. Tumor Mutational Burden (TMB) was calculated based on all somatic mutations and used the 35 megabases (Mb) as the total number of bases in human beings.

2.8 Estimation of response to immunotherapy and drug sensitivity

The predictive response to immune checkpoint blockade (ICB) therapy was estimated by using the TIDE and IPS algorithms. Peng Jiang et al. [30] developed TIDE, a computational method to model two primary mechanisms of tumor immune evasion: the induction of T cell dysfunction in tumors with high infiltration of CTL and the prevention of T cell infiltration in tumors with low CTL levels. Then the correlation between the expression of GZMs and the predictive immune response was explored. The IPS could be calculated by analyzing gene expression profiles in the tumor samples comprising four major categories of genes that determine immunogenicity effector cells, immunosuppressive cells, MHC molecules, and immunomodulators [31]. Then the IPS was calculated with a range of 0–10 based on the z-score for gene expression of representative cell types. A prophetic algorithm [32] was used to establish a ridge regression model to predict the potential sensitivity of 33 drugs based on the expression profile in the TCGA database.

2.9 Analysis of single-cell RNA sequencing data

One scRNA-seq dataset was analyzed to investigate the expression of the GZMs and the potential intercellular communication in single-cell resolution. Unsupervised clustering was based on the gene expression profiles using the R package Seurat [33] and passed to t-Distributed Stochastic Neighbor Embedding (t-SNE) or Uniform Manifold Approximation and Projection (UMPA) for clustering visualization. The clusters were annotated manually or referenced to the inferred results from the package SCINA [34], an automatic cell type detection and assignment algorithm for scRNA-seq. The reference data were downloaded from the supplementary data of a study developed by J Javier Diaz-Mejia [35]. The malignant cells were recognized by the package infercnv [36].

CD8 ${}^{+}$ T cells were divided into five clusters (CD8 ${}^{+}$ Nave-like T cells, CD8 ${}^{+}$ Cytotoxic T cells, CD8 ${}^{+}$ Exhausted T cells, CD8 ${}^{+}$ GZMK high T cells, CD8 ${}^{+}$ Proliferating T cells) according to the expression of signatures genes and their expression of GZMs manually.

We conducted the enrichment analysis based on the differential genes between different CD8 ${}^{+}$ T cell subsets. The intercellular communications were analyzed by using the package CellChat [37]. CellChat was a tool that quantitively inferred and analyzed intercellular communication networks from scRNA-seq data. CytoTRACE [38] (Cellular Trajectory Reconstruction Analysis using Gene Counts and Expression) and Monocle3 [39] were used to infer the cell differentiation trajectories. Package SCENIC [40] was used for simultaneous gene regulatory network reconstruction and cell-state identification. Package scMetabolism [41] was used to quantify the metabolism activity of CD8 ${}^{+}$ T cell subsets at the single-cell resolution.

2.10 Statistical analysis

All statistical analyses were conducted by R version 4.1.0 (R package: pheatmap, ggplot2, rms, glmnet, forest plot, limma, GSVA, survminer, survival ROC, beeswarm, clusterProfile, Seurate, SingleR, etc.). Unless otherwise stated, a two-tailed $p<$ 0.05 was regarded as statistical significance for all analyses. The Kruskal-Wallis test examined the statistical difference of distribution in three or more groups and the Wilcoxon test compared that of two groups in most parts of this study. The correlation analysis in our study was conducted based on the Pearson method unless otherwise stated.

3. Results

3.1 The methylation analysis of GZMs in human cancers

DNA methylation is related to changes in gene expression in cancers. Thus, the DNA methylation differences of GZMs in various cancers between tumors and adjacent normal tissues were evaluated by the GSCALite platform. The results indicated that all five GZMs methylation was significantly down-regulated in liver hepatocellular carcinoma (LIHC), Kidney Clear Cell Carcinoma (KIRC), Lung Squamous Cell Carcinoma (LUSC), breast carcinoma (BRCA), Head and Neck Squamous Cell Carcinoma (HNSC) (Fig. 1A). Next, we evaluated the correlation between the DNA methylation and mRNA expression of GZMs in pan-cancer. The results exhibited that the GZMs mRNA expression was mainly negatively correlated with its DNA methylation in most tumors, especially in Kidney Papillary Cell Carcinoma (KIRP), Thyroid Carcinoma (THCA), KIRC, LUAD, BRCA, and Testicular Germ Cell Tumors (TGCT) (Fig. 1B).

Figure 1.

The methylation analysis and mutation landscape of GZMs in human cancers. (A) DNA methylation differences of GZMs between tumors and normal tissues in various cancers with GSCALite platform. (B) Correlation between the DNA methylation and mRNA expression of GZMs in 33 cancer types. (C) Landscape of genomic aberrations in the GZMs genes in cancer. Each column represents a patient. The frequency of alterations in five genes of GZMs are presented. Only patients with genomic alterations in the indicated genes are shown. Alteration rates per gene are displayed in the left labels. (D) SNV percentage of GZMs in different cancers. (E) GZMs alteration (mutation and CNAs) frequency in 32 TCGA cancer types. (F) Single Nuclear Tide Polymorphisms of GZMs in pan-cancer. The right panel presents the counts of mutation of each gene respectively. (G) CNV profiles of GZMs in 33 cancer types with the GSCALite platform. Hete Amp, heterozygous amplification; Hete Del, heterozygous deletion; Homo Amp, homozygous amplification; Homo Del, homozygous deletion. (H) Correlation between the CNV and mRNA expression of GZMs in different cancers.

Figure 2.

Expression and clinical analysis of GZMs in pan-cancer. (A) Bubble map of the differential expression of GZMs between tumor and normal samples. Blue and red dots represented down-regulated and up-regulated expression in tumors, respectively, the darker the color, the greater the difference. The size of the dots represented statistical significance. The black outline border indicates the FDR $<$ 0.05. (B) Bubble map of the survival analysis. Hazard ratios and cox $p$ -value were presented through bubble color and size. The bubble color from blue to red represented the hazard ratio from high to low, and bubble size was positively correlated with the –log ${}_{10}$ ( $p$ -value). The black outline border indicates Cox $p$ -value $<$ 0.05. (C) Relationship between expression of GZMs and clinical features. Color-blocks indicated that clinical parameters were associated with expression ( $p<$ 0.05), and white blocks were irrelevant. (D) The violin plot showed the expression of GZMs in five immune classifications. (E) Heatmap of correlation coefficient among GZMs and GZMs score. The numbers were the correlation coefficients. The darker the color, the greater the correlation coefficient. (F) Representative IHC images of GZMs expression in breast cancer, renal cancer, renal cancer tissue, lymphoma, normal breast tissue, normal kidney tissue, normal kidney tissue, and normal lymph node tissue.

Figure 3.

Systematically analysis of GZMs scores in pan-cancer. (A) Box plot presents the differential GZMs scores between the tumor and normal tissues. (B) The forest plot showed survival analysis of the high GZMs score compared with the low GZMs score. The heat map of clinical features showed the correlations between GZMs score and clinical features. Color-blocks indicated that clinical features are associated with GZMs score (Kruskal-Wallis test, $p<$ 0.05), and white blocks are irrelevant. The heat map in the right panel presented the correlation coefficient between GZMs score and cell infiltration or IPS/TIDE values (the dots indicated $p<$ 0.05).

3.2 The genetic alteration of GZMs in human cancers

Analysis of GZMs genetic alterations in pan-cancers was explored with the cBioPortal website and GSCALite platform. The results showed that genetic alterations of GZMs occurred in 1.2%–1.8% of patients (Fig. 1C). It seemed that GZMA and GZMK have the same state of CNV, and so do the GZMB and GZMH (Fig. 1C). Figure 1D provided the percentage of single nucleotide variants (SNV) of GZMs in the pan-cancers, indicating that SNV of GZMs mainly appeared in skin cutaneous melanoma (SKCM) and Uterine Corpus Endometrial Carcinoma (UCEC). The genetic alteration profiling of GZMs showed that its mutation was one of the most important factors for alteration in SKCM and UCEC, while the amplification was the main factor in Sarcoma (SARC) and Adrenocortical Carcinoma (ACC), so do the deep deletion in ovarian serous cystadenocarcinoma (OV) (Fig. 1E). The oncoplot showed that the GZMs mutations were mainly consistent with missense mutations (Fig. 1F).

GSCALite platform provided heterozygous and homozygous CNV profiles in pan-cancers. The pie chart exhibited that the heterozygous CNVs of GZMs were more frequent in various cancers (Fig. 1G). The frequency of amplification and deletion was differing from cancers. Furthermore, the correlation between CNV and gene expression was investigated (Fig. 1H). The result showed that GZMA and GZMK expression was positively correlated with CNV in Stomach Adenocarcinoma (STAD), LUSC, and HNSC. The GZMB and GZMH expression were negatively correlated with CNV in Prostate Adenocarcinoma (PRAD) and THCA.

3.3 Expression and clinical analysis of GZMs in pan-cancer

We performed differential gene expression analysis of GZMs in different types of cancers (Fig. 2A). Differential expression between tumor tissues and adjacent normal tissues was identified (FDR $<$ 0.05). In KIRC, KIRP, and Glioblastoma (GBM), the tumor tissue had higher expression of GZMs compared to normal tissue. The correlations between the expression level of each granzyme and OS by using univariate Cox regression analysis in different cancer types were exhibited (Fig. 2B). Across most cancer types, there was no significantly converse prognosis among different granzymes. However, the GZMs expression indicated different prognosis in cancers.

A comprehensive expression-clinical features analysis of GZMs in 33 types of TCGA cancer (Fig. 2C) was carried out. Except for immune classification, it did not show significant relevance between clinical features and the expression of each granzyme in most cancer types. The subtype C2 (IFN-gamma Dominant), C3 (Inflammatory), and C6 (TGF-beta Dominant) had a higher expression of GZMs, while the C5 (Immunologically Quiet) had a lower expression (Fig. 2D), indicating GZMs expression might be associated with the immune microenvironment. Besides, the expression of GZMs at the protein level also shows the differential expression between the tumor and normal tissue and between different cancers (Fig. 2F).

3.4 Systematically analysis of GZMs score in pan-cancer

We calculated the GZMs score for each sample in pan-cancer as described in the methods. The heatmap of correlation coefficients among GZMs expression and GZMs score presented a significantly strong positive correlation (Fig. 2E). GZMs scores in different types of cancers also showed distinction (Fig. 3A). We explored the prognostic prediction value of the GZMs score and evaluated the associations between the GZMs score and clinical features, immune classifications, estimated components of immune cells, and the propensity of response to immunotherapy. High GZMs scores were associated with significantly better OS in BRCA, head and Neck squamous cell carcinoma (HNSC), LIHC, OV, and SKCM, but worse OS in brain lower grade glioma (LGG) and uveal melanoma (UVM) (Fig. 3B).

Furthermore, we explored the correlations between immune-related features and GZMs scores. We applied the CIBERSORT algorithm to assess the components of the tumor-infiltrating immune cells. In most cancers, each GZM gene expression (Fig. S1A) and GZMs score (Fig. 3B) positively correlated with CD8 ${}^{+}$ T cells and M1 Macrophages. The strong correlation between the GZMs score and IPS, TIDE, TMB, and PDCD1 expression also indicated that the GZMs score and GZMs gene expression might have a good predictive value for immunotherapy. The higher GZMs score generally be associated with higher IPS and lower TIDE scores (Fig. 3B, Fig. S2A and B). In LUAD, LUSC, and Colon Adenocarcinoma (COAD), the GZMs score was positively associated with TMB (Fig. S2C). Furthermore, the higher GZMs score was associated with higher expression of PDCD1 in most cancers (Fig. S2D), including LUAD, LUSC, SKCM, Bladder Urothelial Carcinoma (BLCA), etc. These findings also indicated that the GZMs score might have a good predictive value for immunotherapy and had the potential to be applied in more cancers.

To determine patients’ responses to anti-cancer drugs, we also inferred the predictive drug sensitivity of the 33 drugs in TCGA data. High GZMs scores were usually associated with lower predictive value in most cancer types, which meant that tumor tissues with high levels of GZMs might be more sensitive to certain drugs (Fig. S2E). In LUAD, the results presented that patients with higher GZMs scores may be more sensitive to Cisplatin, Docetaxel, Etoposide, Paclitaxel, and Erlotinib (Fig. S2F–K).

3.5 Gene set enrichment analysis based on GZMs score

Figure 4.

Gene set enrichment analysis based on the rank of the correlation coefficients. (A) Heatmap of pathways that were positively associated with GZMs score. The colors represented the normalized enrich score (NES) of this pathway in certain cancer types. (B) Heatmap of pathways that were negatively associated with GZMs score. The colors represented the NES of this pathway in certain cancer types. (C) GSEA on certain gene sets in BRCA and UVM, including 4 positively associated pathways and 2 negatively associated pathways. (D) The bar plot showed the pathways behaved inversely in BRCA and UVM. The NES was plotted on the X-axis.

Figure 5.

Survival analysis of GZMs score in immunotherapy datasets. (A) Kaplan-Meier curves estimated the OS between the high and the low GZMs score group (divided by the cut-off value based on the lowest $p$ -value) in eight immunotherapy datasets. (B) Forest plot visualized the hazard ratios (95% confidence intervals) of OS of high and low GZMs scores in melanoma. (C) The relative fraction of PR/CR or SD/PD of immunotherapy in high and low GZMs score group.

Figure 6.

Characters of GZMs expression in the single-cell resolution. (A) scRNA-Seq data from 15 LUAD samples were projected onto the t-SNE and then annotated (left part), the immune cells and the CD8 ${}^{+}$ T cells were identified and projected onto the t-SNE (top right) and UMAP (bottom right), respectively. (B) Bubble map of the differential expression of GZMs among different types of immune cells and CD8 ${}^{+}$ cell subsets, respectively. The colors represented the average expression of each granzyme in different cell types. The size of the dots represented the proportion of cells expressing corresponding granzymes. (C) CytoTRACE was used to predict the trajectory of the differentiation state of CD8 ${}^{+}$ cells. The predicted differentiation state was visualized in UMAP plots based on the differentiation scores with the dotted arrow (left part). The color of the dots presented the pseudo-time. Pseudo-time analysis of CD8 ${}^{+}$ T cells revealed the trajectory of CD8 ${}^{+}$ T cells from the naïve-like to exhausted CD8 ${}^{+}$ T cells (red arrow) and the trajectory from proliferating CD8 ${}^{+}$ T cells to exhausted cells (blue arrow). The kinetics plot showed the relative expression of GZMB/K across developmental pseudo-time (right part). The dots were colored according to the cell types. (D) Summary of ligand-receptor interactions from GZMB-high subset, GZMK-high subset, and Naïve-like subset to other clusters. $P$ -values (permutation test) are represented by the size of each circle. The color gradient indicates the level of interaction.

Figure 7.

Heterogeneities of CD8 ${}^{+}$ T cell subsets. (A) The UMAP representation of CD8 ${}^{+}$ T cell subsets, including the Naïve-like subset, GZMB-high subset, and the GZMK-high subset. Subsequent analyzes were based on the three subsets. (B)–(C) Kaplan-Meier curves estimated the OS between the high and the low fraction of the GZMK/GZMB-high subsets group, respectively (divided by the cut-off value based on the lowest $p$ -value). (D)–(F) Volcano plots of differential genes in three different CD8 ${}^{+}$ T cell subsets. Red plots represented FDR $<$ 0.05 and log2 (Fold Change) $>$ 0.75, while blue plots represented FDR $<$ 0.05 and log2FC $<-$ 0.75. (G)–(I) Enrichment analysis based on GSEA. The pathways enriched in different subsets were labeled.

We performed the GSEA analysis to explore the pathways associated with the GZMs score (Fig. 4A), positively related pathways were mainly associated with cytotoxicity (positive regulation of leukocyte mediated cytotoxicity, natural killer cell-mediated cytotoxicity, etc.), immune response (cellular defense response, interferon-alpha response, interferon-gamma response, Inflammatory response), activation, and the chemotaxis of immune cells (monocyte chemotaxis, lymphocyte chemotaxis, chemokine signaling pathway). The high GZMs score was associated with the activation of immune-related pathways involved in anti-tumor in most cancers. It was noted that these positively related pathways generally had lower Normalized Enrichment Scores (NES) in thymoma (THYM), UVM, and acute myeloid leukemia (LAML). In these cancers, the high GZMs score was associated with worse survival. As for universally negatively related pathways (Fig. 4B), the higher GZMs score was associated with the inhibition of the AMPK signaling pathway, mTOR signaling pathway, Wnt signaling pathway, and Hippo signaling pathway. The UVM showed apparent differences compared to other cancer types (for example, BRCA), so we further identified differential pathways between UVM and BRCA (Fig. 4C and D). Differential pathways mainly focused on the vesicle-related pathways and co-translational proteins targeting to membrane, which might be related to prognosis.

3.6 Validation of GZMs score for predicting immunotherapy response

We validated the prediction ability of GZMs score for immunotherapy response in eight immunotherapy datasets. High GZMs scores group was associated with better OS in urothelial carcinoma (Invigor210, HR $=$ 0.56, 95% CI [0.38, 0.82], $p=$ 0.003) and melanoma (Gide, HR $=$ 0.35, 95% CI [0.16, 0.75], $p=$ 0.005; GSE91061, HR $=$ 0.27, 95% CI [0.12, 0.64], $p=$ 0.002; Nathanson, HR $=$ 0.12, 95% CI [0.02, 0.91], $p=$ 0.015, Van, HR $=$ 0.32, 95% CI [0.14, 0.71], $p=$ 0.002) (Fig. 5A). In all melanoma datasets, the HR values (HR $=$ 0.33, 95% CI [0.19–0.58]) pooled by meta-analysis demonstrated the significant correlation between GZMs score and OS in the context of immunotherapy (Fig. 5B). In renal cell carcinoma dataset (RCC, Braun2020), the high GZMs score group had a poor OS (HR $=$ 1.36, 95% CI [1.04–1.78], $p=$ 0.025). Except for Braun 2020 dataset and the GSE78220 dataset, patients with high GZMs scores presented higher objective response rates (Fig. 5C).

3.7 Characters of GZMs expression in the single-cell resolution

In this part, we analyzed characters of GZMs gene expressions in a scRNA-seq dataset (GSE131907), which included deep scRNA-seq from 15 samples of lung adenocarcinoma (LUAD). 11 early-stage and 4 advanced-stage tumor samples were included in our analyzes. After quality control, we applied principal component analysis on variably expressed genes across all cells to explore the cellular composition. In the first step, we divided these cells into 31 clusters and annotated these clusters manually (Fig. S3A and C). The malignant cells were identified by the inferCNV method (Fig. S3B). In the second step, we clustered and annotated immune cells by the SCINA method (Fig. 6A). Then we clustered the CD8 ${}^{+}$ T cells into five subsets (CD8 ${}^{+}$ Naïve-like T cells, CD8 ${}^{+}$ Cytotoxic T cells, CD8 ${}^{+}$ Exhausted T cells, CD8 ${}^{+}$ GZMK high T cells, CD8 ${}^{+}$ proliferating T cells) according to signature genes expressions (Fig. S3C). The expression of GZMs in immune cells and CD8 ${}^{+}$ T cell subsets were shown in Fig. 6B. All five GZMs were mainly expressed in CD8 ${}^{+}$ T cells, NK cells, and Gamma Delta T cells. CD8 ${}^{+}$ Naïve-like T cells expressed low-level GZMs, CD8 ${}^{+}$ GZMK-high T cells expressed high-level GZMK and low-level GZMB, while CD8 ${}^{+}$ Cytotoxic T cells and CD8 ${}^{+}$ Exhausted T cells were inverse.

3.8 Pseudo-time analysis of CD8 ${}^{+}$ T cells

To investigate the developmental progression of CD8 ${}^{+}$ T cells, we ordered these cells along the differentiation stages by carrying out the trajectory analysis using Monocle 3 and CytoTRACE packages. CytoTRACE analysis presented the differentiation status of CD8 ${}^{+}$ T cells (Fig. S3D). The estimated pseudo-time from Monocle 3 by assigning the CD8 ${}^{+}$ Naïve-like T cells clusters as the root nodes elucidated the single-cell trajectory path in the order of clusters Naïve-like, Cytotoxic, GZMK-high, and Exhausted. It is noted that the CD8 ${}^{+}$ proliferating T cells seemed to involve a transition to CD8 ${}^{+}$ Exhausted T cells (Fig. 6C, Fig. S3D). With the development of the pseudo-time constructed by Monocle3, the GZMB increased gradually and maintained a relatively high level even in the exhausted state. The GZMK decreased sharply when the state transited from the GZMK-high state to the Exhausted state. The expression of these two granzymes presented different tendencies (Fig. 6C). The expression of GZMA, GZMM, and GZMH didn’t show apparent changes (Fig. S3E–G), so our subsequent analysis focused on GZMB and GZMK.

3.9 Heterogeneity analysis of GZMs-based CD8 ${}^{+}$ T cell subsets

According to the distribution of GZMs in CD8 ${}^{+}$ T cells, we set the CD8 ${}^{+}$ Cytotoxic, the CD8 ${}^{+}$ Exhausted T cells, and the CD8 ${}^{+}$ proliferating T cells as the CD8 ${}^{+}$ GZMB-high T cells subset, which was with the character of high expression of GZMB and low expression of GZMK (Fig. 7A). The following analyses were based upon the three CD8 ${}^{+}$ T cell subsets (Naïve-like, GZMB-high, and GZMK-high) according to the expression of GZMs. We used the BayesPrism algorithm and estimated the fractions of cell subsets in the TCGA-LUAD cohort based on the bulk RNA-seq data. Because of the low fraction and the proliferation of CD8 ${}^{+}$ proliferating T cells, which may affect the accuracy of the deconvolution algorithm, we removed this subset in the analysis of estimating cell fractions. Patients with higher CD8 ${}^{+}$ GZMK high T cells had significantly better prognosis (HR $=$ 0.65, 95% CI [0.48, 0.89], $p=$ 0.006), so did the CD8 ${}^{+}$ GZMB high T cells (HR $=$ 0.67, 95% CI [0.49, 0.92], $p=$ 0.014) and the Naive-like subset (HR $=$ 0.73, 95% CI [0.54, 0.97], $p=$ 0.033). The differential genes among GZMs-based CD8 ${}^{+}$ T cell subsets were shown in Fig. 7D–F. The GNLY, ZNF683, and CXCL13 were highly expressed in the GZMB-high subset, while the LMNA, DUSP2, IGKC, and GPR183 were highly expressed in the GZMK-high subset. The differential genes were ranked by log ${}_2$ (Fold Change) between groups and were inputted into the GSEA analysis (Fig. 7G–I). The GZMK-high subset tended to be associated with complement activating, positive regulation of B cell activation, and immune response compared to the GZMB-high subset. The GZMB-high subset played its roles in homotypic cell-cell adhesion, ATP synthesis coupled electron transport, and hydrogen peroxide catabolic process. Compared to the Naïve-like subset, the GZMK-high, and GZMB-high subsets both presented their functions in antigen processing and presentation. The Naïve-like subset was involved in the secretion and transport of arachidonate, cytoplasmic translation, calcium-dependent cell-cell adhesion via plasma membrane cell adhesion molecules, etc.

Table 1

The information of regulon RUNX3/ETS1 and part of their targets

TF	Gene	highConfAnnot	nMotifs	Best motif	NES	spearCor
ETS1	CD69	TRUE	30	homer__AACAGGAAGT_Ets1-distal	6.63	0.32
	CXCL13	TRUE	4	hocomoco__ETS1_HUMAN.H11MO.0.A	4.01	0.38
	GZMA	TRUE	23	factorbook__ETS1	6	0.60
	GZMB	TRUE	10	factorbook__ETS1	6	0.37
	GZMH	TRUE	1	tfdimers__MD00009	3.07	0.43
	GZMK	TRUE	13	cisbp__M1922	5.66	0.57
	ITGA1	TRUE	39	homer__AACAGGAAGT_Ets1-distal	6.63	0.20
	PRDM1	TRUE	42	homer__AACAGGAAGT_Ets1-distal	6.63	0.07
	ZNF683	TRUE	12	cisbp__M1922	5.66	0.39
RUNX3	GZMA	TRUE	54	transfac_pro__M08866	12.6	0.68
	GZMB	TRUE	75	transfac_pro__M08866	12.6	0.55
	GZMH	TRUE	46	transfac_pro__M08866	12.6	0.57
	GZMK	FALSE	1	hocomoco__PEBB_HUMAN.H11MO.0.C	3.29	0.63
	GZMM	TRUE	22	dbcorrdb__RUNX3__ENCSR000BRI_1__m1	8.7	0.67
	ZNF683	TRUE	1	tfdimers__MD00496	3.46	0.53

The table presents the number of motifs, the best motif, and the spearman correlation coefficients between the transcription factor and target genes.

Figure 8.

The analysis of regulons and the correlations between the RUNX3 and the CD8 ${}^{+}$ T cell subsets. (A) scRNA-Seq data from all cells (average by 20 cells per dot) were projected onto the t-SNE. (B)–(E) Rank of regulons in CD8 ${}^{+}$ Naïve-like subset, CD8 ${}^{+}$ GZMK-high subset, CD8 ${}^{+}$ Exhausted subset, and CD8 ${}^{+}$ Cytotoxic subset based on regulon specificity score (RSS), respectively. (F) The activity of regulon RUNX3 and ETS1 in different cell types based on the t-SNE. (G) The expression of RUNX3/ETS1 gene in different immune cells. (H) The correlations between the expression of RUNX3/ETS1 gene and the infiltration of CD8 ${}^{+}$ T cell subsets. (I) Enrichment analysis based on the genes from the RUNX3 regulon.

We further performed the scMetabolism method to identify the heterogeneity of CD8 ${}^{+}$ T cell subsets in metabolism (Fig. S4A). The activity of TCA cycle, oxidative phosphorylation pathways, and Glycine, serine, and threonine metabolism was higher in the GZMB-high subset (Fig. S4B–D). The GZMK-high subset was associated with higher activity of one carbon pool by folate and part of the lipid metabolism pathways (Fig. S4E). The tryptophan metabolism pathway was up-regulated in both GZMB-high and Naïve-like subsets (Fig. S4F). Consistent with the results before, the Naïve-like subset was associated with arachidonic acid metabolism (Fig. S4G).

3.10 Inference and analysis of cell-cell communications

We utilized the package CellChat to estimate the communications between GZMs-based CD8 ${}^{+}$ T cell subsets and other cells. We identified the significant interactions (L-R pairs) and visualized the differential interactions from three CD8 ${}^{+}$ T cell subsets to our target cell clusters (Fig. 6D). In the L-R pairs level, the MIF-(CD74 $+$ CXCR4) was identified from GZMB-high subsets and Naive-like subsets to CD4 $+$ T cells, Follicular Helper T cells, $\gamma \delta$ T cells, B cells, and NK cells. However, this interaction was lacking in the GZMK-high subset. The CLEC2D-KLRB1 pair had a stronger connection from the GZMB-high subset to the $\gamma \delta$ T cells and NK cells compared to the other two subsets.

3.11 Cell-type-specific regulon activity analysis

We conducted the gene regulatory network and identified the cell state in all cells by SCENIC. We systematically identified critical transcription factors for cell identity and defined a regulon specificity score (RSS) based on Jensen-Shannon divergence [42] (Fig. 8B–E). The results of t-SNE dimensional reduction based on the regulons were presented in Fig. 8A, displaying the distinctness of CD8 ${}^{+}$ T cell subsets in the activity of regulons. We found that the CREM, THAP11, RUNX3, and SRF transcription factors with corresponding regulons had the highest activity in CD8 ${}^{+}$ Naive T cells (Fig. 8B). The GZMK-high subset was with the high activity of SRF, RUNX3, THAP11, CREM, and ETS1. CD8 ${}^{+}$ Cytotoxic T cells had the highest RSS of ASCL2.

We then searched for the transcription factors that regulate the expression of GZMs, which might furtherly affect the fractions of GZMK-high and GZMB-high subsets. Thronging the SCENIC method, we identified the regulon RUNX3 consisted of GZMA/B/M/H while the regulon ETS1 included the GZMA/B/H/K and both of them included the key T ${}_{\text{RM}}$ hallmark gene ZNF683, which presented noticeable and positive correlativity (Spearman’s corr. $>$ 0.3, Table 1). Compared to the RUNX3, the regulon ETS1 also consisted of CD69, ITGA1, PRDM1, and CXCL13, which encoded the proteins related to the tissue-resident memory programs (Table 1). The method used the GENIE3 to infer the co-expression modules between transcription factors and candidate target genes. Then the RcisTarget identified those modules where the regulator’s binding motif is significantly enriched across the target genes and creates regulons with only direct targets. So the RUNX3 had both positive correlations and direct binding motifs to GZMA/B/M/H but not GZMK, which might cause the heterogeneity between the GZMB-high subset and the GZMK-subset. Both the regulon RUNX3 and regulon ETS1 had higher activity in the CD8 ${}^{+}$ T cells and CD4 ${}^{+}$ T cells (Fig. 8F). The RUNX3 and ETS1 gene were mainly expressed in CD8 ${}^{+}$ T cells, but RUNX3 was also highly expressed in $\gamma \delta$ T cells while ETS1 was expressed in endothelial cells (Fig. 8G). The correlation analysis was conducted between the expression of RUNX3/ETS1 and the fraction of GZMK-high/GZMB-high subsets in the TCGA-LUAD cohort, presenting significantly positive correlations, respectively (Fig. 8H). However, it did not show correlations between these two genes and the fractions of Naïve-like subsets (Fig. S5B). We used the genes from the regulon RUNX3/ETS1 to further analyze their functions. The results showed the roles of RUNX3 in activating CD8 ${}^{+}$ T cells or CD4 ${}^{+}$ T cells, regulation of T cell differentiation, mediated immunity cytotoxicity effector, and positive regulation of cell-cell adhesion (Fig. 8I). As for the regulon ETS1, it was also involved in cell-cell adhesion, in addition to mononuclear leukocyte lymphocyte proliferation, and $\alpha -\beta$ mononuclear cell differentiation (Fig. S5D).

4. Discussion

Pan-cancer analysis can reveal the similarities and heterogeneities of different cancers, which may provide theoretical support for cancer therapeutic, target design, and potential therapeutic drug screening. In the present study, we performed a comprehensive bioinformatics analysis of the GZMs and their associations with immune infiltration, biological processes, immunotherapy response, and prognostic value, which vary significantly among different cancers.

Both genetic mutation and epigenetic modification can induce abnormal gene expression during tumorigenesis. We found that gene methylation levels can lead to abnormally low GZMs expression. However, the correlations of CNV with GZMs expression vary from different cancers. Furthermore, we calculated the GZMs score based on granzymes gene expressions in pan-cancer and indicated its potentially predictive performance in cancer immunotherapy, which was validated in eight immunotherapy datasets. We further analyzed why the GZMs score was associated with different prognoses, such as tumor microenvironment (TME), differential pathways, and drug sensitivity. We also illustrated the heterogeneity of GZMs-based CD8 ${}^{+}$ T cell subsets in differential genes, functions, metabolisms, and their effects on the TME in LUAD. Finally, we identified the potential target transcription factors RUNX3 and ETS1, which may affect the cell infiltration and the heterogeneity of GZMs-based CD8 ${}^{+}$ T cell subsets, providing a new perspective for cancer therapeutic strategies.

More and more studies have established the core position of immunotherapy in altering the therapeutic landscape for many cancers. Still, the variability of immunotherapy response highlights the need for identifying predictive biomarkers [43]. A recent study demonstrated that TMB, in concert with PD-L1 expression was a valuable biomarker for ICB selection across some cancer types [44]. The role of granzyme B as a predictive biomarker of immunotherapy was demonstrated before [45, 46]. Therefore, we further explored the predictive values of GZMs for immunotherapy efficacy. The IPS and TIDE scores showed possible predictive value of GZMs score to immunotherapy efficacy in pan-cancer, but it should be validated on a large scale because of the heterogeneity of different cancer types. Our analysis of immunotherapy datasets demonstrated that the GZMs score could predict immunotherapy response in RCC, urothelial carcinoma, and melanoma. However, RCC patients in the high GZMs score group had a worse prognosis. This may be because the GZMs score was strongly positively associated with the infiltration of CD8 ${}^{+}$ T cells, and previous research had reported that in RCC, baseline CD8 ${}^{+}$ T cell infiltration was associated with a worse prognosis [47]. The TME of RCC is different from that of other tumor types, and the subgroup of RCC patients with elevated CD8 ${}^{+}$ T-cell infiltration have particularly poor outcomes, with increased expression of immune checkpoint molecules and a paucity of functional antigen-presenting cells [48].

Furthermore, we explored possible explanations that caused the different prognoses in different cancers. Firstly, the GZMs score could partly represent the infiltrating level of effector CD8 ${}^{+}$ T cells, while the expression of checkpoint molecules of CD8 ${}^{+}$ T cells and the existence of antigen-presenting cells exhibited heterogeneity, such as RCC. The activation level of immune-related pathways (T cell differentiation involved in immune response, cellular defense response, interferon-alpha response, monocyte chemotaxis, lymphocyte chemotaxis, Chemokine signaling pathway, etc.) and the inhibition level of oncogenic pathways may play crucial roles in exerting anti-tumor activity even affecting prognosis. For example, abnormal activation of the Hedgehog pathway was responsible for tumorigenesis and maintenance of multiple cancers [49], which was negatively correlated to GZMs score in most cancer types. The microtubule was one of the critical components of the cytoskeleton and played an essential role in the maintenance of cell shape and the process of signal transduction and mitosis [50]. However, the microtubule-related pathways, the cell cycle pathway, were positively related to the GZMs scores in UVM, which might cause a worse prognosis. In addition, the UVM, and LAML showed apparent differences in enrichment pathways compared to other cancer types such as the mTOR signaling, which promotes cell proliferation and metabolism that contribute to tumor initiation and progression upon hyper-activation [51]. The mTOR signaling pathway was relatively activated in UVM and LAML patients with higher GZMs scores. Similarly, the activation of the Wnt signaling pathway played a crucial role in cell differentiation, polarization, and migration [52]. Its activation may also cause a worse prognosis in LAML.

BASSEZ A [53] had reported PD1-expressing T cells clonally expanded upon anti-PD1 treatment in patients with hormone receptor-positive or triple-negative breast cancer. The expansion was mainly involved in CD8 ${}^{+}$ T cells with pronounced expression of cytotoxic activity (PRF1, GZMB) but undifferentiated pre-effector/memory T cells (TCF7 $+$ , GZMK $+$ ) were inversely correlated with T cell expansion [53]. In clear cell renal cell carcinoma, CD8 ${}^{+}$ T cells upregulate GZMB following nivolumab in responders [54]. Baolin Liu et al. [55] observed increased levels of precursor-exhausted T cells in responsive tumors after anti-PD-1 treatment, characterized by low expression of co-inhibitory molecules and high expression of GZMK, while nonresponsive tumors failed to accumulate precursor exhausted T cells in NSCLC. Based on these studies of the phenotypic properties of different GZMs, we analyzed three patterns of GZMs expressions in CD8 ${}^{+}$ T cells. GZMB and GZMK expressions presented a measure of exclusiveness. In the CD8 ${}^{+}$ T cells differentiation and exhaustion trajectories, GZMK, GZMH, and GZMM increased halfway through the trajectory but then decreased while GZMB increased towards the end of the trajectory (Fig. 6C, Fig. S3E–G). These results indicated that the GZMs-based CD8 ${}^{+}$ T cell subsets may have different responses to immunotherapy. Then we supplemented some contents of the discrepancy between different subsets and tried to find the potential reasons that affect the response to immunotherapy.

Subsequently, we discovered the GZMB-high subset characterized by the higher activity of the TCA cycle, oxidative phosphorylation pathways (OXPHOS), and pentose phosphate pathway (PPP). Previous research demonstrated that upregulated TCA cycle metabolism and OXPHOS are critical aspects of CD8 ${}^{+}$ T cell activation [56]. The PPP is the primary cellular source for NADPH, which is required for fatty acid and plasma membrane synthesis in newly activated CD8 ${}^{+}$ T cells [57]. Both the GZMB-high and Naïve-like subsets displayed the higher activity of tryptophan metabolism, whose availability within the TME is a crucial factor in determining the strength and quality of the T cell response [56]. Lin R et al. [58] demonstrated that reprogramming toward fatty acid oxidation also prolongs TRM cell longevity and the GZMB-high subset presented higher activity of fatty acid degradation, indicating the potential relations between the lipid metabolism and the tissue-resident memory programs. As for the GZMK-high subset, it was associated with complement activation, One carbon pool by folate, arginine, and proline metabolism. Roger Geiger et al. [59] found critical changes in the arginine metabolism that led to a drop in intracellular L-arginine concentration, which directly impacts the metabolic fitness and survival capacity of T cells that are crucial for anti-tumor responses.

In the CellChat analysis, the three subsets also showed discrepancies in MIF- (CD74 $+$ CXCR4), MIF- (CD74 $+$ CD44), and CLEC2D-KLRB1 pairs. Macrophage migration inhibitory factor (MIF) exhibited chemokine-like properties and was identified as a non-cognate ligand of the chemokine receptors CXCR2 and CXCR4 [60, 61]. There was evidence that the surface interaction of MIF with CD74 leads to the activation of ERK1/2-MAPK and AKT/PI3K signaling, and the former was involved in regulating proliferation, differentiation, and transformation while the latter was involved in cell growth, cell cycle entry, cell migration, and cell survival [62, 63, 64]. CD74 is also vital for the regulation of B cell survival by MIF [65]. MIF-triggered ERK1/2 signaling requires the recruitment of CD44 and a non-receptor tyrosine kinase [66, 67]. These L-R pairs indicated the potential effects that the GZMB-high subset associated with the proliferation, differentiation, and survival of immune cells. However, these pairs were lacking between the GZMK-high subsets and other cells. For the other L-R pair, NKRP1A also called killer cell lectin-like receptor B1 (KLRB1), clearly functions as an inhibitory receptor on NK cells, and NKRP1A engagement had been reported to costimulate T cell proliferation and cytokine secretion by activated T cells [68, 69, 70], but also to reduce the release of TNF by CD8 T cells [71]. The GZMB-high subset may have more profound effects on the $\gamma \delta$ T cells and NK cells through this connection. Generally, the higher fractions of all three CD8 ${}^{+}$ T cell subsets were associated with a better prognosis, but the mechanisms had discrepancies in metabolism and cell-cell communications.

Moreover, we established a network of regulons that could be postulated as likely candidates for identifying or regulating states of CD8 ${}^{+}$ T cells throughout the different periods of differentiation and exhaustion. Our research found significant positive correlations between the transcription factor RUNX3/ETS1 and the GZMs. Furthermore, the direct binding motif indicated that both the RUNX3 and ETS1 might be the critical factor in affecting the expression of GZMs and fractions of the GZMs-based subsets. Previous research has found that RUNX3 functions as a tumor suppressor and is frequently deleted or transcriptionally silenced in cancer. Milner JJ et al. [72] identified the transcription factor Runx3 as a key regulator of CD8 ${}^{+}$ T ${}_{\text{RM}}$ cell differentiation and homeostasis in mice. Besides, Runx3-overexpression enhanced TIL abundance and expression of granzyme B. Increasing the number of T ${}_{\text{RM}}$ cells by overexpression of Runx3, moreover, results in better control of melanoma in mice. T ${}_{\text{RM}}$ cells upregulate their expression of the transcriptional regulators Hobit and Blimp-1 in a manner dependent on T-bet and Runx3, respectively [73]. A recent clinical trial (NCT02919683) demonstrated that in oral cancer patients, the CD8 ${}^{+}$ T cell subsets with high expression of GZMB also expressed typical T ${}_{\text{RM}}$ genes (ITGAE and ZNF683) and the tissue-resident memory and cytotoxicity programs, supporting the capacity for rapid response in neo-adjuvant immunotherapy [74]. This was consistent with our results in LUAD and the RUNX3 may also upregulated the expression of ZNF683 (Spearman’s corr. $=$ 0.529, Table 1). Our results suggested that the RUNX3 may cause both the expression of GZMB, the infiltration of GZMB-high CD8 ${}^{+}$ T cell subset, and the tissue-resident memory program, which may suppress tumors and improve the prognosis in LUAD patients. Similar results also appeared in transcription factor ETS1. ETS1 also regulated the expression of GZMs and other T ${}_{\text{RM}}$ mark genes. The proteins encoded by CD69, ITGA1, PRDM1, ZNF683, and CXCL13 were all involved in the tissue-resident memory programs and be regulated by transcription factor ETS1 (Table 1). CD69 was often used to identify T ${}_{\text{RM}}$ cells [75, 76]. Tissue residence is actively promoted by the expression of adhesion molecules such as the integrin CD49a (encoded by ITGA1), which binds to the extracellular matrix components collagen and laminin [77, 78]. T ${}_{\text{RM}}$ cells expressing CD49a were potent killers among T cells in a mouse model of melanoma [79]. T ${}_{\text{RM}}$ cells also down-regulate molecules associated with tissue egress. Blimp-1 (encoded by PRDM1) and Hobit (encoded by ZNF683) instructed tissue retention and bolstered unresponsiveness to tissue-egress signals by suppressing the expression of S1pr1, Tcf7, and Ccr7 and directly or indirectly promotes transcription of GZMB [80]. Yang L et al. [81] demonstrated that CD8 $+$ T ${}_{\text{RM}}$ cells might promote tertiary lymphoid structure generation through CXCL13 secretion. Previous research did not focus on the relationships among the ETS1, GZMs, and tissue-resident memory programs. According to these results, we inferred that the transcription factor ETS1 controls the effector function of CD8 ${}^{+}$ T cells through up-regulating the expression of GZMs and has a role in the differentiation and maintenance of T ${}_{\text{RM}}$ cells, and its overexpression may enhance the number of these cells. Because of the characteristic of TRM cells that locate in non-lymphoid tissues and their capacity for rapid response displayed, ETS1 has the potential to better control solid tumors in immunotherapy, and our survival analysis also indicated that the higher expression of ETS1 tended to be associated with the better prognosis in TCGA-LUAD cohort (Fig. S5C). The expression of GZMs as the effector molecule collaborated with the tissue residence program caused by RUNX3 or ETS1 might improve the prognosis of LUAD patients.

In conclusion, we conducted a comprehensive assessment of GZMs, unveiling complicated and comprehensive roles of granzyme family members in cancer progression and clinical outcome. We also revealed that the higher GZMs score may indicate better response and OS outcome for immunotherapy in melanoma and urothelial carcinoma but worse OS in RCC, indicating that the GZMs score is a potential prognostic biomarker of diverse cancers. RUNX3 and ETS1 may be the potential targets to regulate the infiltration of GZMs-based CD8 $+$ T cell subsets and affect the tissue-resident memory programs in LUAD, which may affect the prognosis of LUAD patients and the response to immunotherapy.