This article is mentioned in:

Introduction

Renal cell carcinoma (RCC), which originates from the tubular epithelial cells of the kidney, ranks among the ten most prevalent types of cancer globally and represent >90% of all renal cancer types. The most prevalent subtype of RCC is clear cell RCC (ccRCC), which accounts for 75–85% of cases (1). Recent statistics indicate that there were 434,419 newly diagnosed cases of RCC worldwide in 2022, accounting for 2.2% of all cancer cases, and 155,702 deaths from RCC, representing 1.6% of all cancer deaths (2). Although advancements in treatments have significantly improved the survival and quality of life for RCC patients, the prognosis for patients with advanced RCC remains poor, as evidenced by a 12% 5-year overall survival (2,3). Therefore, it is essential to focus on developing new treatment methods and optimizing personalized treatment strategies to enhance the prognostic outcomes for individuals suffering from advanced RCC.

Transfer RNA (tRNA) is a fundamental type of RNA molecule found in cells. The tRNA precursor produced through transcription is initially non-functional and must undergo a series of modifications and processing steps, facilitated by specific enzymes, to mature and acquire biological functions (4). These modifications, primarily catalyzed by tRNA modification enzymes, are crucial for tRNA's stability, folding and function, thereby affecting the efficiency and accuracy of protein translation (5). Other studies demonstrate that various RNA modifications [such as N1-methyl adenosine (m1A), 5-methylcytosine (m5C) and N7-methylguanosine (m7G)] and tRNA-modifying enzymes are often dysregulated in different types of cancer, such as colorectal and breast cancer (5,6). For example, m7G tRNA modification enhances the translation of oncogenes that serve pivotal roles in cell cycle regulation and the processes associated with the EGFR pathway (7). The dysregulation of tRNA-modifying enzymes, including methyltransferase-like protein 1 (METTL1) and WD Repeat Domain 4, is linked to poor prognosis in breast cancer and nasopharyngeal carcinoma (8,9).

The tumor microenvironment (TME) is essential in tumor progression and evolution (10). The TME comprises diverse cellular entities, including vascular cells, cancer-associated fibroblasts and infiltrating immune cells (11). These cells, in concert, orchestrate a myriad of processes including angiogenesis, the tumor cells invasion and evasion of growth-inhibitory factors, energy metabolism, immune evasion, cell proliferation and apoptosis (12–14). These functions are often executed in a manner that transcends individual cellular autonomy (15). tRNA modifications regulate immune cell function by affecting translation efficiency and accuracy. For instance, TRMT61A-mediated m1A modification of tRNA promotes the translation of key proteins (such as Myc) in CD4+ T cells by regulating codon decoding, thus ensuring a rapid immune response (16). Furthermore, METTL1-mediated m7G tRNA modification enhances the infiltration of cytotoxic immune cells into cancerous cells through the activating effector genes of the IFN signaling pathway, thereby impacting the TME (17). tRNA modification systems directly support cancer cell phenotypes through translational reprogramming, leading to increased proliferation, metastatic potential and cancer stem cell survival (18). Based on the aforementioned studies, it was hypothesized that the dysregulation of tRNA-related modified gene expression affects TME by enhancing the infiltration of immune cells, and thus promotes the occurrence and development of tumors. Tumor progression and patient survival reflect complex cellular and molecular interactions between the tumor and the host immune system (19). However, to the best of our knowledge, there is currently no reliable model to predict the prognosis of patients with ccRCC based on tRNA modification genes.

The present study aimed to apply univariate and multivariate Cox regression analyses to develop a prognostic model incorporating tRNA modification-related genes (TMRGs), in order to enable independent evaluation of ccRCC prognosis using The Cancer Genome Atlas (TCGA) database. Furthermore, the present study aimed to investigate the development of a nomogram, predict chemotherapy responses, perform somatic mutation analysis and functional enrichment analysis to elucidate the probable underlying mechanisms.

Materials and methods

Data acquisition and processing techniques

Data on relevant clinical information, nucleotide variations and mRNA expressions were obtained from the TCGA public database. The dataset GSE29609 (20) was sourced from the Gene Expression Omnibus (GEO) database (ncbi.nlm.nih.gov/). Leveraging the capabilities of the ‘GEOquery’ package (bioconductor.org/packages/GEOquery/,2.72.0), both gene expression and clinical datasets were extracted. After integrating the clinical data with the GEO transcriptome data by sample name, the final sample size consisted of 39 cases. A total of 39 samples derived from GSE29609 were employed as a validation cohort to assess the precision of the prognostic model. Additionally, 100 TMRGs were identified using the Molecular Signatures Database (MSigDB; https://www.gsea-msigdb.org/gsea/msigdb/human/genesets.jsp; Table SI). First, on the MSigDB website, ‘tRNA’ was searched, and ‘GOBP_TRNA_MODIFICATION’ was selected. Next, the gene set tab-separated values metadata was downloaded. Finally, a comparative analysis of the expression levels of tRNA-modified proteins within ccRCC tissues was performed utilizing the comprehensive resources obtained from the Human Protein Atlas (HPA) database (proteinatlas.org/).

Cluster analysis

Univariate Cox regression analysis was utilized to identify tRNA modification genes associated with prognostic outcomes. To investigate the correlation among TMRGs, a protein-protein interaction (PPI) network was created utilizing the capabilities of the STRING database (string-db.org/) as follows: i) Network type, full STRING network active; ii) interaction source: Text mining, experiments, databases, co-expression, neighborhood, gene fusion, co-occurrence; and iii) minimum required interaction score, medium confidence 0.400. Subsequently, the Cytoscape Cytohubba (cytoscape.org/,v3.10.2) plug-in identified 26 hub genes and modules in the PPI network related to tRNA modification genes. Cluster analysis was performed using the package ‘ConsensusClusterPlus’ (bioconductor.org/packages/ConsensusClusterPlus/1.66.0) program to identify molecular subtypes associated with tRNA modification; parameters included maxK, 9; reps, 10; and pItem, 0.8. Kaplan-Meier analysis was used to analyze the differences between the two groups. χ2 analysis was conducted concurrently to generate heat maps to systematically illustrate the associations between distinct clusters and their corresponding clinical features.

