Single-cell landscape and spatial transcriptomic analysis reveals macrophage infiltration and glycolytic metabolism in kidney renal clear cell carcinoma

Analysis of gene expression differences and prognostic significance in KIRC

In this study on differential expression of KIRC genes and their prognostic significance, various databases and tools were used to analyze gene expression levels in KIRC tumors and adjacent normal tissues, and survival analysis was performed as described in a previous study [13, 14]. Briefly, TNMplot, UALCAN (The University of ALabama at Birmingham CANcer data analysis Portal), and Gene Expression Omnibus (GEO) were accessed to retrieve relevant data. The Human Protein Atlas was also consulted to obtain immunohistochemical (IHC) staining images for further analysis. In conducting survival analysis, the Kaplan-Meier plotter tool was utilized to assess the relationship between clinical stages of KIRC and various factors such as immune cell content and tumor mutational burdens. To determine optimal patient groupings for survival analysis, the tool was configured to use the “Auto select best cut off” feature, which selects the best-performing cutoff value. Overall, the various analyses performed in this study provide important insights into the gene expression and prognostic significance of KIRC, and may inform future research and clinical decision-making in this area.

Single-cell transcriptomic and immune profiling analysis

Single-cell RNA sequencing (scRNA-seq) has emerged as a powerful tool to explore cellular heterogeneity and gene expression at the single-cell level. In this study, scRNA-seq data were obtained from the GSE159115 and GSE121636 datasets of the GEO database. Quality control (QC) was performed using the R package Seurat to ensure the inclusion of high-quality cells and reduce the influence of batch effects. To identify distinct cell subpopulations, uniform manifold approximation and projection (UMAP) clustering was performed using the “BiocManager” and Gene Set Variation Analysis (“GSVA”) packages in R. Cell type annotation was performed by comparing the expression profiles of previously recognized cellular marker genes with those of the identified cell subpopulations using the “SingleR” package in R. The association between gene expression and immune cell infiltration/abundance was explored using the “Gene” module of TIMER, focusing on macrophages, dendritic cells (DC), and CD4+ T, CD8+ T, and B cells. Finally, the association between immune infiltration and overall survival of KIRC patients was estimated using several algorithms, including TIMER, EPIC, MCPCOUNTER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, and XCELL.

Processing of spatial transcriptomics data using seurat algorithm

We performed analysis based on previously published articles and established spatial transcriptome data of kidney injury model animals [15]. The spatial transcriptomics data were processed using the R package Seurat (v4.1.0) and log normalization was performed to normalize the data. To integrate Seurat objects into a single ST dataset and remove batch effects, we utilized the functions SelectIntegrationFeatures, FindIntegrationAnchors, and IntegrateData. To reduce dimensionality, we performed the function RunPCA, and then used the functions FindNeighbors and FindClusters to cluster similar ST points. Initially, different clusters were divided based on hematoxylin-eosin stained (H&E) sections and annotated by unsupervised cluster analysis. Upon annotating clusters with cell markers, we observed that some clusters highly expressed multiple cell markers. To address this, we employed the ssGSEA algorithm to score common cell types based on the average expression matrix of different clusters, which was confirmed to be more effective in ST [16].

Analysis of human KIRC specimens

Tissue microarray (TMA) slides (CL2) consisting of human renal cancer, metastatic, and normal tissues were procured from SuperBioChips Laboratories (Seoul, Republic of Korea). To perform an immunohistochemistry (IHC), the protocol described in a prior study [17] was followed.

Colocalization analysis of immunofluorescence staining

We employed a method that has been described previously, which is summarized as follows [18]. Immunofluorescence analysis was performed by fixing the tissue with 4% paraformaldehyde at room temperature for 10 minutes, permeabilizing it with 0.2% Triton X-100 for 5 minutes, blocking it with 5% BSA for 30 minutes, and incubating it overnight at 4° C with the specified primary antibodies PGAM1_A4170 and CD163_A22619 purchased from Abclonal (Woburn, MA, USA). After washing with PBST, the samples were incubated with secondary antibodies for 30 minutes at room temperature. The BX61VS® Fully Motorized Fluorescence Microscope (Olympus Corporation, Tokyo, Japan) at ×20 magnification was used to digitize all glass slides with high precision. The whole-slide images were viewed and analyzed with Olympus VS-ASW® software at Li-Tzung Pathology Laboratory (Kaohsiung, Taiwan).

Pharmacogenetic prediction model

We developed a pharmacogenetic prediction model to analyze drug sensitivity based on PGAM1 expression, using the knockdown-screen data repository of the Genomics of Drug Sensitivity in Cancer (GDSC) algorithm in Q-omics v.1.0. Pearson correlation coefficient analysis was performed to assess the correlation between PGAM1 expression levels and drug dose levels [19].

Statistical analysis

