Computational construction of TME-related lncRNAs signature for predicting prognosis and immunotherapy response in clear cell renal cell carcinoma.

BACKGROUND: The tumor microenvironment (TME) is closely related to clear cell renal cell carcinoma (ccRCC) prognosis, and immunotherapy response. In current study, comprehensive bio-informative analysis was adopted to construct a TME-related lncRNA signature for immune checkpoint inhibitors (ICIs) and targeted drug responses in ccRCC patients.
METHODS: The TME mRNAs were screened following the immune and stromal scores with the data from GSE15641, GSE29609, GSE36895, GSE46699, GSE53757, and The Cancer Genome Atlas (TCGA)-kidney renal clear cell carcinoma (KIRC). And the TME-related lncRNAs were recognized using correlation analysis. The TME-related lncRNAs prognostic model was constructed using the training dataset. Kaplan-Meier analysis, principal-component analysis, and time-dependent receiver operating characteristic were used to evaluate the risk model. The immune cell infiltration in TME was evaluated using the single-sample gene set enrichment analysis (ssGSEA), ESTIMATE, and microenvironment cell populations counter algorithm. The immunophenoscore (IPS) was used to assess the response to immunotherapy with the constructed model.
RESULTS: In the current study, 364 TME-related lncRNAs were selected based on the integrated bioinformatical analysis. Six TME-related lncRNAs (LINC00460, LINC01094, AC008870.2, AC068792.1, and AC007637.1) were identified as the prognostic signature in the training dataset and subsequently verified in the testing and entire datasets. Patients in the high-risk group exhibited poor overall survival and disease-free survival than those in the low-risk group. The 1-, 3-, and 5-year areas under the curves of the prognostic signature in the entire dataset were 0.704, 0.683, and 0.750, respectively. The risk score independently predicted ccRCC survival based on univariate and multivariate Cox regression. GSEA analysis suggested that the high-risk group was concentrated on immune-related pathways. The high-risk group were characterized by high immune cell infiltration, high TMB and somatic mutation counters, high IPS-PD-1 + CTLA4 scores, and immune checkpoints expression upregulation, reflecting the higher ICIs response. The half inhibitory concentrations of sunitinib, temsirolimus, and rapamycin were low in the high-risk group.
CONCLUSION: The TME-related lncRNAs signature constructed could reliably predict the prognosis and immunotherapy response and targeted ccRCC patients' therapy.

INTRODUCTION

Clear cell renal cell carcinoma (ccRCC), the most common pathological subtype, accounts for 70%–80% of renal cell carcinoma (RCC). Nephrectomy accounts for a primary way to treat locoregional ccRCC. Nevertheless, in 30% of patients, ccRCC eventually develops into metastasis, leading to higher mortality and systemic therapy. Sunitinib, the first‐line targeted strategy, has been recommended for advanced ccRCC. However, sunitinib is unsatisfactory due to the emerged drug resistance. Immunotherapy, the essential antitumor treatment activating the tumor cell‐killing activity of immune system, has been considered as a possible method to treat or cure certain cancers. The immune checkpoint inhibitors (ICIs) that target immune checkpoints like CTLA‐4, PD‐1, and PD‐L1 have rapidly progressed during ccRCC treatment. Nevertheless, only a few ccRCC patients benefit from immunotherapy. Currently, there is a lack of effective biomarkers to predict immunotherapeutic response and guide therapy for ccRCC patients.

The tumor microenvironment (TME) consists of malignant cells and non‐transformed cells (cells of hematopoietic origin and cells of mesenchymal origin) mixed with stromal components. The interaction between tumor cells and TME contributes to cancer development and influences the prognosis and response to anti‐cancer therapies. , Understanding the potential functional mechanism in TME might help select treatment strategies to improve the survival of ccRCC patients.

Long noncoding RNAs (lncRNAs) have crucial roles in regulating diverse pathophysiological processes, including immune responses and cancer. , , LncRNA SNHG12 induces the progression and sunitinib resistance of RCC by upregulating CDCA3. Moreover, the expression pattern of lncRNAs is associated with cancer immunity and TME. Tumor‐infiltrating B lymphocyte‐related lncRNAs indicate immune cell infiltration in the TME and potential predictive biomarkers of immunotherapy response. HOTAIR can modulate cellular and non‐cellular components in the TME and contribute to tumor evolution and progression. According to our knowledge, the TME‐related lncRNAs model to predict prognostic outcome and response to immune as well as targeted therapy within ccRCC is inadequate.