Identification and validation of the model

The least absolute shrinkage and selection operator (LASSO) and multivariable Cox regression analyses were performed to identify genes associated with ccRCC and establish prognostic signatures. The GraphPad software (version 9.5; Dotmatics) was employed to display the coefficients in the selected genes. The risk score was calculated as follows: ∑nn coef(i) × Expr(i) (where i represents genes). Kaplan-Meier analysis and receiver operating characteristic (ROC) curves were generated to assess the predictive value of each attribute. Cox regression analyses were applied to verify whether a signature was an independent risk factor. Correlation analysis, stratification analysis and the establishment of a nomogram, based on risk scores and pertinent clinical attributes, were performed according to clinicopathological standards. Calibration plots were created for the 1-, 3- and 5-year survival rates to evaluate the alignment between the predicted probabilities and the actual survival results.

Enrichment analysis

Gene Set Enrichment Analysis (GSEA; gsea-msigdb.org/gsea/index.jsp) was performed to examine pathway enrichment in high-risk groups. The internal reference gene sets included three categories: C2 KEGG, HALLMARK and C5GO. Significant results were indicated by normalized enrichment scores (NES) >1, a nominal P<0.05 and a false discovery rate (FDR) q-value <0.25.

Immunologic landscape analysis

A quartet of immune-related algorithms, including single sample (ss) GSEA, ESTIMATE, TIDE and CIBERSORT algorithms (RStudio IDE; Posit Software, PBC) were applied to assess and compare the immunological profiles of high-risk and low-risk cohorts. The ssGSEA methodology was utilized to elucidate the activity of immune cells, immune functions and the pathways associated with immunity for each sample. Marker genes for various immune cells were identified from previous studies, such as FOXP3/CTLA4, CD56, which are the marker genes of Tregs and Natural killer cells (21,22). The ESTIMATE method was utilized to calculate immune, stromal, estimate score and tumor purity by analyzing the ratio of immune cells to stromal cells. The composition of immune cell populations infiltrating each tumor specimen was predicted utilizing the CIBERSORT algorithm. Following cluster analysis and subsequent characterization, a comparative examination of the expression of the major histocompatibility complex (MHC) and immune checkpoint molecules was performed. Higher Tumor Immune Dysfunction and Exclusion (TIDE) scores are associated with longer survival and poorer response to checkpoint blockade therapy. Employing the TIDE database (tide.dfci.harvard.edu/), the TIDE scores for TCGA-kidney RCC (KIRC) cohort were calculated to predict immunotherapy responses in both subpopulations.

Analysis of tumor-related scores

The R package ssGSEA (parameters included: Method, ‘ssgsea’; kcdf, ‘Gaussian’; abs.ranking, TRUE) was applied to assess the angiogenic activity, tumorigenic cytokine score, mesenchymal EMT and stemness score for individual samples within the TCGA-KIRC cohort.

Gene mutation analysis

Based on TCGA-KIRC somatic mutation data, the ‘maftools’ (bioconductor.org/packages/maftools/,2.18.0) package was applied to analyze gene mutations. Tumor mutational burden (TMB) was calculated for each patient and compared between the high-risk and low-risk cohorts. The TMB scores were subsequently utilized to perform a comprehensive survival analysis. The cBioPortal database (cbioportal.org/v4.0) indicated that some genes selected by the signature have somatic mutations.

Chemotherapy response

The Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/) was used to evaluate the impact of predictive signatures on treatment response in ccRCC. The ‘oncoPredict’ (https://bioconductor.org/packages/oncoPredict, 1.0.1) package was employed to retrieve the gene expression profiles from the GDSC dataset with the corresponding drug response data. Sensitivity scores were applied to predict the maximum IC50 of all drugs in patients with KIRC.

Primer design

The primer design process in the present study typically adhered to the following: i) The nucleotide sequence of the target gene was retrieved from gene databases, such as National Center for Biotechnology Information GenBank (ncbi.nlm.nih.gov/genbank/), ii) based on the sequence of the target gene, an oligonucleotide sequence was designed that can anneal and bind effectively; the optimal length for primers generally ranges from 18 to 25 nucleotides, which ensures high specificity and appropriate annealing temperatures, iii) calculation of primer melting temperature (Tm) was performed, aiming for a Tm between 55–65°C as recommended; primers were calculated using online tools such as OligoCalculator (idtdna.com/calc/analyzer) and iv a design software such as Primer3 (tmcalculator.neb.com; v2.6.1) minimized non-specific binding between primer sequences and target sequences. Accession numbers are shown in Table SII.