The statistical analysis in this study followed the methods described in a previous publication [20]. Pearson’s correlation coefficient was used to evaluate gene expression correlation, while a t-test or Fisher’s exact test was used for the comparison between two groups, and one-way ANOVA was used for comparison within one group, all analyzed using GraphPad Prism software (GraphPad Software, La Jolla, CA, USA). A p-value of less than 0.05 was considered statistically significant.

Differential expression and genetic variation analysis of PGAM1 in KIRC

We first conducted an analysis of the expression levels of PGAM family genes in the male reproductive system and observed that PGAM1 exhibited higher expression in male reproductive organs of fat fruits (Figure 1A). We subsequently explored the association of PGAM1 mutations with common cancer progression genes, such as VHL, PBRM1, and SETD2 (Figure 1B). Subsequently, we conducted an analysis in which we extracted data pertaining to PGAM1 expression from the TCGA database, focusing on patients diagnosed with Kidney Renal Clear Cell Carcinoma (KIRC). This information was presented using a waterfall plot, offering a visual representation of the top 25 genes that exhibited notable alterations in response to PGAM1 expression changes. This analysis delved into the genetic variations that occur at different levels of PGAM1 expression and their connections with genes that are frequently found to be mutated in KIRC. The insights gained from this examination are depicted graphically in Figure 1C, providing a comprehensive view of these associations and their potential implications in the context of KIRC. The mutation landscape identified various types of mutations, including splice site, missense, frameshift, and nonsense mutations, as well as in-frame ins/dels. Moreover, we used protein-protein interactions to evaluate regulatory network maps and inferred underlying mechanisms. Our analysis revealed that PGAM1 was associated with VHL, which was consistent with the results from Figure 1C. Furthermore, we found that the intermediate factor connecting these two genes was HIF1A (Figure 1D). We conducted univariate and multivariate Cox regression analyses to evaluate the association of PGAM1 expression with OS. Our findings indicated that stage, sex, age, and PGAM1 (high vs low) were significantly associated with OS (Figure 1E). Finally, we applied multivariate Cox regression analysis to the same variables and found that the risk score could serve as an independent prognostic factor (P < 0.05) (Figure 1F).

Subsequently, we conducted a comprehensive assessment of various kidney cancers and discovered that the levels of phosphoglycerate mutase-1 (PGAM1) were elevated in both clear cell renal cell carcinoma (KIRC) and papillary renal cell carcinoma (KIRP) but reduced in chromophobe renal cell carcinoma (KICH) (Figure 2A). Upon scrutinizing the expression of PGAM1 across diverse sample types, cancer stages, metastasis stages, and ccRCC subtypes in The Cancer Genome Atlas (TCGA) (Figure 2B–2E), we observed that PGAM1 expression was significantly amplified in the first stage of metastatic KIRC, and both subtypes of ccRCC exhibited an upward trend. Furthermore, analysis of immunohistochemical staining data from the Human Protein Atlas verified that the PGAM1 protein was upregulated in KIRC tumor tissues relative to normal tissues (Figure 2F). Strikingly, Kaplan-Meier analysis revealed that patients with elevated PGAM1 immunohistochemical scores had a shorter overall survival time (Figure 2G), underscoring the potential clinical significance of PGAM1 in ccRCC.

Spatial transcriptomics evaluation of PGAM1 expression in injured kidneys

To assess the changes and distribution of PGAM1 expression in renal injury, we conducted an analysis of previously published data [15] that utilized two models of injured kidneys for spatial transcriptomic sequencing. The intact kidney tissue was sectioned coronally to expose each major physiologic region, cryosectioned and placed on a tissue capture area on a dedicated 10× Genomics Visium slide embedded with oligonucleotide sequences. Marker gene expression was visualized using Seurat, and distinct regions were identified by biomarkers in proximal tubules (Lrp2) and collecting ducts (Aqp2) (Figure 3A). Our results showed that PGAM1 expression was higher in the proximal tubules of renal injury compared to the control group. Furthermore, we also analyzed the gene levels of common kidney injury biomarkers as controls, including lipocalin/lipocalin 2 (Lcn2), kidney injury molecule 1 (Havcr1), insulin-like growth factor binding protein 7 (Igfbp7), all of which have significantly high expression in kidney injury. Therefore, we concluded that the metabolic gene PGAM1 was highly expressed in the damaged kidney, which may be related to the damage of this region and the change of metabolic pathway.

Correlation analysis between PGAM1 and infiltrating immune cells

In light of the role of the tumor immune microenvironment (TIME) in promoting kidney damage and recruiting immune cells, we further investigated the relationship between metabolism and the immune response in tumors. Our analysis showed that overexpression of PGAM1 resulted in Branch amino acid degradation, as reflected in the KEGG pathway analysis (Figure 4A). Moreover, the activation of the Immunoglobulin complex was observed in the GO-cellular component analysis (Figure 4B), and the same result was reflected in the GO-molecular function analysis, where the ratio of Immunoglobulin receptor binding and Immune response-mediated biological function was greatly increased (Figure 4C, 4D). Based on these findings, we hypothesize that kidney damage, which requires high energy, not only leads to metabolic conversion and abnormality but also promotes infiltration of immune cells, playing a significant role in the tumor immune microenvironment. The activation of the Immunoglobulin complex and increased Immunoglobulin receptor binding observed in our analysis suggests that immune cells may be recruited to the damaged kidney, contributing to the immune response in tumors.

