A training set of TCGA-KIRC data was downloaded from the Cancer Genome Atlas (TCGA) database, which contained gene expression data from 541 tumor tissues and 72 normal tissue samples, as well as corresponding clinical information. The E-MTAB-1980 validation set, which contains clinical and gene expression data from 101 patients with ccRCC, was generated from ArrayExpress. Our next step was to remove genes with raw counts below 10 in more than 25% of the samples. The TPM data were transformed into log2 (TPM+1). There are 1223 genes associated with lactylation modification according to Zhang et al. (12).

The 1223 lactylation modification-related genes were subjected to Pearson correlation analysis with 23 m6A genes to obtain crosstalk genes with screening criteria of correlation>0.5 and padj<0.01. To identify differentially expressed genes (DEGs) between cancer and paracancerous tissues, we used the “DESeq2” package [padj <0.05, |log2fold change (FC)|>1]. Subsequently, the crosstalk genes were merged with DEGs and the intersection was taken to finally obtain the differentially expressed crosstalk genes (DECGs). The correlation results between DECGs were visualized using the “circlize” package.

Each sample’s tumor mutation load (TMB) was calculated using the “Maftools” package. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of DECGs were performed using the “clusterProfiler” package.

We identified ccRCC clinical subtypes using the consensus clustering R package “ConsensusClusterPlus” (17). The ability of DECGs to discriminate between subtypes was assessed using principal component analysis (PCA). After that, we compared the differences between subtypes in terms of clinical variables (age, gender, grade, and stage) as well as overall survival (OS).

Every sample was analyzed for the TME score using the “estimate” package. In each TCGA-KIRC sample, immune cell infiltration was assessed using the online analysis tool TIMER2.0 (CIBERSORT algorithm). Based on our knowledge of the close relationship between immune-inhibitory, immune-stimulatory, and human leukocyte antigen (HLA) genes and TME, a comparison of expression levels between subtypes was made (18, 19).

Univariate COX regression analysis was first performed on the DECGs to screen for genes associated with OS according to p<0.05, and they were used to perform a least absolute shrinkage and selection operation (LASSO) analysis to screen for genes most associated with prognosis for constructing the risk model and to derive a risk coefficient for each gene. The expression of each modeled gene was multiplied by the risk coefficient to calculate the risk score for each patient. Each group of patients was categorized according to the median risk score. To test whether the model was able to discriminate between patients at different risks, PCA and Kaplan-Meier (K-M) survival analyses were performed. In addition, the relationship between clinical variables and risk scores for different clinical characteristics in the high- and low-risk groups was assessed.

Based on the above results, we evaluated the differences between the two risk groups in terms of TME scores, immune cell infiltration, and immune-related gene expression. Then, we calculated immune-related function scores using the “ GSEABase” and “GSVA” packages based on the “ immune. gmt” file. The DEGs were determined using the “limma” package, followed by GO, KEGG, and gene set enrichment analysis (GSEA) using “ClusterProfiler”. The “enrichplot” and “GseaVis” packages were used to visualize the enrichment analysis results.