Cell culture and reverse transcription-quantitative (RT-q) PCR

HK2 and ccRCC lines (786-O, 769-P, ACHN, Caki-1 and OS-RC-2) were purchased from the Shanghai Institutes for Biological Sciences. HK2 and ACHN cells were cultured in MEM (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (FBS; Gibco; Thermo Fisher Scientific, Inc.), while 786-O, 769-P and OS-RC-2 cells were cultured in RPMI-1640 (Gibco; Thermo Fisher Scientific, Inc.) medium supplemented with 10% FBS. Caki-1 cells were maintained in a McCoy's 5A medium (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% FBS. All cell lines were maintained in a sterile incubator (Thermo Fisher Scientific, Inc.) at 37°C with 5% CO2. Total RNA was extracted from the tissues or cell lysates of RCC patients utilizing a TRIzol™ reagent kit (Qiagen, Inc.). Subsequent cDNA synthesis employed a Thermo-script RT kit (Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. The 2(−ΔΔ C(T)) Method to analysis the relative gene expression (23). qPCR was conducted utilizing the CFX96™ Real-Time System (Bio-Rad Laboratories, Inc.) with SYBR Green PCR reagent (Takara Bio, Inc.). The thermocycling conditions were as follows: 95°C for 5 min, followed by 40 cycles of 95°C for 15 sec, 60°C for 25 sec and 72°C for 30 sec. GADPH served as a standardized internal reference for normalization. The primer sequences utilized for qPCR are listed in Table SIII.

Statistical analysis

All statistical analyses were performed using R software (version 4.3.1; Posit Software, PBC) and GraphPad Prism (version 9.5; Dotmatics). The data are presented as mean ± standard deviation (SD) or mean ± standard error of the mean (SEM), as specified in the figure legends. P<0.05 was considered to indicate a statistically significant difference. For survival analysis, Kaplan-Meier curves were generated and the log-rank test was applied to evaluate differences between groups. Univariate and multivariate Cox proportional hazards regression models were used to determine independent prognostic factors. The LASSO regression and multivariate Cox regression were performed to construct the prognostic model. Comparisons between two groups were conducted using the Wilcoxon rank-sum test or Student's t-test, depending on the normality of the data distribution assessed by the Shapiro-Wilk test. For comparisons among multiple groups, one-way ANOVA followed by Tukey's post hoc test was applied when the data were normally distributed, while the Kruskal-Wallis test followed by Dunn's post hoc test was used for non-normally distributed data. χ2 analysis was used to examine associations between categorical variables. GSEA was conducted using the clusterProfiler package (bioconductor.org/packages/clusterProfiler, 4.12.0). P<0.05 and FDR <0.25 were considered to indicate a statistically significant difference in terms of enrichment. The TIDE score differences between groups were assessed using the Mann-Whitney U test. The correlation between continuous variables was evaluated using Spearman's or Pearson's correlation analysis, depending on the distribution of the data. The ROC curve and area under the curve (AUC) were used to assess the predictive accuracy of the prognostic model. Differences in TMB scores were analyzed using the Wilcoxon test. Drug sensitivity prediction was performed using the oncoPredict package and IC50 values were compared using the Wilcoxon test. All experiments were performed with at least three biological replicates. Statistical significance for each analysis is indicated in the corresponding figure legends.

Results

Identification of two tRNA modification-associated clusters

A total of 100 genes (Table SI) associated with tRNA modification were identified in the MSigDB. Univariate Cox analysis identified 24 genes significantly associated with ccRCC prognosis, including six protective factors, leucine carboxyl methyltransferase 2 (LCMT2), SEPSECS, TRNA methyltransferase O (TRMO), TRNA methyltransferase 1L (TRMT1L), TRMT5, TRMT61B and 18 risk factors (Fig. 1A). The expression levels of the 24 genes related to prognosis in ccRCC tissues were compared with that of normal kidney tissues (Fig. 1B). The PPI network and correlation analysis illustrated the complex associations among the 24 prognosis-related genes (Fig. 1C and D). Clustering studies using these 24 genes (Table SIV) indicated that the optimal classification was to stratify patients with ccRCC into two categories (k=2) based on the expression patterns of 24 tRNA modification-associated genes (Fig. 2A-C). The cases represented by clusters 1 and 2 are shown in Table SV. Survival analysis of the two subgroups demonstrated that patients classified within cluster 1 exhibited a significantly more favorable prognosis compared with that of the patients in cluster 2 (Fig. 2D) The heat map significant associations between the clusters and clinicopathological parameters including clinical stage, T, M stage and grade, advanced stage (III/IV) or T4/N1/M1are enriched in Cluster2 (Fig. 2E).

Figure 1. Identification of two tRNA modification-associated clusters. (A) Univariate Cox analysis demonstrated 24 genes significantly associated with ccRCC prognosis. (B) The expression of the 24 prognosis-related genes in ccRCC tissues and normal kidney tissues. (C) The protein-protein interaction network of the 24 prognosis-related genes. (D) Correlation analysis demonstrated that numerous genes were interrelated (red indicates positive correlation, blue indicates negative correlation, and the depth of the color represents intensity). Error bars:Represent 95% confidence intervals. tRNA, transfer RNA; ccRCC, clear cell renal cell carcinoma.

Figure 2. Evaluation of tumor immune microenvironment in two clusters. Clustering analysis using the 24 prognosis-related genes stratified patients with ccRCC into two categories. (A) Consensus CDF in consistent clustering (k=2-9). (B) Relative change in area under the CDF curve from k2-9. (C) Consensus heatmap defining the two clusters (k=2). (D) Survival analysis of the two subgroups. (E) The heat map of gene expression differences was significantly correlated with clinicopathological parameters such as grade, clinical stage, T stage and M stage. (F) Cluster 2 was associated with lower expression levels of several major histocompatibility complex molecules. (G) Expression levels of immune checkpoints in Cluster 2. (H) Immune cell infiltration using CIBERSORT. (I) The ESTIMATE method showed that Cluster 2 exhibited higher tumor purity and lower stromal and ESTIMATE scores. (J) Cluster 2 exhibited significantly increased TIDE scores compared with cluster 1. CDF, cumulative Distribution Function, *P<0.05, **P<0.01, ***P<0.001.

Evaluation of tumor immune microenvironment in two clusters

In addition, cluster 2 exhibited significant association with decreased expression levels of various MHC molecules (Fig. 2F). Considering the disparities in immune infiltration between the two cohorts, the relationship with immune checkpoints was examined. Cluster 2 demonstrated significantly increased expression levels of immune checkpoints such as TNF receptor superfamily member 7, programmed cell death protein 1 (PD-1), T-cell immunoreceptor with Ig and ITIM domains (TIGIT), lymphocyte-activation gene 3 (LAG3), CD40 and cytotoxic t-lymphocyte antigen 4 (CTLA-4; Fig. 2G). CIBERSORT algorithm showed that cluster 2 exhibited a markedly increased infiltration of immune cells, encompassing CD8+ T lymphocytes, activated memory CD4+ T lymphocytes, follicular helper T cells and regulatory T cells (Tregs; Fig. 2H) Using the ESTIMATE methodology, cluster 2 demonstrated significantly increased tumor purity, accompanied by decreased stromal and ESTIMATE scores (Fig. 2I). Prolonged survival and diminished efficacy of checkpoint blockade therapy are linked to elevated TIDE scores (24). Cluster 2 exhibited significantly increased TIDE scores compared with that of cluster 1 (Fig. 2J).

Construction and validation of tRNA modification-related prognostic model

LASSO analysis in conjunction with multivariate Cox regression analysis was applied to refine the gene set of the present model, which resulted in the inclusion of six genes, FtsJ RNA 2′-O-methyltransferase 1 (FTSJ1), LCMT2, METTL6, pseudouridine synthase 1 (PUS1), TRMO and TRMT5, within the signature (Fig. 3A). The coefficients of the six genes in the signature are shown Fig. 3B. The relationships between risk score and six genes are illustrated in Fig. 3C. The risk score was calculated according to the formula: ∑nn coef(i) × Expr(i) and the model was then used to classify patients with ccRCC into low- and high-risk groups using optimized cut-off values (1.053). The risk score was correlated with clinicopathological parameters such as grade and clinical stage (Fig. 3D). At the 1-, 3- and 5-year AUCs (0.751, 0.726 and 0.739), patients that exhibited increased risk scores demonstrated significantly poorer prognosis compared with those with low-risk scores (Fig. 3E-G). The results of the survival analysis were validated through the utilization of the GSE29609 dataset (Fig. 3H-J). Additionally, changes in risk assessments among subgroups associated with numerous clinicopathological parameters were observed. Patients within grades 3 and 4, and stages T3-T4 and III–IV exhibited significantly increased risk scores, which suggested that more advanced tumors were associated with higher risk scores (Fig. 4A-C). Finally, Cox regression analyses demonstrated that the signature constituted an independent risk factor (Fig. 4D and E). From the Cox regression analysis results, age, score, clinical stage and signature were included when constructing the nomogram, with the gene signature as the key component of significance (Fig. 4F). According to the calibration plot, the survival times at 1-, 3- and 5-years exhibited alignment with the projected survival estimates (Fig. 4G).

Figure 3. Construction and validation of tRNA modification-related prognostic model. The 6 genes chosen for the signature by LASSO analysis and multivariate Cox regression analysis are shown in (A). (B) The coefficient of each gene in the signature. (C) The relationship between the risk score and the six genes. (D) The risk score was correlated with clinicopathological parameters such as grade and clinical stage. (E) Survival status, (F) heatmap and the (G) Receiver operating characteristic curve analysis together with the risk score in the TCGA. (H) Survival status, (I) heatmap and (J) ROC analysis together with the risk score were validated using the GSE29609 dataset. *P<0.05, **P<0.01, ***P<0.001. tRNA, transfer RNA; ccRCC, clear cell renal cell carcinoma.

Figure 4. Construction and validation of tRNA modification-related prognostic model. Patients in stage (A) T3-4, (B) stage III–IV and (C) grade 3–4 had higher risk scores. (D) Univariate and (E) multivariate Cox regression analyses. (F) Nomogram based on score, clinical stage, age and signature. (G) Nomogram calibration plots to predict the likelihood of survival for 1-, 3- and 5-years. tRNA, transfer RNA.

Functional enrichment and tumorigenesis score analyses for low- and high-risk ccRCC

GSEA was employed to examine pathways that govern tumorigenesis within the high-risk cohort. Numerous tumor-associated pathways were significantly enriched within the high-risk cohort (Fig. 5B). The findings demonstrated that the high-risk cohort exhibited a significant enrichment in the EMT pathway [NES=1.72; nominal (NOM) P=0.029; FDR q-value=0.12], IL6-JAK-STAT3 signaling pathway (NES=1.75; NOM P=0.027; FDR q-value=0.18), P53 pathway (NES=1.80; NOM P=0.018; FDR q-value=0.14) and the interactions of cytokine-cytokine receptors (NES=1.81; NOM P=0.014; FDR q-value=0.17; Fig. 5A). The scores of angiogenic activity, mesenchymal EMT, tumorigenic cytokines and stemness were assessed in patients with ccRCC. The scores for EMT were significantly decreased, while the assessments of tumorigenic cytokines, angiogenic activity and stemness exhibited a significant increase within the high-risk cohort compared with that of the low-risk cohort (Fig. 5C). The risk score demonstrated a positive correlation with the tumorigenic cytokines scores (R=0.22; P=8.4×10−07), angiogenic activity scores (R=0.25; P=2.4×10−08) and stemness scores (R=0.12; P=0.0092), whereas the risk score negatively correlated with mesenchymal EMT scores (R=−0.11; P=0.015; Fig. 5D).

Figure 5. Functional enrichment and tumorigenesis score analysis for low- and high-risk ccRCC groups. (A) The high-risk group showed significant enrichment of the EMT, IL6-JAK-STAT3, P53 and cytokine-cytokine receptor interaction pathways. (B) Gene set enrichment analysis of numerous tumor-related regulatory pathways enriched in the high-risk group. (C) The scores of angiogenic activity, mesenchymal EMT, tumorigenic cytokines and stemness in the high-risk group. (D) The risk score was positively correlated with tumorigenic cytokines scores (R=0.22; P=8.4×107), angiogenic activity scores (R=0.25; P=2.4×108) and stemness scores (R=0.12; P=0.0092), and negatively correlated with mesenchymal EMT scores. EMT, epithelial-mesenchymal transition; NES, normalized enrichment scores; FDR, false discovery rate; ccRCC, clear cell renal cell carcinoma; NOM, nominal. *P<0.05, ***P<0.001

Estimating immune cell infiltration and immune checkpoint inhibitors

Previous studies showed that the TME serves a pivotal role in tumorigenesis, and GSEA conducted in the present study indicated that several immune-related pathways were associated with the high-risk group. Therefore, the association with the tumor immune microenvironment was investigated. The ssGSEA algorithm demonstrated that the high-risk group showed significantly enhanced immune-related activities and pathways, as well as a more pronounced infiltration of immune cells in comparison with the low-risk group (Fig. 6A and D). Additionally, MHC expression levels were significantly increased in the high-risk group, compared with that of the low-risk group (Fig. 6B). The high-risk cohort exhibited significantly increased expression levels of immune checkpoint inhibitors, including TNFRSF9, PDCD1, TIGIT, LAG3, CD40 and CTLA-4 (Fig. 6C). The CIBERSORT algorithm indicated that plasma cells, CD8+ T cells, activated memory CD4+ T cells, Tregs, follicular helper T cells and M0 macrophages emerged as the predominant immune cells infiltrating the high-risk cohort (Fig. 6E and G). High TIDE scores in high-risk subtypes suggested that patients with high-risk ccRCC may exhibit a diminished response to immunotherapy (Fig. 6F). The ESTIMATE algorithm demonstrated that the high-risk cohort exhibited decreased tumor purity, and increased immunological and estimated scores compared with that of the low-risk group (Fig. 6H).

Figure 6. Immune cell infiltration and immune checkpoint inhibitors were estimated based on characteristics. (A) Analysis of immune-related activities or pathways and between the high-risk and low-risk group. (B) Major histocompatibility complex molecules expression level. (C) Expression of immune checkpoint inhibitors. (D) Analysis of greater immune cell infiltration between the high-risk and low-risk group. (E) The analysis of immune cell infiltration using CIBERSORT. (F) The TIDE scores between high-risk and low-risk group. (G) Heatmap of immune cell infiltration. (H) Analysis of tumor purity, stromal scores, immunological scores and estimated scores between the high-risk and low-risk group. *P<0.05, **P<0.01, ***P<0.001.

Somatic mutations and TMB scores in the signature

Primary nucleotide variation data pertaining to ccRCC were acquired from the TCGA database to investigate the genomic mutation disparities between high-risk and low-risk cohorts. Notably, the five genes exhibiting the highest mutation frequency within the low-risk group were identified as von Hippel-Lindau tumor suppressor (VHL; 47%), polybromo 1 (PBRM1; 43%), titin (TTN; 19%), mTOR (7%) and SET domain containing 2 histone lysine methyltransferase (SETD2; 7%; Fig. 7A). In the high-risk group, the top five genes with the highest mutation frequency were VHL (46%), PBRM1 (39%), SETD2 (20%), TTN (20%) and BRCA1-associated protein 1 (15%; Fig. 7B). Using a cut-off value is 0.925). The cases stratified into a low-TMB cohort exhibited a significantly prolonged survival duration compared with of the high-TMB cohort (Fig. 7C). The present model demonstrated that the high-risk and high-TMB cohort exhibited a markedly poorer prognosis when compared with that of the low-risk and low TMB group (Fig. 7D). Investigation into the mutation rates of marker genes demonstrated that FTSJ1 exhibited amplification mutations, while METTL6 was characterized by a predominance of profound deletion mutations. By contrast, LCMT2 displayed a higher frequency of missense mutations (Fig. 7E).