In the current study, a set of TME‐related lncRNAs was identified, and a novel TME‐related lncRNAs prognostic and treatment‐related signature was constructed and applied for prognostic together with immune‐related and targeted therapeutic decisions in ccRCC patients.

MATERIALS AND METHODS

Data preparation

The Series Matrix Files of GSE15641 with 32 ccRCC samples, GSE29609 with 39 ccRCC samples, GSE36895 with 29 ccRCC samples, GSE46699 with 67 ccRCC samples and GSE53757 with 72 ccRCC samples were downloaded from the gene expression omnibus (GEO). The fragments per kilobase million (FPKM) values of 539 ccRCC and 72 healthy kidney specimens from The Cancer Genome Atlas (TCGA)‐kidney renal clear cell carcinoma (KIRC) were transformed into transcripts per kilobase million (TPM) values, which was considered as similar to GEO values. We secured 530 KIRC samples after eliminating duplicate samples. The batch effects in the TCGA and GEO datasets were corrected by the “ComBat” algorithm of the “sva” package. Finally, 507 ccRCC samples with overall survival (OS) of more than 30 days in TCGA‐KIRC were divided into training (356 cases) and testing (151 cases) datasets with a similar clinical characteristic distribution by the R package “caret” (Table S1). There were 505 ccRCC samples with disease‐free survival (DFS) information.

ESTIMATE‐based TME analysis

ESTIMATE was adopted for computing immune/stromal scores based on expression levels of specific genes.

Identification of TME‐related LncRNAs

The differentiated genes between the high and low immune/stromal scores were elected according to the screen criteria of |log2 fold change| > 0 and p < 0.05. Then, model genes closely related to immune/stromal scores were selected using the weighted gene co‐expression network analysis (WGCNA). We acquired 4061 TME‐related genes in the published research. A Venn diagram was used to select the intersected genes. Pearson correlation analysis with the standard as the correlation coefficient of more than 0.6 and p < 0.05 was adopted for identifying the potential lncRNAs related to TME‐related mRNAs, and the candidate TME‐related lncRNAs were used for performing further analysis.

Gene ontology (GO) as well as kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis

GO as well as KEGG analysis was conducted to explore the fundamental mechanism of the TME‐related mRNAs with R package “clusterProfifiler.”

Risk model construction and verification for KIRC

The differentiated TME‐related lncRNAs between 539 ccRCC and 72 healthy kidney species in TCGA‐KIRC were screened according to |log2 fold change| > 1 and p < 0.05. The prognostic TME‐related lncRNAs in the training data were detected by univariate Cox regression. The slightly contributory TME‐related lncRNAs were deleted by least absolute shrinkage and selection operator (LASSO) analysis. Finally, an optimal prognostic signature was developed using multivariate Cox regression. The risk scores were determined by risk score formula = in the training dataset (356 cases), independent testing dataset (151 patients), and the entire TCGA‐KIRC cohort (507 patients). The terms expi and coefi represent the expression and coefficient of the gene, respectively. KIRC patients in different datasets were separately divided into low‐ or high‐risk group according to median risk score. The separating capacity of the risk model was further verified using principal component analysis (PCA). Time‐dependent receiver operating characteristic (ROC) and Kaplan–Meier (K‐M) curve evaluated the risk model performance in predicting prognosis.

Association of risk model with clinical features

In total, 248 patients with sufficient clinical characteristics including age, gender, TNM stage, T, N, and M in the entire dataset were enrolled for univariate and multivariate Cox regression. Additionally, the Wilcoxon signed‐rank and chi‐square tests were applied to analyze relationship between the risk score and clinical features. The K‐M curve was used to detect the survival differences between both groups based on clinical characteristics‐stratified subgroups.

Construction of nomogram

The significant prognostic variables in multivariate Cox regression were applied for nomogram construction. This nomogram was later utilized for predicting outcome. Besides, calibration curves were utilized to evaluate the consistency between the predicted and actual survival results.

Gene set enrichment analysis (GSEA) enrichment analysis

