Integrative analysis of LAG3 immune signature and identification of a LAG3-related genes prognostic signature in kidney renal clear cell carcinoma

LAG3 is upregulated in KIRC and correlates with a poor prognostic phenotype

The TIMER database was used to analyze differential expression of LAG3 at pan-cancer level. Compared with normal tissues, LAG3 was highly overexpressed in certain cancer types including BRCA (Breast Invasive Carcinoma), ESCA (Esophageal Carcinoma), GBM (Glioblastoma), HNSC (Head and Neck Squamous Cell Carcinoma), LUAD (Lung Adenocarcinoma), LUSC (Lung Squamous Cell Carcinoma), PCPG (Pheochromocytoma and Paraganglioma), and KIRC, while down-regulation was observed in COAD (Colon Adenocarcinoma), KICH (Kidney Chromophobe), PRAD (Prostate Adenocarcinoma), LIHC (Liver Hepatocellular Carcinoma), READ (Rectum Adenocarcinoma), THCA (Thyroid Carcinoma), and UCEC (Uterine Corpus Endometrial Carcinoma) (Figure 1A). Although the assessment of LAG3 expression in KIRC was limited to 533 tumor samples and 72 adjacent normal tissue samples in the TIMER database, notable differences in LAG3 expression were observed between KIRC tumor tissues and normal tissues.

The TISIDB database was utilized to explore the association between LAG3 expression and clinical stage and grade of cancers. Notably, significant correlations between LAG3 expression and stage and grade were identified specifically in renal cancer (Figure 1B, 1C). Therefore, we focused on clarifying the functional role of LAG3 in KIRC. According to TCGA KIRC RNA-seq data, LAG3 was overexpressed in both paired and unpaired cancer tissues as compared to healthy controls (Figure 1D). To explore the clinical significance of LAG3 in KIRC, we initially investigated its expression pattern in relation to various clinical and pathological characteristics. We constructed Sankey plots to demonstrate trends in high and low LAG3 gene expression in KIRC patient samples with different clinical features and survival conditions (Figure 1E). The findings revealed elevated expression levels of LAG3 in advanced stages and higher grades KIRC (Figure 1F). Additionally, LAG3 overexpression was observed in lymph node metastasis and distant metastatic tissues when compared to normal tissues (Figure 1G). The same outcomes are demonstrated by detailed information regarding particular thresholds, namely that LAG3 high and low expression groups are statistically distinct in the T-stage, N-stage, M-stage, and AJCC-stage (Table 1).

High LAG3 expression was associated with more dead events (Figure 1H). Subsequent survival analysis to assess the prognostic value of LAG3 in KIRC revealed an association between high LAG3 and lower survival compared with that in the low expression group (P = 0.00018; Figure 1I). In contrast, there was no difference in the RFS between the high and low LAG3 expression groups (Supplementary Figure 1A). We then conducted univariate and multivariate Cox regression analysis to examine the influence of LAG3 on KIRC survival. In the univariate regression analysis, LAG3 expression, AJCC stage, tumor pathological stage (T), and the presence of metastasis (M) were all identified as significant prognostic risk factors for KIRC (Figure 1J). However, after all relevant factors were included in the multivariate regression analysis, LAG3 expression was not identified as an independent prognostic predictor of KIRC (hazard ratio (HR) = 1,054, P = 0.238; Supplementary Figure 1B).

Based on these discoveries, we applied receiver operating characteristic (ROC) curve analysis to evaluate the diagnostic accuracy of LAG3 expression in discriminating KIRC tumor tissue from normal tissue, as well as different clinicopathological characteristics. LAG3 exhibited remarkable discriminatory potential in distinguishing tumor tissue from normal tissue, as demonstrated by an impressive area under the curve (AUC) of 0.854 (95% confidence interval (CI): 0.874–0.939). These results highlight the potential of LAG3 as a promising diagnostic marker for KIRC (Figure 1K). In the comparisons of the G3&G4 versus (vs.) G&G2 groups, the T3&T4 vs. T1&T2 groups, the M1vs. M0 groups, and the N1 vs. N0 groups, the AUCs were 0.640 (95% CI: 0.592–0.689), 0.624 (95% CI: 0.577–0.671), 0.628 (95% CI: 0.558–0.697), and 0.697 (95% CI: 0.567–0.826), respectively (Supplementary Figure 1C–1F).