Figure 7. Somatic mutation and TMB in the signature. The somatic mutations in (A) the high-risk group and (B) the low-risk group are depicted in waterfall maps. (C) Analysis of the difference in overall survival between the group with high and low TMB levels. (D) Overall survival between in patients with high/low TMB levels and high/low risk. (E) The mutation rates of six marker genes. TMB, tumor mutational burden; FTSJ1, FtsJ RNA 2′-O-methyltransferase 1; LCMT2, leucine carboxyl methyltransferase 2; METTL6, methyltransferase-like protein 6; PUS1, pseudouridine synthase 1; TRMO, TRNA methyltransferase O; TRMT5, TRNA methyltransferase 5.

Drug sensitivity prediction

GDSC was used to predict the therapeutic response of high- and low-risk patient cohorts to widely utilized chemotherapeutic agents. Significant differences were observed in the responsiveness of the high- and low-risk cohorts to an array of chemotherapy agents, including tyrosine kinase (Fig. 8A), mTOR (Fig. 8B), AKT and Erk inhibitors (Fig. 8C). Additionally, the three-dimensional structures of potential drugs (sorafenib, axitinib, vistusertib, dactolisib, taselisib and afuresertib) were displayed using the PubChem database (Fig. 8D).

Figure 8. Drug sensitivity prediction. Drug sensitivity in low- and high-risk groups, including (A) tyrosine kinase inhibitors, (B) mTOR inhibitors and (C) AKT inhibitors and Erk inhibitors. (D) The 3D chemical structures of six potential target medicines. *P<0.05, ***P<0.001.