GSEA was used to reveal GO and KEGG pathways in TCGA‐KIRC cases of diverse risk groups according to the gene sets “c5.go.v7.4.symbols.gmt” and “c2.cp.kegg.v7.4.symbols.gmt” using the R package “clusterprofiler.” , Adjusted p < 0.05 stood for statistical significance.

Immune infiltration analysis

The single‐sample gene set enrichment analysis (ssGSEA) was used to quantify the different enrichment degrees of 29 immunocytes as well as immune functions of two risk groups using the R package “GSVA.” ESTIMATE algorithm was utilized to determine stromal/immune/estimated scores and tumor purity for all cases in TCGA‐KIRC. In addition, differences of cytotoxic lymphocytes, CD8+ T cells, neutrophils, B cells, myeloid dendritic cells, monocytic cells, fibroblasts, and endothelial cells between two risk groups were subsequently analyzed by the Microenvironment Cell Populations (MCP) counter.

Somatic variant analysis

Genetic mutation data of TCGA‐KIRC calculated using the “Varscan” software were downloaded from the GDC, and the somatic variant landscape in the risk groups was analyzed and depicted by the R package “maftools.” TMB of all cases were counted based on the formula: (total mutation ÷ total covered bases) × 106.

IPS analysis

The immunophenoscore (IPS) positively related to tumor immunogenicity was considered a better predictor for ICIs response. The IPSs of KIRC patients were obtained from The Cancer Immunome Atlas (TCIA, https://tcia.at/).

Responses to targeted drug therapy

To evaluate the predicting value of the risk model in targeted drug sensitivity, TCGA‐KIRC samples were applied to analyze the half maximal inhibitory concentration (IC50) between two risk groups by the “pRRophetic” package.

Statistical analysis

R software (ver. 4.1.2) was employed for statistical analysis. Wilcoxon's signed‐rank test or Student's t test was utilized to compare differences between the two groups. For the variable comparison in more than two groups, either the Kruskal–Wallis test or one‐way ANOVA was employed. The chi‐square test was performed to analyze the frequency differences. The survival discrepancies were evaluated using the log‐rank test. The differences were regarded as significant when the bilateral p‐value was <0.05.

RESULTS

Identification of TME‐related lncRNAs

The immune/stromal components were quantitatively analyzed by ESTIMATE‐derived immune/stromal scores based on transcriptomic RNA‐seq data from 769 KIRC specimens derived from TCGA and GEO. The scope of immune scores was between −963.55 and 3632.46, whereas that of stromal scores was between −1390.33 and 2247.83 (Table S2). In the next step, these samples were classified into two clusters according to the median value. Under thresholds of |log2 fold change| > 0 and p‐value <0.05, the DEGs amounted to 5524 for high vs low immune scores and to 6039 DEGs for high vs. low stromal scores (Figure 1A‐D). The primary modules relevant to these immune/stromal scores were further identified using WGCNA. According to module trait heatmap, green model that contained 708 genes together with purple module that contained 269 genes showed higher relevance to immune and stromal scores, respectively (Figure 1E,F; R2 = 0.86 and p = 5e−219 in the green module, R2 = 0.72 and p = 1e−120 in the purple model). We obtained 4061 TME‐related genes from a previous study. Finally, 439 genes were screened out by intersection analysis (Figure 1G; Table S3). GO enrichment assessment revealed that the TME‐associated genes were closely associated with T‐cell activation, regulation of T‐cell activation, positive regulation of cytokine production, leukocyte cell–cell adhesion, and regulation of leukocyte cell–cell adhesion (Figure 1H). KEGG enrichment evaluation verified that these genes were related with the hematopoietic cell lineage, cytokine–cytokine receptor interaction, type I diabetes mellitus, allograft rejection, and viral protein interaction with cytokine and cytokine receptor (Figure 1I). Results suggested that these genes could be TME‐related genes. Finally, 364 TME‐related lncRNAs were selected with a threshold coefficient of correlation >0.6 and p < 0.05 (Table S4).

FIGURE 1

Identification of TME‐related lncRNAs in KIRC. Volcano plot (A) and heatmap (B) of DRGs between high and low immune scores. Volcano plot (C) and heatmap (D) of DRGs between high and low stromal scores. (E) The scale‐free fit index (left) and the mean connectivity (right). (F) Heatmap of the relationships between the state of immune and model eigengenes. (G) Venn diagram of TME‐related mRNAs. (H) The results of gene ontology analysis. (I) The top 15 most significant KEGG pathways

Establishment of a risk model

Under thresholds of |log2 fold change| > 1 and p‐value <0.05, we screened out 232 significantly and differentially expressed lncRNAs associated with TME for 539 ccRCC vs 72 healthy kidney specimens from TCGA‐KIRC (Figure 2A,B). The prognostic signature was constructed in the training dataset. 127 TME‐related lncRNAs were closely related with OS by univariate Cox regression (Table S5). 14 TME‐related lncRNAs screened by LASSO Cox regression analysis were further applied for multivariate Cox regression analysis (Figures 2C). Finally, a risk model was constructed based on six TME‐related lncRNAs to assess its value in predicting the prognosis of ccRCC patients. The univariate and multivariate Cox regression analyses suggested that the six TME‐associated lncRNAs were closely associated with OS except for AC084876.1 (Figures 2E,F). The risk score was estimated based on the coefficient and expression value of six TME‐related lncRNAs: Risk score = (0.471 expression of AC008870.2) + (0.634 * expression of AC068792.1) + (0.217*expression of LINC01094) + (0.075*expression of LINC00460) + (−0.386 *expression of AC007637.1) + (0.256*expression of AC084876.1).

FIGURE 2

Construction of the TME‐related lncRNAs risk model. (A) Volcano plot showing that 203 upregulated genes (red) and 29 downregulated genes (blue) between KIRC and normal kidney specimens. (B) Heatmap showing the expression levels of the top 20 upregulated and downregulated TME‐related lncRNAs. (C) The optimal values of the penalty parameter. Univariate (D) and multivariate (E) Cox regression analysis of the six TME‐related lncRNAs in the risk model

Validating the performance of the risk model

Depending on the median risk score of 0.926, the samples in the training dataset were separated into low‐ and high‐risk groups. PCA analysis showed that the signature had good performance in clustering (Figure 3A). The distribution of the risk score and survival status of the six TME‐related lncRNAs between two risk groups in the training dataset are shown in Figures 3B,C. Clearly, the high‐risk group contained more death samples. According to K–M plots, patients in low‐risk group showed better OS and DFS than those in high‐risk group (p < 0.001, Figure 3D,E). Additionally, the1‐, 3‐, 5‐year AUCs of the risk model were 0.738, 0.668, and 0.757, separately (Figure 3F). Furthermore, the prognostic model in the testing dataset and the entire dataset were also tested. Similarly, the low‐risk group showed a better clinical outcome than the high‐risk group in the testing dataset (Figure 3G,H) and entire dataset (Figure 3I,J), respectively. The AUCs of 1, 3, and 5 years were 0.637, 0.724, and 0.736 in the testing dataset and 0.704, 0.683, and 0.750 in the entire dataset (Figure 3K,L). Additionally, the expression levels of the AC008870.2, AC068792.1, LINC01094, LINC00460, and AC084876.1 in the TCGA‐KIRC were upregulated in the ccRCC samples, whereas the AC007637.1 was downregulated (Figure S1).

FIGURE 3

Validating the performance of the risk model. (A) PCA analysis for the risk model in the training dataset. Distribution of risk score (B) and survival status (C) in the training dataset. Kaplan–Meier curves of OS (D) and DFS (E) for the training dataset. The AUCs of the time‐dependent ROC curves for the training dataset (F). Kaplan–Meier curves of OS and DFS for the testing dataset (G, H) and entire dataset (I, J). The AUCs of the time‐dependent ROC curves for the testing dataset (K) and entire dataset (L)

Relationships between the risk model and clinical characteristics

The risk score and clinical features were analyzed using univariate and multivariate Cox regression over the whole TCGA‐KIRC dataset, revealing that our risk model independently predicted ccRCC prognosis (Figures 4A,B). In the clinical heatmap, tight correlations of grade, T, and stage with the risk model by the χ2 test were observed (Figure 4C). As revealed by scatter diagrams, stage, grade, and T had high‐risk scores based on Wilcoxon's signed‐rank results (Figures 4D‐F). The OS stratified by the clinical subgroups varied prominently between two risk groups, except N1 (Figure S2). To test the clinical practicability of the TME‐related lncRNAs signature, the two independent prognostic indicators yielded by the multivariable Cox analyses, the risk score and age, were incorporated to develop a nomogram for the prognostic prediction. Patients were assigned a total point based on each prognostic parameter in the nomogram (Figure 4G ). Higher total points presented a worse outcome. Calibration curves indicated that the nomogram exhibited a similar performance in predicting prognostic probability to an ideal model (Figures 4H ).

FIGURE 4

Assessment of the independent prognostic value of the risk model. Univariate (A) and multivariate (B) Cox regression analyses of the risk score and clinical characteristics. (C) Distribution landscape of clinical characteristics and the expression profiles of six TME‐related lncRNAs between the high‐ and low‐risk groups. Discrepancies in risk scores by grade (D), stage (E), and T (F). (G) The nomogram combining the risk score with the age in predicting 1‐,3‐, and 5‐year OS for KIRC patients given a total risk score. (H) Calibration curves. ***p < 0.001; **p < 0.001; *p < 0.05

GSEA analysis

The underlying mechanism between the two risk groups was interpreted using the GSEA. According to GO assessment results, enrichment of “activation of immune response,” “adaptive immune response based on somatic recombination of immune receptor built from immunoglobulin superfamily domains,” “antigen receptor mediated signaling pathway,” “B cell activation,” and “B cell mediated immunity” was observed in the high‐risk group (Figures 5A,B). Additionally, KEGG results revealed a high pathway enrichment, such as “autoimmune thyroid disease,” “chemokine signaling pathway,” “cytokine‐cytokine receptor interaction”, “primary immunodeficiency,” and “type I diabetes mellitus,” in the high‐risk group (Figures 5C, D).

FIGURE 5

Gene set enrichment analyses between the high‐ and low‐risk groups. GO enrichment analyses in the high‐risk group (A) and the low‐risk group (B). KEGG pathway analyses in the high‐risk group (C) and the low‐risk group (D)

Relationship between immune cell infiltration and the risk model

As shown in the GSEA results, a tight linkage of the risk model to the immunity was observed. To confirm the above‐mentioned conclusion, 29 immunocytes plus immunity functions were systematically evaluated through ssGSEA utilizing the TCGA‐KIRC cohort. A heatmap revealed the immune infiltration landscape (Figure 6A). The high‐risk patients scored prominently higher for ten immunocytes (e.g., activated dendritic cells, macrophages, CD8+ T cells, TILs, and T‐helper cells), whereas the low‐risk patients scored higher for mast cells (Figure 6F). In addition, the high‐risk patients scored higher for immunity functions, such as APC co‐activation or co‐repression, CCR chemotaxis, immunocheckpoints, Type I INF reaction, HLA (human leukocyte antigen), inflammation promotion, MHC class I, parainflammation, and T cell co‐activation or co‐repression. Low‐risk patients scored high for the Type II INF reaction (Figure 6G). In contrast, the high‐risk patients scored prominently higher ESTIMATE‐derived stromal/immune/estimate scores and lower tumor purity scores than the low‐risk patients (Figures 6B‐E). MCF counter was utilized for quantifying the cytotoxic lymphocytes, CD8+ T cells, mononuclear cells, neutrophils, B cells, myeloid dendritic cells, endotheliocytes, and fibroblasts. The CD8+ T cells, cytotoxic lymphocytes, B cells, and mononuclear cells had positively associated with the risk score, whereas neutrophils and endothelial cells exhibited negative association with the risk score (Figure 6H). Moreover, the high‐risk patients demonstrated higher numbers of B cells, cytotoxic lymphocytes, CD8+ T cells, and monocytic cells, whereas the low‐risk patients exhibited higher neutrophil and endotheliocyte counts (Figure 6I‐N).

FIGURE 6

The relationships between immune cells and the risk signature. (A) The heatmap of 29 immune‐related cells and functions, immune score, stromal score, ESTIMATE score, and tumor purity between two risk groups. (B‐E) The differences of immune score, stromal score, ESTIMATE score and tumor purity between two risk groups. (F) The scores of immune cells comparing high‐ and low‐risk groups by ssGSEA. (G) The scores of immune functions comparing high‐ and low‐risk groups by ssGSEAS. The correlation of the immune cells and the risk scores by MCP counter. (I‐N) The differences of immune cells between high‐ and low‐risk groups by MCP counter. ***p < 0.001; **p < 0.01; *p < 0.05; ns: no significance

Correlations between somatic mutation and risk model

The Varscan mutation annotation files of the TCGA‐KIRC cohort were used to analyze the total TMB and mutation distribution in two risk groups. The patients with TMB and somatic mutation information were assigned into high‐ and low‐risk groups. The TMB values were significantly higher in high‐risk group than those in low‐risk group and positively related with the risk scores (Figures 7A,B). Next, these samples were assigned into the H‐TMB group and the L‐TMB group using the median TMB. According to the Kaplan–Meier plot, the OS was inferior among the H‐TMB population (Figure 7C). Further, combination of TMB and the risk model had excellent risk stratification (p < 0.001, Figure 7D). Moreover, the somatic mutation counts in high‐risk group were significantly higher than those in low‐risk group and positively related to the risk scores (Figures 7E,F ). The waterfall plot implied the differences in mutation landscapes and the frequency of the first 20 mutated genes between two risk groups (Figures 7I,J). VHL (Von Hippel–Lindau Tumor Suppressor), PBRM1 (Polybromo 1), and TTN (Titin) had more than 10% mutation rate in both risk groups. The mutation rate of BAP1 was 15% and 5% among the high‐ and low‐risk populations, respectively (p = 0.004, Figure 7K). Patients with BAP1 mutation had significantly shorter OS than those with BAP1 wild (p < 0.001, Figure 7G). The risk model combined with BAP1 mutation status had excellent risk stratification (p < 0.001, Figure 7H).

FIGURE 7

Somatic mutational analyses between high‐ and low‐risk groups. (A) Difference in TMB between the high‐ and low‐risk groups. (B) Correlation between TMB and risk scores. (C) Survival analysis of OS between H‐ and L‐TMB groups. (D) Survival analysis of OS stratified by TMB and risk groups. (E) Difference in somatic mutation counts between the high‐ and low‐risk groups. (F) Correlation between somatic mutation counts and risk scores. (G) Survival analysis of OS between BAP1 mutation and wild groups. (H) Survival analysis of OS stratified by BAP1 mutation status and risk groups. The top 20 frequently mutated genes in the high‐risk group (I) and low‐risk group (J). (K) Proportion of BAP1 mutation status in high‐ and low‐risk groups. (L) Correlation between the risk score and common immune checkpoints. (M) Expression levels of the common immune checkpoints between the high‐ and low‐risk groups. ***p < 0.001; *p < 0.05; ns: no significance

The risk model for ICIs as well as targeted therapeutic responses

To validate whether higher TMB among high‐risk group implied superior outcomes with ICIs, we evaluated IPS, IPS‐CTLA4 (Cytotoxic T lymphocyte antigen‐4), IPS‐PD‐1, and IPS‐PD‐1 + CTLA4 to unravel the response to ICIs therapy in KIRC patients. The high IPS‐PD‐1 + CTLA4 scores among the high‐risk population suggested that the patients were better responders to ICIs (p = 0.0082, Figure 8A‐D). In addition, immune checkpoints correlated with an individual's response to cancer immunotherapy. The risk scores positively correlated with the PD‐1, LAG3, PD‐L1, TIGIT, and CTLA‐4 levels (Figure 8E). Besides, the high‐risk group also exhibited higher PD‐1, LAG3, TIGIT, CTLA‐4, and IDO1 levels (Figure 8F ). The association of risk model with the anti‐ccRCC targeted drug sensitivity were analyzed. High‐risk population showed lower IC50s of anti‐tumour drugs, like sunitinib (p = 1.7e‐09, Figure 8H), rapamycin (p = 4.7e‐08, Figure 8G), and temsirolimus (p = 4.3e‐12, Figure 8L). The above results demonstrated that our constructed risk model could predict the sensitivity to anti‐ccRCC ICIs and targeted therapeutics.

FIGURE 8

Response to immunotherapy and sensitivity to targeted therapy between the high‐ and low‐risk groups. (A‐D) The association between IPS and the risk signature of KIRC patients. (E) Correlation between the risk score and common immune checkpoints. (F) Expression levels of the common immune checkpoints between the high‐ and low‐risk groups. IC50 values between the high‐ and low‐risk groups for axitinib (G), sunitinib (H), sorafenib (I), rapamycin (G), pazopamib (K), and emsirolimus (L). ***p < 0.001; *p < 0.05; ns: no significance

DISCUSSION

The TME includes tumor cells and surrounding nontumor components, providing an optimum environment for tumor cell development. Research reveals that the immunosuppressive TME helps in tumor cell evasion from immunity and influences the therapeutic effect of antitumor drugs. Despite the immunogenicity of ccRCC, various TKI and ICI resistances are elicited using the exceptionally dynamic, heterogeneous, and adaptive TME. Like other immunogenic carcinomas, the total immune infiltration and tumor mutation load can serve as therapeutic response predictors, but their significance in clinical decision settings is unclear. The prognostic biomarker recognition for TME is probably conducive to define the potential checkpoint blockade responders and discovering novel therapeutic targets.

The study has demonstrated a crucial function exerted by lncRNAs in the TME on multiple tumor‐related processes. Many lncRNAs have critical role in regulating T‐cell function and immune response. Some stromal lncRNAs regulate cancer‐associated fibroblast (CAF) metabolism, enhancing the metastasis of tumor cells. LncRNAs can be considered as tumor biomarkers due to their dysregulation during the development and progression of cancer. Another study established that a TME‐related gene signature might help develop personalized immunotherapy. However, the TME‐associated lncRNAs prognostic signatures for the outcome and immunotherapeutic response forecasting among patients with ccRCC have been explored scarcely.

In the current study, the fused mRNA expression data were exploited to derive the immune/stromal scores using the ESTIMATE algorithm. We filtered TME‐associated mRNAs based on immune/stromal score disparities combined with the WGCNA result. As revealed by the GO enrichment assessment, these TME‐associated genes were linked tightly to the activation of T cells and regulation, positive regulation of cytokine generation, and leukocyte intercellular adhesion and regulation. The KEGG results identified enrichment of these genes in the hematopoietic cell lineage and the interplay between cytokine and cytokine receptors. Subsequently, 364 TME‐related lncRNAs were selected with the threshold of correlation coefficient >0.6 and p < 0.05. Based on the TME‐associated lncRNAs, we then developed a prognostic signature by exploiting the training dataset. Compared with the high‐risk population, the OS and DFS were superior among the low‐risk population. A similar conclusion was obtained in the test and entire datasets. Univariate and multivariate Cox regression analysis suggested that the risk score and age were individual outcome predictors for the survival from KIRC. The nomogram, including age and risk score, signified an excellent clinical application value in estimating ccRCC survival.

Among the six TME‐related lncRNAs, the expression levels of AC008870.2, AC068792.1, LINC01094, LINC00460, and AC084876.1 were upregulated in the ccRCC samples and that ofAC007637.1. LINC01094 facilitates ccRCC cell growth, metastasis, and radio resistance by the ceRNA mechanism. , LINC01094 enhances ccRCC cell growth, invasion and migration through downregulating miR‐184 and upregulating SLC2A3. High LINC01094 expression predicts dismal survival for gastric cancer (GC) cases with infiltration of macrophages. Studies have demonstrated that LINC00460 could affect cell proliferation, invasion, and migration in GC, head and neck squamous cell carcinoma (HNSCC) and colorectal cancer (CRC). , , , LINC00460 upregulation is related to clinical stage and inferior OS in ccRCC patients and accelerated cell proliferation, migration, and invasion through the ceRNA mechanism. The study combined with present studies has demonstrated that LINC01094 and LINC00460 might be oncogenic lncRNAs in ccRCC. No research is present about the AC008870.2, AC068792.1, AC007637.1, and AC084876.1 in ccRCC. So, further study must be performed to probe their functions in ccRCC.

Additionally, molecular mechanism between two risk models was further analyzed. GSEA‐based GO annotation indicated that the high‐risk group exhibited more “activation of immune response,”, “adaptive immune response based on somatic recombination of immune receptor developed from immunoglobulin superfamily domains,” “antigen receptor‐mediated signalling pathway”, “B cell activation,” and “B cell‐mediated immunity.” According to KEGG analysis, the “autoimmune thyroid disease,” “chemokine signalling pathway,” “cytokine‐ cytokine receptor interaction,” “primary immunodeficiency,” and “type I diabetes mellitus” were enriched in the high‐risk group. Therefore, TME‐related lncRNAs might influence the development and immune response in ccRCC via these immune pathways.

Using our risk model, a difference in infiltration of immunocyte was observed between two risk groups. The higher immune score for high‐risk patients might contribute to dismal survival. This observation was consistent with the former research that demonstrated that patients with high immune scores have poorer OS than those with low immune scores. CcRCC, as an immunogenic tumor, induces immune dysfunction by infiltrating immunosuppressive cells in TME. In our study, the infiltrations of six immunocyte types markedly increased in high‐risk group. Such conclusion is similar to the study, which has demonstrated that high‐risk glioma cases with increased infiltration of immunocytes presented a poor prognosis. The large amount of CD8+ T cells and cytotoxic T lymphocytes infiltration in the ccRCC TME cannot prolong survival due to immune evasion via T cell exhaustion. , CC motif chemokine receptors (CCR) and dendritic cells were high in the high‐risk group. The CCL/CCR axis has pro‐cancer properties in proliferation, migration, invasion, and angiogenesis. It can also recruit dendritic cells to increase the anti‐cancer response. In conclusion, the obtained results indicate that our constructed prognostic signature may have the ability to predict prognosis and infiltration of immune cells in ccRCC patients with clinical application value.

Currently, several KIRC clinical trials are in progress to estimate the efficacy of ICIs. The value of PD‐L1 for predicting response to ICIs in ccRCC is controversial. , TMB as a biomarker for predicting the efficacy of ICIs in ccRCC need further investigation. In the present study, the high‐risk cases showed increased TMB compared with low‐risk counterparts. BAP1 mutation in 10%–20% of ccRCC exerted a crucial function on genome stability. The high‐risk population displayed the prominently increased BAP1 mutation rate compared with low‐risk population. Additionally, BAP1 mutation contributed to poor survival among patients with ccRCC, in agreement with the prior work. Incorporating BAP1 mutation could improve the abilities of our signature in risk stratification. The IPS‐PD‐1 + CTLA4 score was remarkably increased among the high‐risk group. The risk score was positively associated with the PD‐1, PD‐L1, LAG3, TIGIT, and CTLA4 levels. The high‐risk population exhibited higher PD‐1, LAG3, CTLA‐4, IDO1, and TIGIT levels. These immune checkpoints directly inhibit the effector of T‐cell function and are responsible for exhausted T cells. , Thus, the high‐risk patients might have superior benefits from immunotherapy. Additionally, the IC50s of rapamycin, temsirolimus and sunitinib decreased in high‐risk population, which signified that these patients showed higher drug sensitivity. Taken together, our constructed risk model helped to predict the sensitivity to targeted therapeutics and immunotherapy in ccRCC patients.

Certain limitations should be noted. Firstly, the retrospective study and data bias on algorithm may influence the conclusion reliability. So, the prospective and multicentre trials are required to demonstrate the clinical availability. Secondly, further study is required to reveal the molecular mechanisms of the six lncRNAs in tumorigenesis and the development of ccRCC. Thirdly, the study of the synergistic effect of cell and spatial domains in antitumoral immunity processes was not mentioned.

In summary, the 6 TME‐related lncRNAs‐based risk model reliably predicted sensitivity to targeted therapeutics and immunotherapy among ccRCC cases.

AUTHOR CONTRIBUTIONS

L Zhou was responsible for study conception and design, data acquisition, data analysis, and drafting and revision of the article. H Fang was responsible for data analysis and drafting. G Fei was accountable for the study guide. M Yin and H Long were responsible for data acquisition and article revision. G Weng was responsible for the revision of the article.

CONFLICT OF INTEREST

The authors declared no conflicts of interest.

Figure S1

Click here for additional data file.

Figure S2

Click here for additional data file.

Table S1

Click here for additional data file.

Table S2

Click here for additional data file.

Table S3

Click here for additional data file.

Table S4

Click here for additional data file.

Table S5

Click here for additional data file.