Taking into consideration LAG3 expression along with relevant clinicopathological characteristics, we employed the Cox regression algorithm to develop nomograms for predicting the 1-year, 3-year, and 5-year overall survival (OS) rates of KIRC patients. When compared to the ideal model, the calibration plots for the 1-, 3-, and 5-year OS rates demonstrated favorable predictive performance in the entire cohort. These results indicated the potential utility of the generated nomograms for accurate prediction of the survival outcomes of KIRC patients (Figure 1L, 1M). As an inhibitory immune checkpoint, high expression of LAG3 is associated with a positive response to immunotherapy due to its impact on the cytotoxic activity of immune cells against tumor cells. Moreover, elevated LAG3 expression may indicate a poorer response to conventional anticancer treatments, such as traditional chemotherapy. Our analysis of the role of LAG3 expression in predicting responses to common anticancer drugs supports this notion. The results demonstrate that most anti-tumor drugs, including sorafenib, exhibit higher drug sensitivity in the LAG3 low-expression group. This suggests that the expression level of LAG3 holds valuable reference significance in the selection of anti-tumor drugs (Supplementary Figure 2).

LAG3 is predominantly expressed on exhausted CD8⁺ T cells and strongly co-expressed with PDCD1

After identifying the elevated expression of LAG3 in KIRC, we further analyzed its expression in subcellular types. Through immunofluorescence experiments, we observed that LAG3 expression is observed in immune cells (CD45+ cells) rather than in tumor cells (PAX8+ cells) (Figure 2A). To precisely identify the immune cell types expressing LAG3, we probed the cell-type localization of LAG3 expression in KIRC using two single-cell datasets. In both datasets, LAG3 was found to be expressed mainly on exhausted CD8⁺ T cells (CD8Tex cells) rather than on tumor cells or other types of immune cells compared and was also co-expressed with PDCD1(Figure 2B–2D). The pie charts show that CD8Tex occupies a non-negligible proportion of the KIRC immune cell composition (Figure 2E). Using the TISIDB database, we found that the expression patterns of LAG3 and PDCD1 showed similar association with OS, clinical stage, and grade in KIRC (Figure 2F–2I).

The early exhaustion of CD8⁺ T cells has been proposed as an indicator of a favorable clinical response to immune checkpoint inhibitor (ICI) therapy. Conversely, the emergence of late-stage exhaustion markers and post-treatment depletion markers has been associated with an unfavorable clinical outcome [13]. CD8Tex cells have the potential to restore their cytotoxic functionality, reverse T cell depletion, and enhance their capacity to eliminate tumors by inhibiting immunological checkpoints present on their cell surface [14]. From this, it can be inferred that LAG-3 and PDCD1 dual blockade may be a highly clinically active immunotherapy for KIRC.

LAG3 expression is strongly positively associated with the majority of immune checkpoints and immune-related molecules

Based on the previous results, we further explored the correlation of LAG3 with other immune checkpoints and immune-related molecules in KIRC. Theoretically, immunoinhibitors tend to be highly expressed in the TME [15]. To gain a deeper understanding of the synergistic impact on immune response, we investigated potential associations between LAG3 and a diverse array of immune checkpoint inhibitors, encompassing PDCD1 (PD-1), CD274 (PD-L1), TIGIT, CTLA4, BTLA, and PDCD1LG2. The results indicated that LAG3 might modulate the antitumor immune response through co-regulation with these immune checkpoint molecules (Figure 3A). Furthermore, LAG3 expression was favorably linked with a wide variety of immunostimulators, including CD27, CD28, CD80, and CD86, indicating that LAG3 and immunostimulators interact to regulate immunological homeostasis (Figure 3B). LAG3 was also positively associated with almost all MHC-related genes, suggesting that LAG3 may interact with MHC molecules in addition to MHCII (Figure 3C). LAG3 exhibited positive associations with key receptors involved in the recruitment of effector tumor-infiltrating immune cells (TIICs), such as CVCR3-CVCR6, as well as chemokines including CCL3-CCL5 and CXCL9-CXCL13 (Figure 3D, 3E). Thus, further research into the connection between LAG3 and chemokines is warranted.