Evaluation of PGAM1 transcriptome using single-cell RNA sequencing database analysis

To assess the transcriptome of PGAM1 in renal cancer at the single-cell level and to explore the heterogeneity of different cell types in the renal cancer microenvironment, we conducted an analysis of two publicly available single-cell RNA-sequencing databases (GSE159115 and GSE121636). A total of 8 samples from GSE159115, comprising 8 types of cells and 32 clusters, and 3 samples from GSE121636, comprising 10 cells and 21 clusters, were included in the study. Using UMAP plots, we identified and visualized 8 and 10 major cell populations in GSE159115 and GSE121636, respectively (Figure 5A, 5B). After quality control and batch effect removal, we analyzed 27,669 (GSE159115) and 25,681 (GSE121636) cells, respectively. We determined cell type-specific markers for each cluster based on the top-ranked differentially expressed genes, which were used for cell type classification. Figure 5C, 5D shows the clinical information for each cell population, where differences in the corresponding proportions of each cluster were observed for different clinical characteristics. We further analyzed the expression and distribution of PGAM1 in these two single-cell sequencing databases. As shown in Figure 6A, 6B, we found that PGAM1 expression levels in different cell types of renal cancer were high in the same region as the high glycolysis model. In addition, oxidative phosphorylation in Figure 6B was generally increased, which was also reflected in the Inflammatory response.

PGAM1 over-expression is correlated with immune infiltration

To investigate the potential relationship between KIRC and PGAM1, we performed gene set enrichment analysis (GSEA) using differentially expressed genes (DEGs) and identified highly correlated clusters with macrophages in two single-cell RNA sequencing databases (Figure 7A, 7B). However, macrophages exhibit tissue residency and possess pro- or anti-inflammatory functions. Given our previous findings of PGAM1 overexpression in renal injury models and human KIRC, and the crucial role of the tumor immune microenvironment (TIME) in tumor growth, metastasis, and immune evasion, we further analyzed the relationship between PGAM1 and tumor-immune interactions.

Based on the hypothesis that PGAM1 expression might be associated with the tumor microenvironment of renal cancer and differ significantly in macrophages (Figure 8A), we assessed the correlation of macrophage subtypes (M0, M1, and M2) with PGAM1. Our analysis showed a positive correlation between PGAM1 and M1 and M2, while M0 showed the opposite trend (Figure 8B). We further analyzed macrophage biomarkers and found that PGAM1 was positively correlated with CD86, CCL15, CXCL10, CD163, TLR1, and TLR8 (Figure 8C). This result indicates that patients with high PGAM1 expression exhibited elevated renal cell score, immune score, and ESTIMATE score.

To validate the connection between PGAM1 and CD163 in KIRC, we investigated the co-localization of these two proteins in TMA through triple-labeling the whole KIRC with immunofluorescence for PGAM1, CD163, and DAPI. We segmented the human KIRC TMAs into early and late stages and analyzed their fluorescence intensity through panoramic tissue scanning. In Figure 9A, the red fluorescent PGAM1 and green fluorescent CD163 demonstrated that the normal group had a higher intensity than the tumor group, and the red fluorescent intensity was higher in the early stage than in the late stage, which was similar to Figure 2F. The green fluorescence was consistent with the late stage of the normal group being lower than that of the tumor group. The 2.5D image reconstruction showed a clear overlapping phenomenon of PGAM1 and CD163 signals. After quantitative analysis by fluorescence, we confirmed that PGAM1 was positively correlated with CD163 (r = 0.3119; p = 0.0275) (Figure 9B). Taken together, these findings indicate that PGAM1 overexpression is associated with the active TIME of macrophages.

Pharmacogenetic screening for potential PGAM1 inhibitors in KIRC

In this study, we aimed to identify drugs with potential efficacy against KIRC by analyzing the GDSC repository for drugs that exhibit enhanced potency in the presence of high PGAM1 expression. We performed a cross-correlation analysis to investigate the effects of 473 drugs on shRNA-mediated PGAM1 in various KIRC cells. Our analysis identified four drugs, namely GNE-317, MS-275, AC45971100, and NSC-35468, that displayed altered potency (Figure 10A). Notably, KIRC cell lines with high shPGAM1 efficiency exhibited increased sensitivity to these drugs (Figure 10B–10E). Collectively, these findings suggest that GNE-317, MS-275, AC45971100, and NSC-35468 hold potential as anticancer agents targeting PGAM1 to regulate KIRC growth.