As a result of the analysis above, we obtained the sample TMB and then downloaded the MSI score file via the “cBioPortalData” package. We defined samples as MSI when the score exceeded 0.3, and MSS if the score was below it. Moreover, the “easier” package calculates an immunotherapy response score based on TME characteristics, where higher scores indicate greater immunotherapy sensitivity (20). We downloaded each patient’s immunophenotype score (IPS) from the Cancer Immunome Atlas (TCIA, https://tcia.at/home) and divided it into <=8 and >8 groups according to IPS. The relationship between these metrics and risk scores was finally evaluated to reflect the predictive value of risk scores on immunotherapy response.

Through multivariate and univariate COX regressions, several independent predictors of OS were identified. A nomogram was then created with the help of the “survival” and “rms” packages. Calibration plots, Receiver Operating Characteristics (ROC), and Area Under Curve (AUC) were used to assess the predictive capability of the nomogram. Decision curve analysis (DCA) was used to determine the clinical value of the nomogram.

Both the metabolism of lactate and the lactylation modification process are closely related to the function of mitochondria. We obtained 1136 mitochondria genes (MTGs) from the MitoCarta 3.0 database (https://www.broadinstitute.org/) and then took intersections with risk model genes to obtain key genes. Following the pan-cancer analysis, key genes were assessed for differential expression and prognostic value across multiple cancer types. An assessment of the association between clinical features and TME of ccRCC was conducted using the TISIDB database (http://cis.hku.hk/TISIDB/). The ChEA3 database (https://maayanlab.cloud/chea3/#top) was used to obtain potential TF regulating key genes, and we selected the “Mean Rank” panel and took the top 10 genes for subsequent analysis and further screened the TF regulating key genes by differential analysis and K-M survival analysis. Patients with ccRCC were categorized into two groups according to the expression of key genes. After performing a differential analysis using the “limma” package, downstream pathways were identified using GO, KEGG, and GSEA analyses.

Cell lines used in this study were obtained from Procell Life Science&Technology Co., Ltd (Wuhan, China). An incubator containing 37°C and 5% CO2 was used to grow HK-2 and ACHN cells. Medium: MEM + 10% fetal bovine serum (FBS) + 1% antibiotics (HK-2 and ACHN), 1640 + 10% FBS + 1% antibiotics (786-O).

The overexpression plasmids and control plasmids of HIBCH were synthesized by Obio Technology (Shanghai) Corp., Ltd. We transfected the plasmid into 293T cells using calcium phosphate transfection to collect the viral fluid, which was then used to infect ACHN and 786-O cells, resulting in elevated levels of HIBCH expression in the cells. The RNA extraction was performed 48 hours after cell infection, and PCR was carried out to determine overexpression efficiency. Meanwhile, further cell phenotyping experiments were carried out.

For RNA extraction, we used TRIzol reagent (Invitrogen, Thermo Fisher Scientific, Inc.) (Eight pairs of ccRCC cancer and paracancerous tissue specimens were obtained from the Human Genetic Resources Center, The First Affiliated Hospital of Nanchang University.), and then reverse transcribed by Takara PrimeScript RT kit (Takara Bio, Inc., Otsu, Japan). The qPCR was performed on a Roche LightCycler96 real-time fluorescent quantitative PCR system using an SYBR premix Ex Taq kit (Takara Bio, Inc., Otsu, Japan). The relative expression of genes was calculated based on the 2^-ΔΔCt method. Primer sequences: HIBCH-F: 5’-GGAGTTGGTCTCTCAGTCCATG-3’, HIBCH-R: 5’-CCAAGTTTTCCTTGGAGTCGTGG-3’.

Cell migration was measured using 24-well transwell chambers, each upper chamber was inoculated with approximately 30,000 cells, and 200ul of FBS-free medium was added; the lower chamber was filled with 600ul of medium containing 20% FBS and counted after 36 hours. Cells were cultured to 80% density in 6-well plates, scratched, and then switched to an FBS-free medium and photographed in the same field of view at 0 h and 24 h, respectively.

Statistical analyses were conducted using R (version 4.2.2) or GraphPad Prism (version 9.0), and p<0.05 was considered statistically significant. Analysis of variance was used to compare categorical variables, and t-tests were used to compare continuous variables. Correlations between continuous variables were examined by Spearman or Pearson correlation analysis. Non-parametric samples comparing two independent samples were compared using Wilcoxon, while multiple independent samples were compared using Kruskal-Wallis.

By correlation analysis, we obtained 604 crosstalk genes, including 105 DEGs ( Figure 1A ). Then we removed 5 genes that were not detected in the validation cohort and finally obtained 100 DECGs. As shown in the volcano plot, there were 27 low and 73 high-expressed genes in the tumor tissue ( Figure 1B ). Subsequently, correlation network plots demonstrated a close association between DECGs ( Figure 1C ). the results of GO and KEGG analysis suggested that DECGs may be involved in biological processes such as protein modification and energy metabolism, and may play an important role in the HIF-1 signaling pathway ( Figure 1D ). Interestingly, the VHL/HIF pathway is linked to lactate production in ccRCC (13), which certainly suggests to us that these DECGs deserve to be studied in depth. According to Figure 1E , DECGs are mutated in 45.83% of samples, with MTOR showing the highest mutation frequency (18%).

Two subtypes of ccRCC were identified ( Figure 2A ), and the PCA confirmed this ( Figure 2B ). A comparison of the clinical characteristics of the two subtypes was conducted following that. In the TCGA cohort, age and gender did not differ between the two groups, whereas the distribution of grading and stage showed significant differences, with the C2 group having a higher nuclear grade and a more advanced clinical stage ( Figure 2C ). Although the E-MTAB-1980 cohort also showed the same trend, the C2 group in the cohort had more male patients ( Figure 2D ). The C2 group suffered a worse prognosis in both cohorts according to the K-M survival analysis ( Figure 2E ). In addition, the Sankey diagram more visually demonstrates the close association between subtypes and clinical features ( Figure 2F ).

The results showed that all scores, except for tumor purity, showed higher levels in the C2 group ( Figure 3A ). Figure 3B shows that the two subtypes infiltrated differently with immune cells, for example, the C1 group had more abundant monocytes and macrophages; the C2 group was infiltrated with more T cell follicular helper (TFH), T cell regulatory (Tregs), and NK cell activated. In addition, it is clear from the heat map that most of the immunoinhibitory, immunostimulatory, and HLA genes were expressed at higher levels in the C1 group ( Figures 3C–E ). In conclusion, all of the above results indicate that the ccRCC subtypes identified by DECGs have distinct TME and clinical characteristics. Therefore, an in-depth analysis of the predictive value of DECGs for the prognosis and immunotherapy of ccRCC is warranted.

For the risk model, LASSO regression analysis identified 11 prognostic genes in the training set ( Figure 4A ). All were associated with a better prognosis, except TTLL3 and CHFR ( Figure 4B ). Risk score per patient = SORL1*(-0.056) + HIBCH*(-0.122) + KDR*(-0.047) + VASH1*(-0.016) + VWA7*(-0.036) + TMEM25*(-0.202) + PLCL1*(-0.116) + PRUNE2*(-0.055) + TTLL3*0.069+CHFR*0.501+ABCG1*(-0.054). According to PCA, these 11 risk genes (RGs) could be assigned to different risk groups of ccRCC patients ( Figure 4C ). K-M survival curves show that the high-risk group has shorter long-term survival times ( Figure 4D ). The validation cohort also demonstrated similar results. Figures 4E, F demonstrate the distribution of survival status and RGs with a risk score, which is almost consistent with the trend in both cohorts.

Subsequently, we analyzed the risk characteristics of the different clinical variables. Among patients of different ages (<=65 vs >65), gender (Male vs Female), grades (G1/2 vs G3/4), and stages (Stage I/II vs Stage III/IV), the cumulative risk was increasing year by year in the high-risk group and was consistently higher than in the low-risk group ( Figures 5A–D ). We also found large differences in risk scores between variables within each variable. Although risk scores did not differ significantly between the two age groups, male, high-core graded (G3/4), and late-stage (Stage III/IV) patients tended to have higher risk scores ( Figures 5E–H ). In the Sankey diagram depicting the correlation between C1/C2 subtypes and risk groupings, we can observe that the C2 group with a worse prognosis is almost exclusively distributed in the high-risk group, whereas the C1 group, with a slightly better prognosis, is in a homogeneous distribution in the two risk groups (SF1 I).

In the TME score, both the immune score and the estimated score were higher in the high-risk group ( Figure 6A ). And the higher abundance of CD8 T cells and TFH clustered in the high-risk group ( Figure 6B ), which usually exerts anti-tumor immune effects. Additionally, only a few HLA genes showed differential expression, but most of the immunoinhibitory and immunostimulatory were different from them. High-risk patients expressed more CD96, CTLA4, IL10RB, LAG3, LGALS9, PDCD1, and TIGIT levels among immunoinhibitory molecules. The same is true for immunostimulatory, most of which are highly expressed in high-wind samples, such as CD70, IL6, and TNFRSF18 ( Figure 6C ). A similar expression pattern was observed in cohort E-MTAB-1980 ( Figure 6D ).

The immune-related function scores were generally higher in the high-risk group than in the other group ( Figures 7A, B ). Furthermore, the results of GO, KEGG, and GSEA also showed a strong association of risk grouping with immune-related functions. Biological pathways and functions related to immunity were enriched in DEGs between high- and low-risk groups ( Figure 7B ). Moreover, primary immunodeficiency and cytokine-cytokine receptor interaction pathways were enriched in the high-risk group ( Figure 7C ). The series of results suggest a close association between risk groups and the immune microenvironment, especially in high-risk groups.

Patients with ccRCC were analyzed based on TMB, MSI, IPS, and easier scores to predict their response to immunotherapy. High-risk samples showed significantly higher levels of TMB and easier scores, and the samples in the MSI group had higher risk scores ( Figures 8A–C ). Similarly, higher scores in the IPS scores about immune checkpoints were associated with higher risk scores ( Figure 8D ). According to the above results, TME differs greatly between the two risk groups, and the immunotherapy response may be more durable and effective in patients at high risk. Interestingly, we got the same results in a bladder cancer immunotherapy cohort. We found that patients in the immune-responsive group had higher risk scores and more high-risk patients ( Figure 8E ).

Using the training set, we performed univariate and multivariate COX regression analysis to identify independent predictors of OS, including age, stage, and risk score ( Figures 9A, B ). Based on the sum of the corresponding scores for each factor, a nomogram was constructed to predict patient survival at 1, 3, and 5 years ( Figure 9C ). Both the training and validation data sets showed the predicted probabilities to be almost in line with the actual probabilities ( Figure 9D ). Furthermore, the results of the ROC analysis also showed strong predictive performance of the model with 1-year AUC=0.863 for the TCGA cohort and 1-year AUC=0.900 for the validation cohort ( Figure 9E ). The DCA demonstrated that the nomogram was superior to the TNM staging for clinical purposes ( Figure 9F ).

From risk model genes, we identified HIBCH, a gene closely related to the mitochondrial function that may play a role in the development of ccRCC ( Figure 10A ). To begin with, HIBCH is expressed at lower levels in tumor tissues, and high levels are associated with better clinical outcomes ( Figures 10B, C ). High levels of HIBCH expression are associated with lower tumor grades and stages, reflecting its relationship to clinical variables ( Figure 10D ). The immune microenvironment and HIBCH also appear closely linked, which classify ccRCC into 6 immune subtypes and may be useful to classify different types of ccRCC according to their immune response ( Figure 10E ). In addition to having a negative correlation with immune cell infiltration, HIBCH expression was also found to be negatively correlated with TME scores ( Figures 10F, G ). More interestingly, the immune checkpoints CTLA4 and PDCD1 showed a negative correlation with the expression level of HIBCH as well ( Figure 10H ). According to the pan-cancer analysis, HIBCH displayed similar effects in numerous cancers (SF2). Based on these findings, there may be a mechanism through which HIBCH interacts with ccRCC’s immune microenvironment, which may influence tumor development and treatment response.

To further explore the regulatory mechanisms of HIBCH in ccRCC, we identified several potential TFs that regulate HIBCH expression, as well as potential pathways that inhibit ccRCC development. The co-expression heat map demonstrated the relationship between HIBCH expression and the top 10 TFs, most of which had high correlation coefficient values ( Figure 11A ). Difference analyses revealed significant differences between tumors and normal tissues in the expression levels of SPI1, GATA1, NR1H3, FLI1, SP2, MYBL2, and TFAP2C ( Figure 11B ). Subsequent K-M survival analysis of these DEGs revealed that FLI1, SP2, MYBL2, and TFAP2C were associated with OS ( Figure 11C ). Moreover, SP2 was lowly expressed in tumor tissues and associated with a good prognosis, while MYBL2 was highly expressed in tumor tissues and associated with a poorer prognosis. It is more likely that these two genes play a role in the regulation of HIBCH, but more experiments are needed to verify this conjecture. Following enrichment analysis, DEGs between high and low HIBCH expression groups had enriched biological functions associated with immunity ( Figure 11D ). More importantly, the GSEA results suggested that HIBCH is closely associated with FCGR-related pathways ( Figure 11E ). This gene family encodes the receptor for the Fc portion of immunoglobulin G, which is involved in a range of immune processes. This further suggests a complex mechanism of interaction between HIBCH and the immune microenvironment and has the potential to be a relevant biomarker for immunotherapeutic response.

Our in vitro studies revealed that HIBCH was higher expressed in HK2 than in ACHN and 786-O ( Figure 12A ). Moreover, by extracting RNA from kidney cancer and paraneoplastic tissues for qPCR, the same results were obtained, and the cancer tissues expressed lower levels of HIBCH ( Figure 12B ). Subsequently, we further explored the effect of abnormal expression of HIBCH on the migration ability of kidney cancer cells to elucidate its role in the metastasis of kidney cancer. Figure 12C showed that we successfully overexpressed HIBCH in ACHN and 786-O. As compared to the overexpression group, the number of cells in the control group was significantly higher ( Figure 12D ); in the scratch assay, the control cell migration rate was also higher ( Figures 12E, F ), and these results were statistically significant.