Association between LAG3 expression and the cancer-immunity cycle

Previous reports have shown that various chemokines and immunomodulators modulate the activity of the cancer-immunity cycle [16]. Therefore, we investigated a potential association of LAG3 with the cancer-immunity cycle in KIRC. We found that LAG3 expression was positively correlated with most immunity cycle activities, especially the trafficking of immune cells to tumors (step 4). LAG3 expression was strongly correlated with recruitment of CD8⁺ T cells (r = 0.74), NK cells (r = 0.82), macrophages (r = 0.70), and Th1 cells (r = 0.68), intermediately correlated with recruitment of dendritic cells (r = 0.58), and weakly correlated with recruitment of CD4⁺ T cells (r = 0.37), Th2 cells (r = 0.12), Treg cells (r = 0.12), B cells (r = 0.27). In accordance with our hypothesis, a negative correlation was observed between LAG3 expression and the efficacy of cancer cell eradication (step 7) (Figure 4), which confirmed its role as an inhibitory immune checkpoint that suppresses immune surveillance. Consistently, these results also revealed the role of LAG3 in promoting TIICs infiltration and shaping the inflammatory TME in KIRC.

LAG3 expression is correlated with immune infiltration in KIRC

Analysis of the relationship between LAG3 and six main types of TIICs using TIMER revealed that LAG3 expression correlated positively correlates with CD8+ T cells, B cells, and dendritic cells (Figure 5A). Analysis of the abundance of infiltrating immune cells based on TCGA-KIRC data revealed that high LAG3 expression tended to imply richer infiltration by immune-activated TIICs, such as CD8+ T cells, CD4+ T cells, and Treg cells (Figure 5B). Additional exploration of LAG3 expression in immune cell markers demonstrated that LAG3 was significantly correlated with CD8+ T cells, Treg cells, and Th1 cells, and only moderately correlated with M1 macrophage cells, M2 macrophage cells, and Th2 cells (Figure 5C). Subsequent investigations revealed notable distinctions between the high- and low-LAG3 expression groups. Specifically, the high-LAG3 group exhibited elevated immune scores, stromal scores, and ESTIMATE scores when compared to the low-LAG3 expression group (Figure 5D). Correlation analysis further demonstrated a robust association between LAG3 expression and immune score (p < 0.001, r = 0.724), a moderate correlation with ESTIMATE score (p < 0.001, r = 0.579), and a weaker correlation with stromal score (p < 0.001, r = 0.179) (Figure 5E). In consideration of the significance of mutational burden in immunotherapy, we analyzed the correlation between LAG3 expression and Tumor Mutational Burden (TMB) as well as Microsatellite Instability (MSI). LAG3, as an immune checkpoint marker gene, is infrequently mutated. Similar to PD-1, MSI, and TMB, it serves as an independent factor for evaluating the efficacy of immunotherapy. Therefore, the expression level of LAG3 is not directly associated with MSI or TMB. Our results confirm this, as the correlation analysis indicates a weak association between LAG3 expression and TMB (p = 0.006, r = 0.14), and no statistically significant correlation with MSI scores (p = 0.924, r = 0.01) (Figure 5F).

LAG3 interaction network construction and enrichment analysis

Using STRING, GeneMANIA, BioGRID, and HitPredict to screen and collect LAG3 co-expressed genes, we identified eight associated genes in the intersection to better characterize the function of LAG3 in KIRC (Figure 6A–6D, Supplementary Table 1). The co-expressed genes associated with LAG3 were further examined by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses, providing valuable insights into the underlying biological mechanisms (Figure 6E). The results revealed significant enrichment in key biological processes such as T-cell activation and cell-cell adhesion, indicating that LAG3 cannot only blocks T cell activation but is also responsible for cell adhesion in KIRC.