Expression levels of six signature genes in ccRCC tissues and cell lines

The protein expression levels of the six signature genes in renal cancer were derived from immunohistochemical staining data from the HPA database (Fig. 9A). To validate the results, total RNA was extracted from an array of ccRCC cell lines (786-O, 769-P, ACHN, Caki-1 and OSRC-2) and HK2 cells. Subsequently, the mRNA expression levels FTSJ1, LCMT2, METTL6, PUS1, TRMO and TRMT5 of mRNA were assessed. The RT-qPCR results showed that the mRNA expression levels of FTSJ1, METTL6 and PUS1 were markedly increased in ccRCC cells, while the mRNA expression levels of LCMT2, TRMO and TRMT5 were significantly decreased ccRCC cells compared with that of HK2 cells (Fig. 9B).

Figure 9. Expression of the six signature genes in ccRCC tissues and cell lines. (A) Immunohistochemical analysis of six signature genes. (B) Reverse transcription-quantitative PCR analysis of the six signature genes. ***P<0.001 and ****P<0.0001. The content displayed by the error bars represents the degree of dispersion of the data. The shorter the error bars, the more reliable the data. FTSJ1, FtsJ RNA 2′-O-methyltransferase 1; LCMT2, leucine carboxyl methyltransferase 2; METTL6, methyltransferase-like protein 6; PUS1, pseudouridine synthase 1; TRMO, TRNA methyltransferase O; TRMT5, TRNA methyltransferase 5.

