- Research
- Open access
- Published:
Single-cell combined bioinformatics analysis: construction of immune cluster and risk prognostic model in kidney renal clear cells based on CD8+ T cell-associated genes
European Journal of Medical Research volume 29, Article number: 89 (2024)
Abstract
Background
Kidney cancer is an immunogenic solid tumor, characterized by high tumor burden and infiltration of CD8+ T cells. Although immunotherapy targeting the PD1/CTLA-4 axis has demonstrated excellent clinical efficacy, clinical outcomes in most patients are poor.
Methods
We used the RNA sequencing data from the GEO database for KIRC GSE121636 and normal kidney tissue GSE131685, and performed single-cell analysis for cluster identification, pathway enrichment, and CD8+ T cell-associated gene identification. Subsequently, the significance of different CD8+ T-cell associated gene subtypes was elucidated by consensus clustering, pathway analysis, mutated gene analysis, and KIRC immune microenvironment analysis in the TCGA–KIRC disease cohort. Single gene analysis identified LAG3 as the most critical CD8+ T-cell-associated gene and its function was verified by cell phenotype and immunohistochemistry in KIRC.
Results
In the present study, CD8+ T-cell associated genes in KIRC were screened, including GZMK, CD27, CCL4L2, FXYD2, LAG3, RGS1, CST7, DUSP4, CD8A, and TRBV20-1 and an immunological risk prognostic model was constructed (risk score = − 0.291858656434841*GZMK − 0.192758342489394*FXYD2 + 0.625023643446193*LAG3 + 0.161324477181591*RGS1 − 0.380169045328895*DUSP4 − 0.107221347575037*TRBV20-1). LAG3 was identified and proved as the most critical CD8+ T cell-associated gene in KIRC.
Conclusion
We proposed and constructed an immunological risk prognostic model for CD8+ T cell-associated genes and identified LAG3 as a pivotal gene for KIRC progression and CD8+ T-cell infiltration. The model comprehensively explained the immune microenvironment and provided novel immune-related therapeutic targets and biomarkers in KIRC.
Introduction
Kidney cancer is a common malignancy of the urinary system that originates from the renal tubular epithelium. The most common histological subtype of kidney cancer is renal cell carcinoma (RCC), accounting for approximately 90% [1, 2], including kidney renal clear cell carcinoma (KIRC, 70%), papillary RCC (10–15%), and chromophobe RCC (5%) [2]. The Global Cancer Statistics (2020) stated that RCC accounted for 2.2% of newly diagnosed cancers annually, of which 25 to 30% of patients were diagnosed as advanced or metastatic with a 5-year survival rate of 10%, and 20% to 30% of patients having a propensity for recurrence and metastasis after local operation [3,4,5].
Depending on tumor immunogenic characteristics, the systemic treatment of RCC has witnessed significant changes in the last 20 years [6]. Traditional immunotherapy was predominantly based on interferon (IFN)-α and interleukin (IL)-2; however, IFN-α exhibits poor efficacy and IL-2 displays high toxicity [7, 8]. The subsequent application of vascular endothelial growth factor (VEGF), tyrosine kinase inhibitors (TKIs), and mammalian target of rapamycin (mTOR) inhibitors has improved the efficacy and safety of RCC systemic treatment [6]. Recently, immunotherapy agents targeting programmed death-1 (PD-1)/programmed death-ligand 1 (PD-L1) axis alone or a combination with anti-cytotoxic T lymphocyte-associated protein 4 (CTLA-4) monoclonal antibodies or antiangiogenic agents has greatly expanded the systemic treatment options for RCC [6, 9].
Although the application of immune checkpoint inhibitors (ICIs) as tumor therapeutics has led to major improvements in the RCC systemic treatment, most patients fail to achieve a durable complete response (CR). This could be because RCC is significantly different from other solid tumors in immunogenic features, has a high mutational burden and CD8+ T cell infiltration, and is associated with poor prognosis [10, 11]. CD8+ T cells constitute the major anti-tumor effective cells in the tumor microenvironment (TME) and exert cytotoxic effects. However, their function is impaired by immunosuppressive cells or molecules in the TME [12]. Up-regulation of co-suppressor molecules, including PD1 and CTLA-4, on the surface of CD8+ T cells bind to relevant ligands, ultimately causing CD8+ T cell exhaustion [13, 14]. Therefore, blocking PD-1-mediated inhibitory signaling by monoclonal antibodies could reverse CD8+ T cell exhaustion, thereby hindering tumor progression. However, this contradicted CD8+ T cell characteristics in RCC [10, 11, 15]. Recent studies have reported that the timing of PD-1 inhibition could negatively affect T-cell priming and memory CD8+ T cell formation, thus contributing to more appropriate timings in RCC immunotherapy [16, 17].
No relevant studies are available to explore the function of CD8+ T cell-associated gene sets in the KIRC and its TME. The high accuracy and specificity of single-cell RNA sequencing allows analysis of the functional status of CD8+ T cells and the expression of its associated genes using the single-cell sequencing data of immune-infiltrated KIRC (GSE121636 [18]) and normal kidney (GSE131685 [19]). In the present study, we developed a CD8+ T cell-associated immunological risk prognostic model using CD8+ T cell-associated marker genes obtained from single-cell sequencing analysis to predict the survival status, tumor immune microenvironment, and immunotherapy responsiveness of KIRC patients, thereby providing a potential target and predictive evidence for immunotherapy.
Materials
Single-cell data filtration, pre-processing, and cluster identification
The single-cell sequencing data of KIRC GSE121636 [18] and normal kidney GSE131685 [19] were screened using the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). See Table 1 for details. The Seurat package was used to generate objects and filter out poor-quality cells. The standard data pre-processing processes were performed and percentages of gene numbers, cell counts, and mitochondrial sequencing counts were calculated. The filtering criteria were genes with less than only three cells detected and disregarded cells with less than 200 detected gene numbers. Cells with less than 200 or more than 2500 genes detected and those with high mitochondrial content (> 10%) were filtered out as well. We scaled the UMI counts using a scale factor of 10,000 to normalize the library size effect of each cell. Following the log transformation of data, other factors, including “percent.mt,” “nCount_RNA,” and “nFeature_RNA” were corrected for variation regression using the ScaleData function in Seurat (v3.0.2). The corrected normalized data metric was applied to standard analysis. The top 50 variable genes were extracted for principal component analysis (PCA). The top 10 principal components were retained for UMAP visualization and clustering. Cell clustering was performed using the FindClusters function (resolution = 0.5) included in the Seurat R package.
Differential expression and survival prognosis analysis
The “Survival R package” was used to analyze differential expression, overall survival (OS), disease-specific survival, and progression-free survival of CD8+ T cell-associated genes based on single-cell sequencing and assays. Furthermore, we established the correlation between key genes and clinicopathological features of KIRC and constructed a prognostic nomogram and calibration curve in the TCGA–KIRC cohort.
Consensus clustering
A cluster analysis was performed using the “Consensus Cluster Plus R Package” [20] and agglomerative PAM clustering with a 1-Pearson correlation distances and resampling 80% of the samples for 10 repetitions. The optimal number of clusters was determined using the empirical cumulative distribution function plot.
Identification of differentially expressed genes (DEGs)
Differential mRNA expression between subtypes was studied using the Limma package (version: 3.40.2). The screening criteria were adjusted as p < 0.05 and |fold change| > 2.
Functional enrichment analysis and gene set enrichment analysis (GSEA)
The “clusterProfiler” package was used to evaluate the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of different CD8+ T cell-associated subtypes, with q-value and p-value thresholds of < 0.05. Moreover, the difference in gene sets between high and low expression of CD8+ T cell-associated genes was further assessed using the gene set enrichment analysis (GSEA) software (http://www.broadinstitute.org/gsea/index.jsp).
Somatic mutation and immune landscape analysis
The KIRC somatic mutation data were obtained from the TCGA GDC data portal, and waterfall plots were created using the “Maftools” package. The relative percentages of different immune cell types were determined using the ESTIMATE and xCell algorithms, and the relative percentages of immune cell types between the two subtypes were compared by landscape plots. Moreover, TISIDB (http://cis.hku.hk/TISIDB/index.php) was used to further evaluate the correlation between two subtypes and different immune indicators, including lymphocytes, immune inhibitors, immunostimulators, major histocompatibility (MHC) molecules, chemokines, and chemokine receptors.
Construction of CD8+ T cell-associated risk signature
Regression analysis was performed using the Lasso–Cox method via the glmnet package. A tenfold cross-validation was set up to obtain the optimal model. The prognostic significance of genes involved was assessed using the Cox method, and the relationship between different risk scores and patient follow-up time, events, and changes in the expression of individual genes was analyzed.
Cell lines and cell culture
Human renal cancer cell lines (786-O and ACHN) were cultured in Dulbecco’s modified Eagle medium (DMEM; Gibco Thermo Fisher Scientific, USA), containing 10% fetal bovine serum (Lonsera, Uruguay), and 1% penicillin–streptomycin solution (Keygen, China). All cell lines were purchased from the Shanghai Institutes for Biological Sciences and incubated in 95% humidified air at 37 °C and 5% CO2.
RNA extraction and RT-PCR
The RNA was extracted using the RNA extraction kit (Takara Kusatsu, Japan), and the Hiscript II First-Strand cDNA Synthesis Kit was used to synthesize the complementary DNA (Vazyme, China). Reverse transcriptase-polymerase chain reaction (RT-PCR) was performed using the MonAmpTM SYBR Green qPCR Mix (Monad Biotech, China).
Small interfering RNA
The small interfering RNAs (siRNAs) against the LAG3 gene were designed and synthesized by GenePharma Co. (China).
Cell proliferation and colony formation assays
For the cell proliferation assay, 1000 cells were seeded into 96-well plates for 0 h to 120 h, and 10 µL of the cell counting kit-8 (Keygen, China) solution was added per well. After a 2 h incubation at 37 °C, optical density at 450 nm (OD 450 nm) was measured using a microplate reader (Bio-Tek, USA). For the colony formation assay, cells were seeded into 6-well plates at a density of 1–2 × 10 [3] cells/well and incubated for 10 to 14 days at 37 °C. Next, the cells were washed using phosphate-buffered saline, fixed with 4% polyformaldehyde (Service Bio, China), and stained with 0.1% crystal violet solution (Keygen, China). Colonies containing > 50 cells were counted using the ImageJ 2X software 2.1.4.7 (Rawak Software Inc., Germany).
Wound healing and transwell assay
For the wound healing assay, cells were inoculated into 6-well plates and treated with si-/nc-LAG3. A straight scratch was created on the plate with a sterilized needle tip when the cell density was approximately 70%. The cell wound edge was marked and imaged under a microscope at the starting time point. After 0 h to 24 h, the migrated distance was measured and the wound closure percentage was calculated. For transwell assays, cells were inoculated into a 24-well transwell cell apical chamber containing the matrix gel (BD, USA) for evaluating invasion and gel-free for migration. The bottom and upper chambers contained the RPMI medium and fetal bovine serum-free medium, respectively. Cells invading the bottom chambers were fixed with 4% polyformaldehyde, stained with 0.1% crystal violet solution, counted, and imaged under a microscope.
Tissue samples and tissue microarrays
Formalin-fixed and paraffin-embedded prostate cancer tissue samples were collected from patients who underwent radical nephrectomy in the affiliated Zhongda Hospital of Southeast University, China, from April 2020 to November 2021. The study samples were from patients with KIRC, and the pathological diagnosis was confirmed by at least two pathologists. With the tumor as the center, normal tissues adjacent to the tumor were used as study materials, and two pairs of tissue microarrays were created with a 0.6 mm diameter.
Immunohistochemistry
Formalin-fixed and paraffin-embedded tissue was dewaxed and dehydrated using xylene and serially-diluted ethanol. The tissue sections were incubated at 121 ℃ in an autoclave for 5 min to extract the antigen, following which these were incubated with anti-LAG3-monoclonal antibody at 4 ℃ overnight, and the bound antibody (Proteintech) was incubated at 37 ℃ for 30 min.
Statistical analysis
The statistical analysis was performed using the R software (version 4.0.2). Multivariate Cox regression analyses were used to evaluate the prognostic significance. When p < 0.05 or log-rank p < 0.05, the difference was significant.
Results
Identification of CD8+ T cell-associated gene subgroups by single-cell analysis
We used the GEO database to screen a single-cell sequencing dataset about KIRC GSE121636 [18], including subsets GSE3440844, GSE3440845, and GSE3440846, as well as the normal kidney single-cell sequencing dataset GSE131685 [18], including subsets, namely, GSE4145204 and GSEE4145205. The correlation between the pre-filtered (Fig. 1B) and post-filtered (Fig. 1C) base data was analyzed by excluding cells with less than 200 or more than 2500 genes detected and cells with mitochondrial content > 10% (Fig. 1A). The analysis suggested that the cells were more highly active and that the sequencing depth was not near saturation after data filtration. Gene clusters with a high degree of variation were obtained by the PCA analysis and displayed PCA plots of PC1 and PC2 after combining to remove the batch effects (Fig. 1D), as well as obtain the top 10 highly variable features (HVGs) in 2000 HVGs including S100A8, S100A9, JCHAIN, APOC1, TPSAB1, DEFB1, TPSB2, C1QB, C1QC, and SEPP1 (Fig. 1E).
The HVGs obtained above were included in the PCA analysis and 50 highly variable gene clusters were obtained using the ElbowPlot function and visualized by JackStrawPlot (Fig. 2A). Subsequently, KIRC cells and normal kidney tissue cells were divided into 14 clusters using the UMAP dimensionality degradation method (Fig. 2B). The top 10 genes in 14 clusters (Additional file 1: Table S1) were automatically matched to each cell type in the kidney tissue [21], ultimately obtaining eight different cell subtypes including the epithelium, CD8+ T cells, B cells, mast cells, NKT cells, neutrophils, Tregs, and plasma cells (Fig. 2C). A comparison of the original UMAP cluster with the artificially annotated cell cluster (Fig. 2B–C) revealed that B cells, NKT cells, neutrophils, and plasma cells were differentially expressed in KIRC and normal kidney tissue cells, with NKT cells being relatively highly expressed in KIRC (Fig. 2D). This analysis further revealed that the epithelial cells and mast cells were specific to the normal kidney tissue, whereas Tregs and CD8+ T cells were specific to KIRC, with CD8+ T cells exhibiting the highest expression (Fig. 2E). Moreover, different metabolic pathways between KIRC and normal kidney tissue were analyzed, with the KIRC tissue mostly enriched in DNA_replication, Tumor_proliferation_signature, Tumor_inflammation_signature, and IL-10_anti-inflammation_signaling_pathway, whereas the normal kidney tissue was enriched in MYC_targets, ECM-related_genes, and Cellular_response_to_hyposis (Fig. 2F). The artificially annotated cell clusters were used to study the metabolic pathways of each cell, with CD8+ T cell metabolic functions mostly enriched in Tumor_proliferation_signature, Tumor_inflammation_signature, G2M_checkpoint, and DNA_replication (Fig. 2G), whose metabolic functions overlapped highly with those of KIRC. The above analysis indicated that CD8+ T cells could be an essential immune cell population in the formation, progression, and immune response to KIRC.
Identification of CD8+ T cell-associated subtypes by consensus clustering
The above single-cell sequencing analysis revealed that the top 10 CD8+ T cell-associated genes differentially expressed in KIRC were GZMK, CD27, CCL4L2, FXYD2, LAG3, RGS1, CST7, DUSP4, CD8A, and TRBV20-1 (Additional file 1: Table S1). Subsequently, we verified the differential expression of these genes in the TCGA–KIRC cohort by heat map (Fig. 3A), paired differential expression analysis (Fig. 3B), and unpaired differential expression analysis (Fig. 3C). The results suggested that their expression was statistically significant in KIRC versus paracancerous tissues. Among these, GZMK, CD27, CCL4L2, LAG3, RGS1, CST7, DUSP4, CD8A, and TRBV20-1 were highly expressed in KIRC and were considered oncogenes, whereas FXYD2 was lowly expressed. The differential expression of the above genes in cancer and normal tissues was further validated using the HPA database (https://www.proteinatlas.org), which was consistent with the TCGA–KIRC cohort analysis. However, the database did not have immunohistochemical plots for CCL4L2, DUSP4, and TRBV20-1 (Additional file 3: Figure S1, Additional file 2: Table S2). Moreover, the number of clusters with the highest average consistency within the group was K = 2 (Fig. 3D), and the distribution curve was the greatest at K = 2 (Fig. 3E). Therefore, K = 2 was selected to perform clustering, and the consistent cluster (K = 2) is depicted in Fig. 3F, with 231 KIRC patients in subtype C1 and 276 KIRC patients in subtype C2. Differential expression and OS analyses of the top 10 CD8+ T cell-associated genes in C1 and C2 subtypes further revealed that GZMK, CD27, CCL4L2, LAG3, RGS1, CST7, DUSP4, CD8A, and TRBV20-1 were highly expressed in the C2 subtype, whereas FXYD2 was highly expressed in the C1 subtype (Fig. 3G). The OS was shorter in the C2 subtype compared to the C1 subtype (Fig. 3H). These results were similar to the differential expression analysis of the TCGA–KIRC cohort, indicating that the C2 subtype was the oncogenic group for KIRC.
Identification of differentially expressed genes and signal pathways in different CD8+ T cell-associated subtypes
The above analysis revealed a poor prognosis in the oncogenic subtype (C2), whereas the low expression subtype (C1) of CD8+ T cell-associated genes had a better prognosis. Consequently, we identified key DEGs and signaling pathways in each subtype to understand the molecular mechanisms. A total of 507 abnormally regulated genes were identified, including 15 up-regulated and 492 down-regulated genes (Fig. 4A–B). Because of a few up-regulated genes, we only performed the KEGG analysis on down-regulated genes and identified them to be mostly enriched in cytokine–cytokine receptor interaction, primary immunodeficiency, T cell receptor signaling pathway, cell adhesion molecules (CAMs), antigen processing and presentation, and chemokine signaling pathway (Fig. 4C). The up-regulated genes were largely enriched in glutamate and leukotriene activities, including hypoglycin A gamma-glutamyl transpeptidase activity, leukotriene C4 gamma-glutamyl transferase activity, and glutathione hydrolase activity (Fig. 4D). However, the down-regulated genes majorly corresponded to immune-related activities, including immune system processes, immune responses, and adaptive immune responses (Fig. 4E). Moreover, the GSEA analysis further suggested that C1 and C2 subtypes significantly differed in gene set enrichment in immune cells (T cell receptor signaling pathway, antigen processing and presentation, natural killer cell-mediated cytotoxicity, B cell receptor signaling pathway, and leukocyte transendothelial migration) (Fig. 4F) and oncogenic pathways (VEGF signaling pathway, P53 signaling pathway, DNA replication, apoptosis, and pathways in cancer) (Fig. 4G). These results suggested significant differences in gene expression and related pathway enrichment between C1 and C2 subtypes, which could be the underlying mechanism leading to the different prognosis of KIRC.
Somatic mutations and tumor microenvironment landscape in C1 and C2 subtypes
Different somatic mutation profiles were observed between different CD8+ T cell-associated subtypes, and the top five genes with high mutation frequencies in the C1 subtype were VHL (55.9%), PBRM1 (50.4%), TTN (24.4%), MTOR (11.0%), and BAP1 (9.4%) (Fig. 4H). The top five genes in the C2 subtype with high mutation frequencies were VHL (67.2%), PBRM1 (50.4%), SETD2 (23.2%), TTN (20.8%), and BAP1 (16.8%) (Fig. 4I). This indicated that the mutation frequency in the oncogenic C2 subtype was higher, especially in VHL and BAP1 than in the C1 subtype.
Because the study was based on CD8+ T cells, which could potentially influence the activation or silencing of tumor immune responses, we analyzed the tumor immune microenvironment in two subtypes. The StromaScore, ImmuneScore, ESTIMATEScore, and MicroenvironmentScore were higher in the C2 subtype compared to the C1 subtype (Fig. 5A–C, 5F–G) although the xCell algorithm’s StromaScore subgroup type was not meaningful (Fig. 5E). Next, the infiltration differential map of 64 immune cell types was assessed by xCell with the LM64 signature matrix (Fig. 5D). The quantitative analysis indicated that the aDC, CD8+ T cells, CD8+ Tcm, macrophages, NKTs, Th1 cells, and Th2 cells in the C2 subtype were significantly higher than those in the C1 subtype, whereas monocytes and smooth muscle cells were lower than those in the C1 subtype (Fig. 5H). Subsequently, a comprehensive analysis of the immune gene correlation between different subtypes revealed that the expression of both immune checkpoint-related genes and MHC molecules was significantly higher in the C2 subtype than in the C1 subtype (Fig. 5I–J). For immunostimulators and chemokines, the expression of CD27, CD48, ICOS, CCL5, CXCL9, and CXCL10 in the C2 subtype was significantly higher than in the C1 subtype (Fig. 5K–L).
Construction and validation of CD8+ T cell-associated risk signature
We used the LASSO regression analysis and Cox univariate analysis to construct a risk prognostic model based on CD8+ T cell-associated genes. After integrating survival time, survival status, and gene expression data, the final six-gene risk-score model was developed premised on 10 CD8+ T cell-associated genes on the LASSO regression analysis using the following model equation: RiskScore = − 0.291858656434841*GZMK − 0.192758342489394*FXYD2 + 0.625023643446193*LAG3 + 0.161324477181591*RGS1 − 0.380169045328895*DUSP4 − 0.107221347575037*TRBV20-1 (Fig. 6A–B). Subsequently, it was finally determined by the Cox univariate analysis that GZMK, FXYD2, LAG3, RGS1, and DUSP4 displayed prognostic significance (Fig. 6C). Next, the relationship between survival status and the five genes above was investigated. Our results demonstrated that the number of survival states was considerably higher in the anti-cancer cohort (C1 subtype) compared to the oncogenic cohort (C2 subtype) (Fig. 6D). The risk scores for the above five genes were further quantified using the KM analysis in the TCGA–KIRC cohort, where high-risk scores corresponded to poorer OS (Fig. 6E) and 1, 3, and 5-year survival rates of 0.67, 0.67, and 0.69, respectively (Fig. 6F).
Identification of LAG3 as the most critical CD8+ T cell-associated gene in KIRC
The correlation of CD8+ T cell-related genes (GZMK, CD27, CCL4L2, FXYD2, LAG3, RGS1, CST7, DUSP4, CD8A, and TRBV20-1) was further explored with prognosis and diagnosis in KIRC. We further validated the OS, disease-specific survival, progression-free interval, diagnostic receiver operating characteristic (ROC), and time-dependent ROC curves of the above 10 genes in the TCGA–KIRC cohort. The results suggested that the area under the ROC curve was > 0.85 for all genes, except FXYD2, and the diagnostic effect was good (Additional file 4: Figure S2A–J). LAG3 exerted a significant pro-carcinogenic effect and KIRC patients with high LAG3 expression had shorter OS (p = 0.008), disease-specific survival (p = 0.006), and progression-free interval (p = 0.043) (Additional file 4: Figure S2H). In contrast, FXYD2 and DUSP4 displayed a significant anti-carcinogenic effect with all p < 0.05 (Additional file 4: Figure S2E–F), whereas the other single genes exerted no great prognostic significance (Additional file 4: Figure S2). Notably, only LAG3 exhibited a relatively good time-dependent curve, AUC1-year = 0.592, AUC3-year = 0.574, and AUC5-year = 0.580 (Additional file 4: Figure S2H). Further analysis of its association with CD8+ T cells revealed that the expression of LAG3 and DUSP4 was positively correlated with CD8+ T cell enrichment (Fig. 7A–B), whereas FXYD2 was not (Fig. 7C). The correlation between key genes and immunotherapy was further explored by TMB scores, suggesting that LAG3 was positively correlated with TMB scores (ρSpearman = 0.14 and p = 0.006), whereas DUSP4 was insignificant (Fig. 7D, F). The study revealed that the expression of LAG3 and DUSP4 was not associated with TIDE scores (Fig. 7E, G). Compared to the expression of DUSP4, LAG3 displayed a higher correlation with CD8+ T cells and a function in immunotherapy. Therefore, we followed up with an in-depth exploration of LAG3. We further verified the correlation between LAG3 expression and different T-cell subtypes in the TCGA–KIRC disease cohort. We found that the LAG3 level was positively correlated with several T cells, including CD8 T cells, cytotoxic cells, T cells, T helper cells, Tcm, Tem, TFH, Th1 cells, Th2 cells, and Tregs, and negatively correlated with Th17 cells, with no correlation with Tgd (Additional file 5: Figure S3). Therefore, LAG3 was identified as the most critical T cell-related prognostic gene in KIRC. Further investigation of the association of LAG3 expression with sunitinib and sorafenib sensitivity indicated that the high LAG3 expression group (G2) had lower sunitinib IC50 and higher sorafenib IC50 (Fig. 7H–I). Moreover, LAG3 was correlated with multiple immune-related genes (Fig. 7K). This study revealed that the LAG3 expression was not associated with stemness scores between the high and low LAG3 expression groups (Fig. 7J).
The single-cell analysis revealed that LAG3 expression was mostly concentrated in zone 1 with scattered distribution in zones 2, 5, 6, 7, and 10 (Fig. 8A). Cluster identification demonstrated that LAG3 was concentrated in CD8+ T cells in KIRC, whereas it was highly scattered in NKT or B cells in both normal kidney tissues and KIRC (Fig. 8B). The GSEA pathway after single-gene differential analysis suggested that its pathway included the PID_CD8_TCR_PATHWAY, WP_INTERACTIONS_BETWEEN_IMMUNE_CELLS_AND_MICRORNAS_IN_TUMOR_MICROENVIRONMENT, WP_CANCER_IMMUNOTHERAPY_BY_PD1_BLOCKADE and REACTOME_IMMUNOREGULATORY_INTERACTIONS_BETWEEN_A_LYMPHOID_AND_A_NON_LYMPHOID_CELL, whose NES > 1.9, p.adj < 0.05, and FDR < 0.05 (Fig. 8C–F). In addition, the combined enrichment scores of each sample on multiple pathways revealed the most significant positive correlation between LAG3 and Tumor_Inflammation_Signature, apoptosis, and IL-10_Anti-inflammatory_Signaling_Pathway (Fig. 8G–I).
Validation of LAG3 as a key CD8+ T cell-associated gene in KIRC in in vitro experiments
An analysis of the CCLE dataset revealed a significantly high expression of LAG3 in A-704, OS-RC-2, SNU-349, KMRC-3, SLR 26, and ACHN kidney cancer cell lines (Additional file 6: Figure S4). Next, the RT-PCR results of existing kidney cancer cell lines in the laboratory found LAG3 to be highly expressed in 786-O and ACHN cell lines, which were selected for subsequent validation (Fig. 9A). The knockdown efficiency of LAG3 small interfering reagents in 786-O and ACHN cell lines by RT-PCR was validated (Fig. 9B). Subsequently, the growth curve of the siLAG3 group was significantly slower in the 786-O and ACHN cell lines than that of the control group in the CCK8 experiment (Fig. 9F). Compared with the control group, the number of cell clones in the siLAG3 group was significantly lower in 786-O and ACHN cell lines by clone experiments (Fig. 9C). The scratch experiment revealed that the siLAG3 group was less efficient at healing scratches at 24 h than the control group (Fig. 9D). Cell migration and invasion were detected to be significantly down-regulated in the siLAG3 group by transwell assay (Fig. 9E). Combined with the preliminary results of LAG3 differential expression in KIRC versus normal kidney tissues in the HPA database (Additional file 3: Figure S1), immunohistochemistry revealed the high expression of LAG3 in KIRC was further confirmed in the tissues of our patients (Fig. 10).
Discussion
Kidney cancer constitutes the 14th most common malignancy worldwide with 431,288 new diagnosed cases and 179,368 new deaths in 2020 [3]. The major etiologies of kidney cancer include hypertension, obesity, and smoking [22]. The ICI-based combination therapy demonstrated excellent clinical efficacy in several large clinical trials and is now the first-line care standard for patients with advanced or metastatic renal cancer with a low OS at first diagnosis [23,24,25,26]. Despite the pivotal function of the PD-1/CTLA-4 axis in the treatment of RCC has greatly improved clinical outcomes compared to previous treatment options, the majority of patients with RCC did not achieve durable clinical benefits after ICI-based combination therapy [23,24,25,26]. Therefore, it is highly essential to investigate the tumor immune microenvironment and explore novel immunotherapeutic targets, and ultimately optimize systemic treatment for RCC.
In response to the immune microenvironment of renal cancer, Finke [6] and Stein [15] provided reviews on the immunology and immunobiology of renal cancer, respectively. In most solid tumors, the degree of CD8+ T cell infiltration was positively correlated with good prognosis for tumor patients [27]. CD8+ T cells exerted a direct cytotoxic effect on target cells and performed a critical role in anti-cancer immunity. However, the expression of suppressive molecules (PD-1 and CTLA-4) on the surface of CD8+ T cells increases in response to sustained stimulation by tumor-specific antigens, and their function decreases and eventually reaches the exhausted state, as demonstrated in multiple cancer models [13, 14, 28]. Therefore, blocking their inhibitory signaling using anti-PD1/CTLA-4 antibodies could rejuvenate exhausted CD8+ T cells, enhance their cytotoxicity, promote tumor cell lysis, and restrict tumor metastasis [29]. However, kidney cancer had a distinctive immune profile such that the degree of CD8+ T cell infiltration is positively correlated with poor prognosis, and the specific mechanism was unclear. Several hypotheses have been proposed to explain this paradoxical phenomenon. First, the activation status and virulence potential of CD8+ T cells were highly specific in kidney cancer, where stem-like TCF1+ or PD-1+ TIM3− LAG3− CD8 + T cell subsets contribute to the anti-cancer immune effect [30,31,32]. Second, the low density of tertiary lymphoid structures generated numerous immature DC cells, causing the infiltration of polyclonal CD8+ T cells that could not recognize tumor-associated antigens [11, 27, 32]. Third, the absence of tumor-specific genes, such as the relative absence of PBRM1 mutations in highly CD8+ T cells RCC, which was often associated with a better prognosis [33]. Finally, specific metabolic dysregulation of CD8+ T cells in RCC restricted CD8+ T cell activation and did not recover through the PD-1 axis inhibition [34].
The high significance of CD8+ T cells in the immunotherapy of kidney cancer and the results of single-cell sequencing analysis (CD8 + T cells were significantly differentially expressed in KIRC and kidney tissue) can be used to construct a prognostic model of CD8+ T cell-associated genes to guide clinical decisions in KIRC. First, the top 10 DEGs in CD8+ T cells of KIRC were obtained by cluster identification in single-cell analysis, including GZMK, CD27, CCL4L2, FXYD2, LAG3, RGS1, CST7, DUSP4, CD8A, and TRBV20-1. Subsequently, the above gene expression was divided into two subtypes by cluster analysis, that is, CD8+ T cell-associated gene low expression (C1) and high expression (C2) subtypes. The cluster analysis results were used for grouping, and the DEGs, pathway enrichment, and mutated genes between C1 and C2 subtypes were comprehensively studied. The up-regulated genes were largely enriched for glutamate and leukotriene activity, whereas the down-regulated genes were enriched for immune-related activities. Furthermore, VHL and STED2, the most commonly mutated genes in primary KIRC, displayed remarkably high mutation rates in the C2 subtype compared to the C1 subtype, predicting a poor prognosis such as metastasis or drug resistance in the C2 subtype [35, 36]. We assessed the KIRC immune microenvironment in different subtypes and revealed that the stromal and immune scores were significantly higher in the CD8+ T cell-associated gene high expression subtype than in the C1 subtype, with significant differences between the two subtypes in immune infiltrating cells and immune-related molecules. Therefore, LASSO regression and Cox univariate analyses were used to construct the CD8+ T cell-associated risk prognostic model: RiskScore = − 0.291858656434841*GZMK − 0.192758342489394*FXYD2 + 0.625023643446193*LAG3 + 0.161324477181591*RGS1 − 0.380169045328895*DUSP4 − 0.107221347575037*TRBV20-1, which will assist the clinicians to assess the prognosis for survival, immune status, and drug selection in patients with KIRC. Relevant studies were investigated to identify CD8+ T cell-related genes in kidney cancer. Genomics, radiology, and artificial intelligence modalities can be used to identify renal cancer differentiation more easily and earlier, and predict its Fuhrman grade and responsiveness to immunotherapy response, thus assisting the clinicians in defining the risk of stratification of the disease, treatment choices, follow-up strategies, and prognosis [42, 43].
Using tumor single-gene related studies as the reference [44], we performed single-gene analysis of the above genes and finally identified LAG3 as the most valuable CD8+ T cell-associated gene in KIRC. Lymphocyte activating gene 3 (LAG3) or CD223 is highly expressed in different T cells, CD8+ T cells, CD4+ T cells, and Tregs, to maintain homeostasis [37]. However, persistent tumor-associated antigen stimulation causes its chronic expression, ultimately promoting T-cell exhaustion in cancers [37, 38]. Therefore, LAG3 serves as the third clinical checkpoint in case of limited or even no response in 60 to 80% of cancer patients in PD1/CTLA-4 axis immunotherapy [39]. Currently, several clinical trials are being conducted on immunotherapies targeting LAG3 in combination with PD1/CTLA-4 axis inhibitors to treat cancers including kidney cancer 40, 41. In the current study, we analyzed for the first time the distribution of LAG3 in renal cancer and renal tissue using single-cell analysis and investigated its expression, prognostic significance, immune microenvironmental relevance, and pathway enrichment using bioinformatics techniques. We confirmed that LAG3 promotes the progression and metastasis of renal cell carcinoma and is positively correlated with CD8+ T cells using cell phenotype studies and immunohistochemistry. Presently, immunotherapy targeting LAG3 is largely used for melanoma, pancreatic cancer, and hematological tumors, with only a few studies on renal cancer. We elucidated the function of LAG3 in KIRC. We believe our findings will provide a preliminary basis and direction for LAG3-targeted immunotherapy and even CAR-T therapy in patients with kidney cancer.
The present study had certain limitations. First, heterogeneity obtained in retrospective studies needs to be verified by conducting prospective studies. Second, we only applied the top 10 CD8+ T cell-associated genes to construct the risk prognosis model, which lacked comprehensiveness and extensiveness. Third, we only validated the function of LAG3 in KIRC at the in vitro level and lacked in vivo experiments, as well as the expression mechanism and CD8+ T cell activity warrant further exploration. In conclusion, more basic and large clinical trials are required to validate these findings (Table 1).
Conclusion
The top 10 CD8+ T cell-associated genes were obtained by single-cell analysis, including GZMK, CD27, CCL4L2, FXYD2, LAG3, RGS1, CST7, DUSP4, CD8A, and TRBV20-1. These genes were divided into low- and high-expression subtypes by cluster analysis, and DEGs, pathway enrichment, mutant genes, and KIRC immune infiltration in different subtypes were studied. The best risk prognosis model was constructed (RiskScore = − 0.291858656434841*GZMK − 0.192758342489394*FXYD2 + 0.625023643446193*LAG3 + 0.161324477181591*RGS1 − 0.380169045328895*DUSP4 − 0.107221347575037*TRBV20-1). Finally, the single-gene analysis identified LAG3 as the most valuable CD8+ T cell-associated gene in KIRC, which was further confirmed by cell phenotype studies and immunohistochemistry.
Availability of data and materials
Renal epithelial cells (HK-2) and renal cancer cell lines (ACHN and 786-O) were obtained from the surgical laboratory of Zhongda Hospital, Southeast University, and were provided by Dr. Weipu Mao.
References
Ricketts CJ, et al. The cancer genome atlas comprehensive molecular characterization of renal cell carcinoma. Cell Rep. 2018;23:3698. https://doi.org/10.1016/j.celrep.2018.06.032.
Moch H, Cubilla AL, Humphrey PA, Reuter VE, Ulbright TM. The 2016 WHO classification of tumours of the urinary system and male genital organs—part a: renal, penile, and testicular tumours. Eur Urol. 2016;70:93–105. https://doi.org/10.1016/j.eururo.2016.02.029.
Sung H, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2021;71:209–49. https://doi.org/10.3322/caac.21660.
Ljungberg B, et al. The epidemiology of renal cell carcinoma. Eur Urol. 2011;60:615–21. https://doi.org/10.1016/j.eururo.2011.06.049.
Athar U, Gentile TC. Treatment options for metastatic renal cell carcinoma: a review. Can J Urol. 2008;15:3954–66.
Diaz-Montero CM, Rini BI, Finke JH. The immunology of renal cell carcinoma. Nat Rev Nephrol. 2020;16:721–35. https://doi.org/10.1038/s41581-020-0316-3.
Coppin C, et al. Immunotherapy for advanced renal cell cancer. Cochrane Database Syst Rev. 2005. https://doi.org/10.1002/14651858.CD001425.pub2.
Floros T, Tarhini AA. Anticancer cytokines: biology and clinical effects of interferon-alpha2, interleukin (IL)-2, IL-15, IL-21, and IL-12. Semin Oncol. 2015;42:539–48. https://doi.org/10.1053/j.seminoncol.2015.05.015.
Topalian SL, et al. Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. N Engl J Med. 2012;366:2443–54. https://doi.org/10.1056/NEJMoa1200690.
Nakano O, et al. Proliferative activity of intratumoral CD8(+) T-lymphocytes as a prognostic factor in human renal cell carcinoma: clinicopathologic demonstration of antitumor immunity. Cancer Res. 2001;61:5132–6.
Giraldo NA, et al. Orchestration and prognostic significance of immune checkpoints in the microenvironment of primary and metastatic renal cell cancer. Clin Cancer Res. 2015;21:3031–40. https://doi.org/10.1158/1078-0432.CCR-14-2926.
Jhunjhunwala S, Hammer C, Delamarre L. Antigen presentation in cancer: insights into tumour immunogenicity and immune evasion. Nat Rev Cancer. 2021;21:298–312. https://doi.org/10.1038/s41568-021-00339-z.
Farhood B, Najafi M, Mortezaee K. CD8(+) cytotoxic T lymphocytes in cancer immunotherapy: a review. J Cell Physiol. 2019;234:8509–21. https://doi.org/10.1002/jcp.27782.
van der Leun AM, Thommen DS, Schumacher TN. CD8(+) T cell states in human cancer: insights from single-cell analysis. Nat Rev Cancer. 2020;20:218–32. https://doi.org/10.1038/s41568-019-0235-4.
Drake CG, Stein MN. The immunobiology of kidney cancer. J Clin Oncol. 2018. https://doi.org/10.1200/JCO.2018.79.2648.
Verma V, et al. PD-1 blockade in subprimed CD8 cells induces dysfunctional PD-1(+)CD38(hi) cells and anti-PD-1 resistance. Nat Immunol. 2019;20:1231–43. https://doi.org/10.1038/s41590-019-0441-y.
Pauken KE, et al. The PD-1 pathway regulates development and function of memory CD8(+) T cells following respiratory viral infection. Cell Rep. 2020;31: 107827. https://doi.org/10.1016/j.celrep.2020.107827.
Borcherding N, et al. Mapping the immune environment in clear cell renal carcinoma by single-cell genomics. Commun Biol. 2021;4:122. https://doi.org/10.1038/s42003-020-01625-6.
Liao J, et al. Single-cell RNA sequencing of human kidney. Sci Data. 2020;7:4. https://doi.org/10.1038/s41597-019-0351-8.
Wilkerson MD, Hayes DN. Consensus cluster plus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26:1572–3. https://doi.org/10.1093/bioinformatics/btq170.
Danaher P, et al. Gene expression markers of tumor infiltrating leukocytes. J Immunother Cancer. 2017;5:18. https://doi.org/10.1186/s40425-017-0215-8.
Bukavina L, et al. Epidemiology of renal cell carcinoma: 2022 update. Eur Urol. 2022;82:529–42. https://doi.org/10.1016/j.eururo.2022.08.019.
Motzer RJ, et al. Avelumab plus axitinib versus sunitinib for advanced renal-cell carcinoma. N Engl J Med. 2019;380:1103–15. https://doi.org/10.1056/NEJMoa1816047.
Rini BI, et al. Pembrolizumab plus axitinib versus sunitinib for advanced renal-cell carcinoma. N Engl J Med. 2019;380:1116–27. https://doi.org/10.1056/NEJMoa1816714.
Motzer RJ, et al. Nivolumab plus ipilimumab versus sunitinib in advanced renal-cell carcinoma. N Engl J Med. 2018;378:1277–90. https://doi.org/10.1056/NEJMoa1712126.
Rini BI, et al. Atezolizumab plus bevacizumab versus sunitinib in patients with previously untreated metastatic renal cell carcinoma (IMmotion151): a multicentre, open-label, phase 3, randomised controlled trial. Lancet. 2019;393:2404–15. https://doi.org/10.1016/S0140-6736(19)30723-8.
Fridman WH, Zitvogel L, Sautes-Fridman C, Kroemer G. The immune contexture in cancer prognosis and treatment. Nat Rev Clin Oncol. 2017;14:717–34. https://doi.org/10.1038/nrclinonc.2017.101.
Zajac AJ, et al. Viral immune evasion due to persistence of activated T cells without effector function. J Exp Med. 1998;188:2205–13. https://doi.org/10.1084/jem.188.12.2205.
Huang Y, Jia A, Wang Y, Liu G. CD8(+) T cell exhaustion in anti-tumor immunity: the new insights for cancer immunotherapy. Immunology. 2022. https://doi.org/10.1111/imm.13588.
Pignon JC, et al. irRECIST for the evaluation of candidate biomarkers of response to nivolumab in metastatic clear cell renal cell carcinoma: analysis of a phase II prospective clinical trial. Clin Cancer Res. 2019;25:2174–84. https://doi.org/10.1158/1078-0432.CCR-18-3206.
Jansen CS, et al. An intra-tumoral niche maintains and differentiates stem-like CD8 T cells. Nature. 2019;576:465–70. https://doi.org/10.1038/s41586-019-1836-5.
Giraldo NA, et al. Tumor-infiltrating and peripheral blood T-cell immunophenotypes predict early relapse in localized clear cell renal cell carcinoma. Clin Cancer Res. 2017;23:4416–28. https://doi.org/10.1158/1078-0432.CCR-16-2848.
Miao D, et al. Genomic correlates of response to immune checkpoint therapies in clear cell renal cell carcinoma. Science. 2018;359:801–6. https://doi.org/10.1126/science.aan5951.
Beckermann K, et al. Targeting metabolic dysregulation of T cells in kidney cancer. J Clinical Oncol. 2020. https://doi.org/10.1200/JCO.2020.38.6_suppl.722.
Jonasch E, Walker CL, Rathmell WK. Clear cell renal cell carcinoma ontogeny and mechanisms of lethality. Nat Rev Nephrol. 2021;17:245–61. https://doi.org/10.1038/s41581-020-00359-2.
Cancer Genome Atlas Research N. Comprehensive molecular characterization of clear cell renal cell carcinoma. Nature. 2013;499:43–9.
Goldberg MV, Drake CG. LAG-3 in cancer immunotherapy. Curr Top Microbiol Immunol. 2011;344:269–78. https://doi.org/10.1007/82_2010_114.
Wherry EJ. T cell exhaustion. Nat Immunol. 2011;12:492–9. https://doi.org/10.1038/ni.2035.
Postow MA, Callahan MK, Wolchok JD. Immune checkpoint blockade in cancer therapy. J Clin Oncol. 2015;33:1974–82. https://doi.org/10.1200/JCO.2014.59.4358.
Andrews LP, Marciscano AE, Drake CG, Vignali DA. LAG3 (CD223) as a cancer immunotherapy target. Immunol Rev. 2017;276:80–96. https://doi.org/10.1111/imr.12519.
Brignone C, Escudier B, Grygar C, Marcu M, Triebel F. A phase I pharmacokinetic and biological correlative study of IMP321, a novel MHC class II agonist, in patients with advanced renal cell carcinoma. Clin Cancer Res. 2009;15:6225–31. https://doi.org/10.1158/1078-0432.CCR-09-0068.
Ferro M, et al. Artificial intelligence and radiomics in evaluation of kidney lesions: a comprehensive literature review. Ther Adv Urol. 2023;15:17562872231164804. https://doi.org/10.1177/17562872231164803.
Ferro M, et al. Radiogenomics in renal cancer managemen—current evidence and future prospects. Int J Mol Sci. 2023. https://doi.org/10.3390/ijms24054615.
Ascione CM, et al. Role of FGFR3 in bladder cancer: treatment landscape and future challenges. Cancer Treat Rev. 2023;115: 102530. https://doi.org/10.1016/j.ctrv.2023.102530.
Acknowledgements
No.
Funding
No.
Author information
Authors and Affiliations
Contributions
HF-G contributed to this manuscript and should be considered as the first authors. He was jointly responsible for the study design, data extraction and analysis as well as the writing and revision of the manuscript. WP-M and JK-Q were considered co-responding authors, who were responsible for the design of the study idea and review of the manuscript. H-S, AF-H, and H-L were considered as co-second authors who were responsible for data collection and analysis, as well as article polishing and revision. ZH-Z and DL-L were considered as co-third authors, responsible for staining and identification of pathological sections of renal carcinoma.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
All subjects gave their informed consent for inclusion before they participated in the study. The study was conducted in accordance with the Declaration of Helsinki.
Consent to publish
All the authors have read and agreed to submit the manuscript to European Journal of Medical Research with equal responsibility.
Competing interests
The authors declare that the research was conducted in the absence of any commercial or financial relationships.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 1: Table S1.
Annotated table of single cell subtypes.
Additional file 2:
Table S2. Basic information on immunohistochemistry patients in the HPA database.
Additional file 3: Figure S1.
Immunohistochemistry of CD8+ T cell-associated genes in normal renal tissue and KIRC tissue in the HPA database. (A) CD8A, (B) CD27, (C) CST7, (D) GZMK, (E) LAG3, (F) RGS1, and (G) FXYD2.
Additional file 4: Figure S2.
Overall survival, disease-specific survival, progression-free interval, diagnostic ROC, and time-dependent ROC analysis of CD8+ T cell-associated genes in the TCGA–KIRC cohort. (A) CCL4L2, (B) CD8A, (C) CD27, (D) CST7, (E) DUSP4, (F) FXYD2, (G) GZMK, (H) LAG3, (I) RGS1, and (J) TRBV20-1.
Additional file 5: Figure S3.
Correlation of LAG3 expression with different T cell subtypes. (A) CD8 T cells. (B) Cytotoxic cells. (C) T cells. (D) T helper cells. (E) Tcm. (F) Tem. (G) TFH. (H) Tgd. (I) Th1 cells. (J) Th17 cells. (K) Th2 cells. (L) TRegs.
Additional file 6:
Figure S4. Expression of LAG3 in renal cancer cell lines predicted using the CCLE database.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Gao, H., Sun, H., He, A. et al. Single-cell combined bioinformatics analysis: construction of immune cluster and risk prognostic model in kidney renal clear cells based on CD8+ T cell-associated genes. Eur J Med Res 29, 89 (2024). https://doi.org/10.1186/s40001-024-01689-8
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s40001-024-01689-8