Construction of a prognostic model of LAG3-related genes

Univariate and multivariate analyses of the eight LAG3-related genes revealed that CENPJ, LAG3, and SNX18 were independent prognostic variables for KIRC patients (Figure 7A, 7B). Based on this information, we developed a predictive nomogram that demonstrated favorable performance when compared to the ideal model in terms of the rates of survival at 1-, 3-, and 5-year survival rates (Figure 7C, Supplementary Figure 3A). A predictive signature consisting of these three prognostic biomarkers (CENPJ, LAG3, and SNX18) was created using LASSO Cox regression analysis based on these potential prognostic biomarkers (Supplementary Figure 3B, 3C). Following the calculation of the risk score for KIRC patients using the equation ((0.4692× CENPJ expression) + (0.1458 × LAG3 expression) + (−0.5678 × SNX18 expression)), the KIRC patient cohort was segregated into high-risk and low-risk groups. The data analysis revealed that KIRC patients with elevated risk scores exhibited an unfavorable prognosis (Figure 7D, 7E). Furthermore, the area under the receiver operating characteristic (ROC) curves at 1, 3, and 5 years were determined as 0.671, 0.653, and 0.67, respectively, thereby illustrating the robust predictive capacity of this feature in predicting outcomes for KIRC patients (Figure 7F). Through further analysis of immune features, we observed differences in the composition of immune cells between the high-risk and low-risk groups (Supplementary Figure 3D). In comparison to the low-risk group, the high-risk group demonstrated higher levels of immune cell infiltration and association with immune pathways (Figure 7G, 7H). The GSEA analysis also revealed enrichment in immune-related pathways for the high-risk group (Figure 7I). This indicates that there is a difference in the immune status between the high-risk and low-risk groups, suggesting that the high-risk group may have potential predictive value for the immunotherapeutic responsiveness of KIRC.

Data acquisition and expression analysis

The RNA-sequencing data of KIRC obtained from The Cancer Genome Atlas (TCGA) database were subjected to normalization, tumor purity assessment, and gene annotation conversion processes. These operations yielded expression data for 610 samples (Tumor: 538, Normal: 72). Corresponding clinical data for tumor samples were also obtained from the TCGA-KIRC dataset. All data were then processed with the “TCGAbiolinks” in the R package. All gene expression values underwent a logarithmic transformation and subsequent normalization process. The R package “ggalluvial” was used to construct the Sankey diagram. The R package “ggstatsplot” was employed to construct the correlation ridge plots. We examined the cell-type localization of LAG3 and PDCD1 expression in KIRC using single-cell sequencing datasets from Tumor Immune Single-cell Hub 2 (TISCH2) [30], an online site that provides detailed cell-type annotations.

Drug sensitivity analysis

Using the R package “pRRophetic” algorithm, we calculated the half-maximal inhibitory concentration (IC50) of commonly used drugs for KIRC patients in the LAG3 high and low expression groups. The results were visualized using the R packages “ggplot2” and “ggpubr”.

Immune feature analysis

The list of immunoinhibitors, immunostimulators, MHC molecules, chemokines, and receptor molecules was obtained from the TISIDB online website, a comprehensive repository portal for tumor immune system interactions [31]. Analysis of LAG3 correlation with tumor-infiltrating immune cells (TIICs) and cell surface molecular markers was done using the TIMER database, a website that comprehensively analyzes immune infiltration in various cancers [32]. Immune infiltration in the high and low LAG3 expression groups in TCGA-KIRC RNA-seq was compared using the “ssGSEA” algorithm provided in the R package “GSVA”. The ESTIMATE algorithm [33] was used to examine the immune score, stromal score, and tumor purity, which provides valuable information about the presence and composition of immune cells and stromal cells within the tumor microenvironment (TME). Furthermore, we investigated the association between LAG3 expression and the seven-step cancer immunity cycle, which represents the intricate interplay between immune responses and the fate of tumor cells, using the Tumor Immune Profiling (TIP) approach [34]. The obtained results were visualized using the “ggpubr” and “ggplot2” in the R packages.