Discussion

There is increasing evidence that the occurrence of RCC is a multi-step process characterized by the interplay among genetic, epigenetic and transcriptional changes (25,26). The development of RCC is often accompanied by mutations of the VHL gene, which inactivates the VHL protein and promote the accumulation of Hypoxia-inducible factor (HIF)-1α and HIF-2α, creating a favorable microenvironment for tumor cell growth (27,28). Concurrently, downstream target genes such as VEGF, platelet-derived growth factor, TGF-α and CXCR4 are activated, leading to increased tumor angiogenesis (29–31). Anti-angiogenic drugs such as sunitinib and sorafenib, which inhibit tumor angiogenesis, have been developed based on this mechanism (32,33). However, a large proportion of patients with RCC exhibit a short response to these targeted drugs and several patients develop drug resistance (34–37). Additionally, some patients experience severe adverse reactions during treatment, such as hand-foot skin reactions, hypertension and diarrhea (38–40). Therefore, studying molecular markers related to RCC metastasis, recurrence and prognosis is important in order to improve RCC treatment and patient survival rates.

tRNA-modifying enzymes serve a pivotal role in the occurrence, progression and treatment resistance in tumors by regulating the modification state of tRNA, which affects protein translation and cellular stress responses (41–43). The present study stratified patients with ccRCC into two distinct clusters, using the molecular subtypes of 24 TMRGs. The two clusters exhibited significant differences in immune infiltrating cells and clinical characteristics. Additionally, these clusters exhibited a strong association with the related immune pathways and tumor-related mechanisms. A tRNA modification-related signature composed of six genes (FTSJ1, LCMT2, METTL6, PUS1, TRMO and TRMT5) were also identified as predictors of clinical outcomes and treatment responses in patients with ccRCC. The present findings could potentially enhance the precision of survival probability forecasts for individuals afflicted by patients with ccRCC.

