The single-cell gene expression matrices including raw count and TPM data were obtained from the GEO database: GSE108989 (CRC), GSE99254(HCC), and GSE98638(NSCLC), and we isolated CD8+ T cells from peripheral blood (P), adjacent normal (N), and tumor tissues (T).

The raw count expression matrices were processed by R package Seurat v3.2.0 (http://satijalab.org/seurat/). To filter out the low-quality cells, we excluded cells with <600 and >10,000 detected genes (17). Counts were log normalized and scaled by linear regression against the number of reads with function NormalizeData and ScaleData. The highly variable genes (HVGs) were generated with FindVariableFeatures. Principal component analysis (PCA) was performed on the top 2,000 HVGs using function RunPCA. The appropriate PCs were selected for graph-based clustering with functions FindNeighbors and FindClusters. For visualization of clustering analysis, we performed uniform manifold approximation and projection for dimension reduction (UMAP) using RunUMAP function in Seurat. To eliminate the obvious effect from different patients, we performed standard normalization and variable feature selection after acquiring the data. Next, the function FindIntegrationAnchors was performed to find a set of anchors, which is used to integrate the data by the function IntegrateData. Also, the function AddModuleScore was used to calculate scores of gene list in different cell types.

Differentially expressed genes (DEGs) of each cluster were identified based on Wilcoxon rank-sum test using function FindAllMarkers compared with the rest of the clusters. In brief, for each cluster, only genes that met these criteria were considered cluster-specific DEGs (1): log2FC >0.25; (2) expressed >25% in either of the two groups of cells; (3) adjusted p-value <0.05. The top DEGs were selected to annotate each cluster based on the canonical markers from previous studies; also, the CellMarker databases (18) and R package SingleR v1.4.0 (19) were performed to further improve the accuracy of the annotation.

To explore the potential functional changes of CD8+ T cell of different clusters for each cancer, we performed development trajectory analysis by R package Monocle v2.18.0 (20) with the cluster-specific genes of each cluster. Dimensional reduction and cell ordering were performed using reduceDimension and orderCells functions with default parameter.

Biological significance was explored by Gene Ontology (GO) term enrichment analysis by R package clusterProfiler v3.18.0 (21) including biological process, cellular component, and molecular function. Adjusted p-value <0.05 was considered statistically significant. Visualization is realized by R package ggplot2 v3.3.3 and ggalluvial v0.12.3 (22).

In order to identify the highly linked genes in specific clusters, weighted gene coexpression network analysis (WGCNA) was performed with functions in the R package WGCNA (23). To attenuate the effects of noise and outliers, we constructed pseudocells (24) which were calculated as averages of 10 cells randomly chosen within each cluster. The function pickSoftThreshold was used to calculate the soft power parameter and blockwiseModules to construct coexpression network. Finally, corPvalueStudent was used to mine modules related to specific cell types.

According to the median of gene expression values, cells were divided into high and low groups, then the survival curve was shown using the Kapla–Meier curve with a log-rank test by GEPIA2 (http://gepia2.cancer-pku.cn/) to illustrate the relationship between differential genes and overall patient survival. The terms with p-value <0.05 were identified as significant. Additionally, multiple hypothesis testing (FDR) was performed to the significant p-value using the function p.adjust.

For generated gene set, we performed gene set variation analysis (GSVA) to the data in the present study and validation dataset using R package GSVA (25) with the default parameter, to evaluate the effect of distinguishing Tex cells.

The RNA-seq data were downloaded from the TCGA database (https://portal.gdc.cancer.gov/) and then calculated into TPM value. According to the annotation, samples were grouped into normal and four stages of cancer groups. The differential expression of mRNA in the different groups was performed by Wilcoxon test with p-value <0.05.

The protein expression level in normal and cancer tissues was analyzed using the Human Protein Atlas (HPA) database (https://www.proteinatlas.org). Immunohistochemistry pictures were downloaded from the Tissue Atlas and Pathology Atlas.

Immunohistochemistry was carried out on human liver cancer tissue microarray (Shanghai Outdo Biotech Co., Ltd., Shanghai, China). The tissue sections were first dried at 63°C for 1 h, dewaxed, and rehydrated before epitope retrieval by heating at 100°C in 10 mM sodium citrate (pH 6.0) for 5 min in EDTA solution for 20 min. The sections were cooled down to room temperature for 30 min. The tissue sections were treated with 3% hydrogen peroxide for 20 min to eliminate the endogenous peroxidase and alkaline phosphatase activity in the tissue. After cooling down to room temperature, the sections were treated by blocking agents for 10 min. The sections were incubated with individual primary antibody (RGS1, Invitrogen, Waltham, MA, USA; Product #PA5-86730) diluted 1:1,000,800 overnight at 4°C, followed by secondary antibodies at room temperature for 30 min. 3,3′-Diaminobenzidine (DAB) was then applied as a substrate to reveal the antigen. Hematoxylin was used for counterstaining. The stained images were counted by ImageJ software (26), and the optical density (OD) value was used for quantification.

We obtained single-cell transcriptome data of human T cells from the GEO database, including 12 patients from CRC (27), six patients from HCC (28), and 14 patients from NSCLC (29). After strict quality control and filtration, we collected 4,010, 1,752, and 4,439 CD8+ T cells from peripheral blood (P), adjacent normal (N), and tumor tissues (T) ( Supplementary Figures S1A, B ; Supplementary Table S1 ).

We divided the cells of each cancer into different cell types annotated with cluster-specific gene expression ( Figure 1A ; Supplementary Table S2 ). Specifically, in CRC, cells were identified as naïve T cells (Tn, cluster 4, marked with selectin L (SELL), lymphoid enhancer binding factor 1 (LEF1), C–C motif chemokine receptor 7 (CCR7), transcription factor 7 (TCF7), effector T cells (Teff, cluster 2 and cluster 5, marked with C-X3-C motif chemokine receptor 1 (CX3CR1), killer cell lectin-like receptor F1 (KLRF1)), fibroblast growth factor-binding protein 2 (FGFBP2), Fc fragment of IgG receptor IIIa (FCGR3A), exhausted T cells (Tex, cluster 3 and cluster 6, marked with hepatitis A virus cellular receptor 2 (HAVCR2), programmed cell death 1 (PDCD1), lymphocyte activating 3 (LAG3), TOX, CXCL13, tissue-resident memory T cells (Trm, cluster 1, marked with CD69, integrin subunit alpha E (ITGAE)), mucosal-associated invariant T cells (MAIT, cluster 7, marked with solute carrier family 4 member 10 (SLC4A10)), RAR-related orphan receptor C (RORC) ( Figures 1B, C ). Cluster 0 was located between cluster 2 and cluster 3, which was represented as Teff and Tex cells, besides, the genes marked Teff and Tex cells showed relative high expression levels in this cluster, therefore, we identified it as Pre_exhaust T cells. Similarly, cells in HCC were identified as Tn (cluster 4), Teff (cluster 0), Pre_exhaust (cluster 1), Tex (cluster 3 and cluster5), MAIT (cluster 2), and cells in NSCLC that were identified as Tn (cluster 4), Teff (cluster 0 and cluster 6), Pre_exhaust (cluster 1 and cluster 7), Tex (cluster 3), Trm (cluster 2), and MAIT (cluster 5). Different cell types showed preference in different tissues. In general, Tn and Teff cells were mainly enriched in peripheral blood, Pre_exhaust and Tex cells were mainly enriched in tumor tissues, and Trm cells existed more in adjacent normal tissues of CRC and HCC, and in tumor tissues of NSCLC, which may be related to tissue specificity. Also, the number of MAIT cells was relative less than the others. Additionally, we performed exhausted cell scoring by the function AddModuleScore in R package Seurat using the exhaustion gene list including HAVCR2, T-cell immunoreceptor with Ig and ITIM domains (TIGIT), LAG3, PDCD1, CXCL13, layilin (LAYN), TOX, cytotoxic T-lymphocyte associated protein 4 (CTLA4), and B and T lymphocyte associated (BTLA) and visualized in Supplementary Figure S1C , showing the accuracy of Tex cell classification ( Supplementary Figure S1C ).

To further explore the function of each cluster, we performed GO enrichment analysis using the cluster-specific genes ( Figure 1D ). In all three cancers, Tn cells were enriched in T-cell activation, T-cell differentiation, and ribosome biogenesis; Teff cells were enriched in cellular defense response, positive regulation of cytokine production, and T-cell-mediated cytotoxicity; Tex cells were enriched in negative regulation of immune system process, T-cell apoptotic process, and response to hypoxia; Trm cells were enriched in antigen receptor-mediated signaling pathway and leukocyte chemotaxis; MAIT cells were enriched in interferon-gamma production; and Pre_exhaust cells were enriched in terms associated with Teff and Tex cells, including cellular defense response and negative regulation of immune system process, showing the transitional characteristics during T-cell exhaustion.

Tex cells, as shown above, played a negative role in immune system process, and it was suggested that Pre_exhaust cells were a transitional stage from Teff to Tex cells. To find out the highly linked genes associated with T-cell exhaustion, we used the R package WGCNA to construct the weighted coexpression network in Teff, Pre_exhaust, and Tex cells ( Figure 2A ). The blue module (R-value 0.96, p-value 2e−145), turquoise module (R-value 0.96, p-value 7e−67), and blue module (R-value 0.96, p-value 2e−172) represented Tex cells in CRC, HCC, and NSCLC, respectively. Combining with the differentially expressed genes (DEGs) upregulated in the Tex vs. Teff cell comparison in three cancers ( Supplementary Table S3 ), 35 genes were found in three Tex cell-related modules and overexpressed in Tex cells ( Figure 2B ; Table 1 ), which were defined as “Candidate” gene set, including exhaustion marker PDCD1, CTLA4, HAVCR2, TOX, and TIGIT. It is noteworthy that the housekeeping gene glyceraldehyde-3-phosphate dehydrogenase (GAPDH), which can catalyze an important energy-yielding step in glycolysis metabolism, was also upregulated and enriched in the highly linked DEGs, thus we speculated that the glycolysis progress was disordered in Tex cells.

As shown in Figure 2C , poor overall survival was correlated with higher Candidate gene set expression in multiple cancers including liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), head and neck squamous cell carcinoma (HNSC), kidney chromophobe (KICH), skin cutaneous melanoma (SKCM), thymoma (THYM). In addition, the GSVA scores of Candidate gene set in Tex cells were much higher than those of the other clusters in CRC, HCC, and NSCLC, so were in the other independent validation datasets ( Figure 2D ). GSE116390 dataset (30) was composed of four types of tumor-infiltrating CD8+ T cells, which were exhausted, memory-like, naïve, and effector memory-like (EM-like) subsets, from B16 melanoma tumor-bearing mice, and GSE123813 dataset (31) contained six types of tumor-infiltrating CD8+ T cells, which were from 11 patients with advanced basal cell carcinoma, showing higher score in the exhausted related clusters.

As indicated above, the Candidate gene set was enriched in exhausted CD8+ T cells with poor prognosis and was able to distinguish Tex cells from the other CD8+ T cells in different cancers, indicating that the GSVA score of these 35 genes might be an effective prognostic marker or a marker to identify Tex cells.

The exhaustion of T cells was gradually formed, so we performed differential expression analysis of Pre_exhaust and Tex cells compared with Teff cells, in order to further explore the critical genes associated with the formation and development of T-cell exhaustion. In the Pre_exhaust vs. Teff cell comparison, 119 DEGs existed in all three cancers, 57 of whom were upregulated and the others were downregulated. Correspondingly, there were 162 up-regulated and 88 down-regulated DEGs in the Tex vs. Teff cells comparison. Furthermore, there were 40 DEGs overexpressed in both Pre_exhaust and Tex cells compared with Teff cells ( Figure 3A ), which may contribute to the origin of T-cell exhaustion, including the canonical exhaustion marker PDCD1 which encodes an inhibitory receptor (32), CD69 with the capability to mediate the cell retention (33), Cbl Proto-Oncogene B(CBLB) whose deletion can inhibit CD8+ T-cell exhaustion (34), hypoxia-inducible factor-1 (HIF1A) which can stably be expressed in hypoxia condition, consistent with the biological process of Tex cells ( Figure 1D ).

Obviously, regulator of G protein signaling 1 (RGS1) showed the greatest fold change in Tex cells across three cancers and also across different patients, showing the greatest fold change in Pre_exhaust cell in HCC and NSCLC ( Figure 3A ; Supplementary Figure S2 ), illustrating its potential roles during T-cell exhaustion. RGS1 was expressed highly in whole Tex cells compared with Teff cells ( Figure 3B ), eliminating that the high fold change of RGS1 in Pre_exhaust and Tex cells was not caused by partial cells with abnormally high expression value but high expression in whole cells. In order to evaluate the role of RGS1 in tumorigenesis, we analyzed the expression levels of RGS1 between tumor and normal tissues in the TCGA database ( Figure 3C ). RGS1 expression was significantly upregulated in multiple cancers including BRCA, cholangiocarcinoma (CHOL), esophageal carcinoma (ESCA), glioblastoma multiforme (GBM), HNSC, kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), LUAD, rectal carcinoma (READ), and stomach adenocarcinoma (STAD) when compared with the normal samples. We also added and compared the expression levels of RGS1 across different cancer stages and found that RGS1 expression was associated with stage in some cancer types, such as RGS1 expression was higher in stages II, III, and IV vs. stage I in STAD. This result revealed that RGS1 was likely a key tumorigenesis regulator in multiple cancers and may be associated with prognosis. Conspicuously, the prognosis analysis was analyzed in 33 TCGA cancer types ( Supplementary Table S4 ). RGS1 expression was significantly correlated to poor prognosis in seven cancers, including LIHC, adrenocortical carcinoma (ACC), pancreatic adenocarcinoma (PAAD), sarcoma (SARC), SKCM, STAD, and THYM ( Figure 3D ), suggesting that RGS1 was a potential prognostic factor in the survival of the above cancers. Apart from that, the protein level of RGS1 in HPA database showed the immunohistochemical (IHC) staining of RGS1 was negative staining in normal tissues and positive in liver cancer tissues, demonstrating that RGS1 was significantly expressed in cancer tissues than in normal liver tissues ( Figure 3E ). Additionally, we performed IHC verification to quantify the protein expression of RGS1 using the local clinical samples of liver cancer and normal tissues ( Figures 3F, G ), similarly, RGS1 protein displayed stronger staining in hepatocarcinoma, in line with the statistical result (p = 7.1e−5). The clinicopathological information of the patient samples and protein expression value are provided in Supplementary Table S5 . These results showed that RGS1 was highly expressed in Tex cells in cancers, upregulated in tumor tissues in mRNA and protein level, and with poor prognosis in multiple cancers, which indicated its potential key role in T-cell exhaustion or cancer progress, and RGS1 might be an effective prognostic marker or a marker to identify Tex cells.

To discover the relationship between RGS1 and Candidate gene set as their potential roles in exhausted T cells, we calculated the correlation coefficient in single cell and tissue level ( Figure 4 ). In three-cancer single-cell datasets, it was obvious that RGS1 showed positive correlation with 35 genes (up to 0.3~0.8), indicating that consistency of coexpression patterns between these genes. In the bulk RNA sequencing datasets, the correlation coefficients were almost positive, especially high with PDCD1, CTLA4, TNF receptor superfamily member 9 (TNFRSF9), HAVCR2, TOX, and TIGIT, which further confirmed the potential key roles of RGS1 in Tex cells. The negative correlations were mainly with the genes involved in cell cycle and DNA replication, including minichromosome maintenance complex component 7 (MCM7), enolase superfamily member 1 (ENOSF1), minichromosome maintenance complex component 5 (MCM5), minichromosome maintenance complex component 3 (MCM3), proliferating cell nuclear antigen (PCNA), which pointed out the different results observed between single cell and bulk RNA sequencing analysis, suggesting we could obtain deeper understanding about T-cell exhaustion by single-cell sequencing. One interpretation of the above inconsistent result is that the bulk transcriptome change (which mixes many immune and nonimmune cells together) likely reflects an overall expression in these cells and does not discriminate specific T-cell states, which further highlights the advantages in defining and studying T-cell exhaustion state by using single-cell sequencing data.