Protein-protein interaction network construction and correlation analysis

Genes co-expressed with LAG3 were obtained by searching the STRING, GeneMANIA, BioGRID, and HitPredict databases [35–37]. After screening the recurring genes, eight LAG3-related genes were identified. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted to investigate the function of the LAG3-dominated gene interaction network using “clusterProfiler” package in R software.

Prognostic analysis and establishment of a LAG3-related genes risk signature

We utilized the Kaplan-Meier plotter [38] to evaluate the association between differential expression of LAG3 and overall survival (OS) in KIRC patients. We also performed univariate and multivariate Cox regression analyses to assess the prognostic significance of LAG3. The results were visualized as forest plot generated using R package “survivor and survminer”. A nomogram was produced using the “rms” software and the total recurrence rate in 1-, 3-, and 5-year OS were forecasted based on the results of multivariate Cox proportional hazards analysis.

The regression coefficient for the three most statistically significant prognostic genes were obtained from a multivariate Cox proportional hazards regression model. To perform the LASSO Cox regression, we utilized the package “glmnet” in R and estimated the penalty parameter through 10-fold cross-validation. Subsequently, we calculated the risk score using the following formula:

$$\text{Riskscore}={\displaystyle \sum \text{Expi}\text{×}\text{βi}}$$

Based on the median risk score, the patients were categorized as high- and low-risk groups. We then conducted Kaplan-Meier survival analysis and employed a log-rank test to assess the statistical significance of survival outcomes between these two groups.

Functional state analysis of LAG3-related gene risk features

Gene Set Enrichment Analysis (GSEA) analysis was performed using the R package “clusterProfiler” to characterize the potential biological functions of genes between the high-risk and low-risk groups. The reference gene set used was c2.cp.all.v2022.1.Hs.symbols.gmt (All Canonical Pathways). Visualization of enrichment analysis results was conducted using the “ggplot2” package. Single-sample Gene Set Enrichment Analysis (ssGSEA) was carried out with the R package “GSVA” to assess the infiltration levels of immune cell types and immune pathways between the high-risk and low-risk groups.

Immunofluorescence (IF) staining

After obtaining the patient’s consent and approval from the Ethics Committee of Qingdao Central Hospital (KY202318401), we obtained pathological sections from KIRC patients. Anti-LAG3 antibody (ab209236, 1:100 dilution), anti-CD45 antibody (ab40763, 1:100 dilution), and anti-PAX8 antibody (ab191870, 1:1000 dilution) were purchased from Abcam (Cambridge, UK). BSA (G5001), PBS (G0002), EDTA (G1203), and DAPI (G1012) were purchased from Servicebio (Wuhan, China). Paraffin-embedded tissue sections underwent xylene deparaffinization, followed by ethanol dehydration. Microwave antigen retrieval (EDTA buffer, pH 9.0, 10 minutes) was performed, and non-specific binding was blocked with BSA for 30 minutes. Anti-CD45, anti-PAX8, and anti-LAG3 antibodies were applied and incubated overnight at 4°C. After three PBS washes, sections were exposed to secondary antibodies for 50 minutes at room temperature under light protection. DAPI staining (10 minutes) for cell nuclei and a 5-minute application of autofluorescence quenching reagent were followed. Sections were then mounted on slides for microscopy, and images were captured using a fluorescence microscope.

Statistical analysis

Differences between groups of data were assessed using the Wilcoxon test or Kruskal–Wallis test as appropriate. Pearson correlation analysis was conducted to determine correlations between variables. Kaplan–Meier analysis and log-rank tests were conducted for survival analysis. Univariate and multivariate Cox regression analyses were employed to identify independent prognostic factors. All statistical analyses were performed using R software (version 4.0.4) and P < 0.05 was considered to indicate statistical significance.