FTSJ1 is a 2′-O-methyltransferase responsible for adding a methyl group to the 2′-O position of tRNA, thereby modifying tRNA (44). In non-small cell lung cancer, FTSJ1 acts as a tumor suppressor by modifying tRNA and downregulating DRAM1 (45). Conversely, in triple-negative breast cancer, it functions as a tumor promoter by facilitating tumor progression and weakening the immune system's attack (46). However, the role of FTSJ1 in ccRCC is currently unclear. The present study used the common resource library (the GEPIA database, which demonstrated that LCMT2 (also known as TYW4) belongs to the highly variable methyltransferase superfamily. LCMT2 is identified as the putative homolog of the carboxy methyltransferase gene PPM2 from Saccharomyces cerevisiae, which serves a key role in the biosynthesis of the hypermodified guanosine known as wybutosine (47). LCMT2 catalyzes the final stage of wybutosine biosynthesis, namely methylation and methoxycarbonylation (48). LCMT2 is involved in DNA methylation, carrying mutations of intratumoral heterogeneity and microsatellite instability in colon cancer frameshift mutations (49). METTL6 is a tRNA methyltransferase that orchestrates the synthesis of 3-methylcytidine at the C32 site of the specific serine tRNA isoreceptor. METTL6 serves an important role in tumor cell proliferation; it was reported that knocking out METTL6 in mice led to a reduction in energy expenditure (50), and the expression levels of METTL6 are increased in luminal breast cancer (51).

PUS1 is a gene encoding an enzyme involved in the synthesis of pseudouridine, a modified base present in RNA (52). The pseudouridine modification has a role in RNA stability, translation and splicing, and is important for cellular functions such as protein synthesis. Disruption of PUS1 causes alterations in RNA processing and protein production, which may result in cancer cell proliferation and survival (53). Previous studies have suggested that PUS1 may serve a significant role in the development of various types of cancer, such as hepatocellular carcinoma, breast cancer and RCC (54–56). PUS1 promotes hepatocellular carcinoma through pseudouridylation of mRNA to enhance the translation of carcinogenic mRNA (54). The upregulation of PUS1 expression leads to an increase in the activity, migration, invasion and colony-forming ability of RCC cancer cells (55). Fang et al (56) reported that PUS1 could be utilized to predict adverse outcomes and triple-negative status in breast cancer, although these findings are preliminary and require further investigation to confirm its role in these specific types of cancer.

TRMO is a gene encoding an enzyme involved in RNA methylation, and its dysregulation may promote this process by affecting the methylation of tRNAs involved in the synthesis of carcinogenic or tumor suppressor proteins (57). TRMO dysregulation has been related to the susceptibility of differentiated thyroid cancer (58). TRMT5 is a nuclear coding protein involved in the post-transcriptional maturation of mitochondrial tRNA and mutations in TRMT5 can lead to complex hereditary neuropathic syndrome (59). A previous study has demonstrated the potential role of TRMT5 in hepatocellular carcinoma, targeting TRMT5 to inhibit hepatocellular carcinoma progression by inhibiting the HIF-1α pathway and enhancing sensitivity to adriamycin (60). The present study performed a comprehensive search using the PubMed database, utilizing the terms ‘GENE and ccRCC’, which indicated that the number of studies on the five characteristic genes associated with tRNA modification (FTSJ1, LCMT2, METTL6, TRMO and TRMT5) in ccRCC is limited. One of these studies reported that, by overexpressing PUS1 and knocking down PUS1 in RCC cells in vitro, PUS1 expression was associated with RCC cell viability, migration, invasion and colony-forming ability (55). Further experimental and clinical studies are required to confirm their association with the development of ccRCC.

In recent years, as the treatment of ccRCC has diversified, immunotherapy has emerged as a popular research area (61–63). The TME is a complex and dynamic milieu that encompasses not only tumor cells but also various extracellular matrix components, blood vessels and immune cells (64). The interplay between tumor cells and immune cells within the TME is instrumental in orchestrating the dynamics of tumor, immune evasion and response to cancer therapy (65). M2 macrophages are generally regarded as tumor promoters since they facilitate tissue remodeling, angiogenesis (formation of new blood vessels) and immunosuppression. M0 macrophages can be stimulated by Th2 cytokines to transform into M2 macrophages, which release cytokines such as IL-10 and TGF-β, inhibit the anti-tumor immune response and promote tumor cell proliferation (66). Tregs are immune cells that suppress the activity of other immune cells and maintain immune tolerance (67).

Previous studies have shown that Tregs are a type of immunosuppressive cell (68–70). In cancer, Tregs are frequently elevated in the TME and promote immune escape by inhibiting the functions of other immune cells that might attack the tumor (71–73). Tregs and macrophages M0 were the main immune cells infiltrating the high-risk group. Regulation of tRNA modifications may be associated with immune cell activation (74) In particular, specific tRNA modifications could facilitate the translation of genes that are crucial for maintaining the suppressive activity Tregs (16). For instance, tRNA modifications can influence macrophages from a resting state (M0) towards pro-inflammatory (M1) or anti-inflammatory (M2) phenotypes (75). Although the direct relationship between tRNA modification genes and the precise regulation of Tregs or M0 macrophages remains an area warranting further investigation, it is plausible that these modifications impact translational control over immune cell differentiation, function and metabolism. Research in this domain has the potential to yield novel insights into immune regulation and identify therapeutic targets for diseases related to immune dysfunction. The relationship between Tregs, M0 macrophages and the 6 genes identified (FTSJ1, LCMT2, METTL6, PUS1, TRMO and TRMT5) is currently unclear. The present results suggested that patients in high-risk subgroups of TMRGs are more susceptible to cancer immunosuppression. Tumor cells and immune cells within the TME typically express immune checkpoint proteins, such as PD-1 and CTLA-4, which inhibit T cell activation and facilitate immune tolerance (76,77). The upregulation of these checkpoints constitutes the primary mechanism through which tumors evade the immune response (78). Immunotherapies, such as immune checkpoint inhibitors (anti-PD-1 and anti-CTLA-4) have demonstrated successful outcomes in the treatment of certain types of cancer by restoring T-cell responses and overcoming immunosuppression within the TME (79). According to the present study, the expression levels of TNFRSF9, PD-1, TIGIT, LAG3, CD40 and CTLA-4 were significantly increased in the high-risk groups. Therefore, the tumor immune environment could be evaluated using immune checkpoint expression and the potential effectiveness of immune checkpoint inhibitors may be anticipated. TMB scores are a molecular marker to determine whether tumor patients are suitable for immunotherapy (80). Patients harboring tumors with a higher TMB score experience a significantly enhanced clinical benefit after receiving immune checkpoint inhibitors (81). TMB also showed strong efficacy in predicting overall survival.

In the present study, the low-TMB group exhibited a significantly prolonged survival compared with that of the high-TMB cohort. Furthermore, the high-risk + high-TMB subgroup demonstrated a significantly poorer prognosis compared with that of the low-risk + low-TMB group. Notably, there were significant differences in sensitivity between high- and low-risk groups to various chemotherapy agents; these results may potentially inform tailored treatment strategies for both high- and low-risk patients in the future. RT-qPCR was performed to measure the expression levels of the model genes in ccRCC cell lines. A total of six model genes were observed, which is consistent with a previous analysis. However, further experiments are needed, both in vivo and in vitro, to verify functional disparities among the model genes. Among the several prognostic genes associated with ccRCC, TMRGs remain understudied. To the best of our knowledge, there is no reliable model to predict the prognosis of patients with ccRCC based on tRNA modification genes. Compared with other prognostic models of ccRCC, the model in the present study is composed of fewer genes, had higher accuracy and may improves feasibility for future clinical and basic studies (82).

However, there are certain limitations of the present study. First, the present study is predominantly limited to publicly accessible databases for bioinformatic methods; therefore, it is necessary that further in vivo and in vitro experiments are conducted to elucidate the impact of these six model genes on the occurrence and development of the disease in patients with ccRCC. Second, the present study only analyzed the association between the risk model and immune cells, immune function, MHC molecules, immune checkpoints and immunotherapy. However, the significance of the model, particularly when integrated with immunotherapy, necessitates an ongoing commitment to accumulate a substantial array of samples for thorough evaluation. Finally, screening tRNA modified genes only based on MSigDB has certain limitations, such as i) data coverage: Although MSigDB contains a large gene set, it mainly focuses on known biological pathways and functional annotations and may not cover all genes associated with tRNA modification, particularly newly discovered genes or genes that have not been adequately studied; and ii) species diversity: MSigDB data mainly focuses on model organisms such as humans and mice, and may lack information on tRNA modified genes in other species, which limited the scope of the present study. However, the tRNA modified genes screened from the MSigDB are comprehensive, and to the best of our knowledge, no other database has been found to date with more comprehensive data. If more comprehensive or additional databases are found in the future, they may be used for additional validation in the future.

In conclusion, the present study identified a new prognostic model in ccRCC containing six genes (FTSJ1, LCMT2, METTL6, PUS1, TRMO and TRMT5), based on TMRGs. Furthermore, the disparities in immune profiles, genetic mutation statuses and pharmacological sensitivities across various molecular subtypes and risk categories were examined. The present findings may improve prognosis and facilitate personalized treatment strategies for patients with ccRCC, thereby enhancing individualized patient management.

Supplementary Material

Supporting Data

Acknowledgements

Not applicable.

Funding

The present work was supported by the Youth Project of Health Commission of Nantong City (grant no. QN2022017), Nantong Science and Technology Bureau (grant no. MS22019009) and Basic Research and Social Minsheng Plan Project (grant no. JC12022008).

Availability of data and materials

The data generated in the present study may be requested from the corresponding author.

Authors' contributions

BZ designed the present study. ZC, CS, XZ and WZ performed experiments. ZC, SYX, YFJ and XZ conducted the data analysis. XZ and ZC wrote the manuscript which was examined and revised by BZ. ZC and BZ confirm the authenticity of all the raw data. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References