The Role of Critical N6-Methyladenosine-Related Long Non-Coding RNAs and Their Correlations with Immune Checkpoints in Renal Clear Cell Carcinoma

Abstract

Purpose

This study aimed to evaluate the functions of critical N6-methyladenosine (m6A)-related long non-coding RNAs (lncRNAs) and their correlations with immunotherapeutic targets in clear cell renal cell carcinoma (ccRCC).

Methods

m6A-related lncRNAs were analyzed using the dataset from The Cancer Genome Atlas database via Pearson correlation analysis. Then, their prognostic functions in patients with ccRCC were determined via univariate Cox analysis. A prognostic m6A-related lncRNA signature (MRLS) in ccRCC was established using the least absolute shrinkage and selection operator (LASSO) Cox regression model. In addition, the correlations between these prognostic m6A-related lncRNAs with immune checkpoints were further evaluated in clinical samples.

Results

MRLS was established by the LASSO Cox regression model on the basis of seven prognostic m6A-related lncRNAs. The risk score for each patient was calculated using the MRLS model, and the patients were further stratified into high- and low-risk subgroups. The MRLS model was validated with a robust prognostic ability by the stratification analysis. On the basis of age, grade, stage, and risk score, a nomogram was developed with a strong reliability in forecasting the overall survival percentages of the patients with ccRCC. Moreover, seven prognostic m6A-related lncRNAs enrolled in the MRLS model were found to be correlated with various immunotherapeutic targets, namely, PD-1, PD-L1, CTLA4, and LAG3, and the expression levels of which in the high-risk subgroup were significantly higher than those in the low-risk subgroup. The significant correlations between LINC00342 and the aforementioned immunotherapeutic targets were also confirmed in clinical samples.

Conclusion

In this study, seven m6A-related lncRNAs were identified as potential biomarkers for forecasting the prognosis of patients with ccRCC and evaluating the efficacy of immunotherapy for these patients. Furthermore, a prognostic and predictive MRLS model with a high reliability was constructed to predict the overall survival probability of patients with ccRCC.

Keywords:

• immunotherapy
• N6-methyladenosine modification
• long non-coding RNA
• renal cell carcinoma
• prognosis

Introduction

Renal cell carcinoma (RCC), one of the most frequent cancer types in the urinary system,¹ accounts for approximately 3% of all kinds of cancer diagnoses and cancer-specific mortalities worldwide.^2,3 The incidence rate of RCC has increased annually by 2% worldwide during the last two decades,^2,3 translating into an average increase of 295,000 newly diagnosed RCC cases and 134,000 cancer-related deaths every year.^1,4 Clear cell RCC (ccRCC), the most common pathological type of RCC, accounts for approximately 75% of RCCs.⁵ Notably, the five-year survival rate of patients with RCC in clinical stage I is about 80–95%, but this rate dramatically decreases to less than 10% in patients with stage IV RCC,⁶ whose average overall survival length is only 10–15 months, indicating the great clinical significance of early diagnosis and accurate prognosis of patients with RCC to formulate the optimal therapeutic regimens and improve the patients’ quality of life. RCC is widely perceived as an immunogenic disease, but it is known to mediate immune dysfunction to a large extent by stimulating the infiltration of immune-inhibitory cells into tumor microenvironment.^7–9 Several underlying mechanisms have been developed to interpret the invalidity of the antitumor immune response caused by these multiple tumor-infiltrating cells, including inactivation of effector T cells and antigen-presenting cells via upregulation of suppressive factors, such as checkpoint molecules.⁸

Targeting immune suppression by inhibiting checkpoints has achieved positive therapeutic effects on some patients harboring ccRCC. To date, combinatorial methods containing checkpoint blockades are the standard choice for patients with advanced ccRCC.⁸ However, a substantial proportion of patients still do not benefit from checkpoint blockades. Under the circumstances, the identification of reliable biomarkers of response to checkpoint inhibition is of great clinical significance to facilitate improvements in the clinical efficacy of these therapies.

Many studies demonstrated that the immune checkpoints, including CTLA4, PD-1, PD-L1, and LAG3, could be effective targets for stimulating T cell-mediated antitumor immunity.⁸ However, up to date, no clinical biomarkers can accurately predict the prognosis of patients with ccRCC and evaluate the efficacy of immune checkpoint inhibitors. Hence, researchers were striving to break ground in intricate biological processes and molecular mechanisms of ccRCC to identify the biomarkers as potential prognosis factors in ccRCC and therapeutic targets.

Many biological processes, including the regulation of RNA post-transcriptional stability,¹⁰ localization,¹¹ translocation, slicing,¹² and translation,¹³ are modulated by various RNA modifications, among which RNA methylation accounts for over 60%.¹³ N6-methyladenosine (m6A) modification, an invertible and dynamical RNA epigenetic process regulated by m6A regulators containing “writers” (methyltransferases), “readers” (signal transducers), and “erasers” (demethylases),¹⁴ is the most prevalent epigenetic methylation of messenger RNAs (mRNAs) and noncoding RNAs (ncRNAs).^15,16 Recent studies reported that m6A modification participates in the regulation of oncogenesis and tumor progression in various kinds of malignancies,^17–19 including RCC.^20,21 For example, lncRNA DMDRMR cooperating with IGF2BP3 exhibits essential oncogenic roles in ccRCC in an m6A-dependent manner.²² Zhang et al²³ found that METTL14-mediated m6A modification negatively regulates the mRNA stability of bromodomain PHD finger transcription factor to drive lung metastasis in patients with RCC. In addition, both Zhou et al²¹ and Zhang et al²⁰ suggested that the genetic alterations of m6A regulators in RCC are strongly related to the potential prognosis of this disease.

The pullulation of transcriptome sequencing over the past decade has indicated that more than 70% of the genome is transcribed into RNAs, most of which are ncRNAs.²⁴ As a major type of ncRNAs equipped with more than 200 nucleotides in length, long ncRNAs (lncRNAs) are also considered to play a vital role in various biological processes with the prosperity of a wide range of researches in lncRNAs.^15,24,25 lncRNAs exhibit diverse biological functions through the regulation of gene expression levels and functions at the transcriptional, translational, and post-translational levels.²⁶ Dysregulation of lncRNAs is widely involved in the pathogenesis of many malignancies,^24,26,27 including ccRCC.²⁸ Accumulating evidences indicate that lncRNAs may serve as both prognostic biomarkers and therapeutic targets in many cancer types, especially in cancers in the urinary tract.²⁴ However, the connection between m6A regulators with the dysregulation of lncRNAs in cancers remains unclear. Thus far, no study has explored the underlying mechanisms by which m6A modifications contribute to lncRNA-dependent ccRCC genesis and evolution.

To close this gap of cogent evidence on the utility of m6A-related lncRNAs in forecasting the prognosis of patients with ccRCC and their potential as therapeutic targets for this disease, we determined the prognostic value of m6A-related lncRNAs via bioinformatics and statistical analyses of the data of patients with ccRCC from The Cancer Genome Atlas (TCGA) database. Univariate Cox regression analysis revealed 39 m6A-related lncRNAs with prognostic value in patients with ccRCC. Afterward, a prognostic m6A-related lncRNA signature (MRLS) was developed in accordance with the ability of seven m6A-related lncRNAs to predict the overall survival of patients with ccRCC. Then, an accurate nomogram was formulated with a high reliability to evaluate the overall survival of patients bearing ccRCC. The correlations of the seven critical m6A-related lncRNAs with immune checkpoints and the expression levels of immune checkpoint molecules in the high- and low-risk subgroups were evaluated to assess the efficacy of immune checkpoint inhibitors for the treatment of ccRCC. Furthermore, in order to evaluate the clinical values of the MRLS model, the focus next turned toward the expressions of some critical m6A-related lncRNAs immersed in the model and the connections between critical m6A-related lncRNAs and immune checkpoints in the TCGA samples and 20 ccRCC tissues and paired adjacent normal tissues of our own cohort.

Materials and Methods

Datasets and m6A-Related Genes

The Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/) and the cBioPortal website (https://www.cbioportal.org/) were employed to acquire the mRNA expression files (fragments per kilobase of transcript per million mapped reads [FPKM] normalized) and the corresponding clinicopathological data, respectively. Only ccRCC patients with data on overall survivals were enrolled in the final analysis to attenuate the impact of statistical bias on the results. The inclusion criteria were as follows: (1) patients were definitely diagnosed with ccRCC; (2) patients were followed up for overall survival. The exclusion criteria were as follows: (1) existences of other primary cancers; (2) patients without any follow-up data. Finally, 530 patients harboring ccRCC were extracted from the TCGA datasets. Thirty-five m6A RNA methylation regulators were selected on the ground of the results from previously published studies.^13,19,29 The expression matrices of 35 m6A-related genes were also collected from the TCGA datasets.

Annotation of lncRNAs

The lncRNAs extracted from the TCGA dataset were annotated in compliance with the lncRNA annotation file of the Genome Reference Consortium Human Build 38, which was acquired from the GENCODE website (https://www.gencodegenes.org/human/). A total of 10,876 lncRNAs were recognized in the TCGA dataset after the ensemble IDs of the genes were identified. Eight transcript forms covering incRNA, antisense, processed transcript, sense intronic, 3ʹ overlapping ncRNA, sense overlapping, and macro lncRNA of lncRNAs were included in the analysis.

Bioinformatics Analysis

The m6A-related lncRNAs were screened from the dataset via Pearson correlation analysis (Pearson R > 0.7 and p value < 0.001). Prognostic m6A-related lncRNAs were identified via univariate Cox regression analysis.

On the basis of the expression levels of the identified prognostic m6A-related lncRNAs, the patients bearing ccRCC were stratified into subgroups via consensus clustering by using the “ConsensusClusterPlus” package.¹³ Survival probability and KEGG pathway analyses were also carried out between the two subgroups by using R software 3.6.1. The PD-L1 expression and immunescore calculated for each patient were compared between different clusters. Furthermore, the relationship of the prognostic m6A-related lncRNAs identified herein with PD-L1 was evaluated via Pearson correlation analysis. The compositions of various immune cells within tumor microenvironment were also compared between different clusters.

Least absolute shrinkage and selection operator (LASSO) Cox regression was performed with the package “glmnet”¹⁵ to establish the optimal prognostic MRLS for patients afflicted with ccRCC that involved the seven m6A-related lncRNAs, which were perceived as indispensable biomarkers. MRLS risk score was calculated using the formula , where Expr_i represents the relative expression of m6A-related lncRNA i in the signature, and coef_i indicates the coefficient of lncRNA I. The entire TCGA cohort was randomly grouped into training and test sets at a ratio of 1:1. Afterward, the patients were divided into high- and low-risk arms at the median cut-off on the basis of MRLS risk scores. The MRLS risk scores were further validated in the test set via survival analysis, risk plot, and ROC curve. In addition, the correlations between risk scores and clinicopathological variables were investigated using R software 3.6.1.

Independent prognostic factors were confirmed by both univariate and multivariate Cox regression analyses with the combination of MRLS and clinical variables in the training and test sets. After that, a predictive nomogram was constructed on the basis of the selected independent prognostic variables by using R software 3.6.1. Three packages containing “rms”, “nomogramEx”, and “regplot” were employed to delineate the nomogram. Furthermore, the prognostic and predictive accuracy of the nomogram was evaluated with calibration curves.

The correlations of the seven critical m6A-related lncRNAs with immune checkpoints, including CTLA4, PD-1, PD-L1, and LAG3, were assessed via Pearson correlation analysis. Furthermore, the difference in the expression of each of these immune checkpoint molecules was analyzed between the high- and low-risk subgroups. In addition, the CIBERSORT algorithm, which is a tool for the deconvolution of the expression matrix of human immune cell subtypes according to the principle of linear support vector regression, was employed to evaluate the immune cell infiltration of the high- and low-risk subgroups.

The relative expression levels of LINC00342 in ccRCC and normal tissues, as well as in different subgroups, were examined and analyzed on the basis of different clinical characteristics by using the RNA sequencing data of TCGA samples in the UALCAN database. In addition, the connections between LINC00342 and the immune checkpoints were assessed by analyzing the RNA sequencing data from GEPIA. Furthermore, the expression levels of LINC00342, PD-1, PD-L1, LAG3, and CTLA4, were detected in the 20 ccRCC tissues and paired adjacent tissues of our own cohort to evaluate the correlations between LINC00342 and the immune checkpoints via RT-qPCR.

RNA Extraction and RT-qPCR

The use of clinical samples was approved by the Ethics Committee of the Second Affiliated Hospital of Nanchang University. The total RNA from the ccRCC tissues and adjacent normal tissues was extracted using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) in compliance with the manufacturer’s instructions. The reverse transcriptions of the RNA samples were completed with PrimeScript RT reagent kit (TaKaRa, Dalian, China). RNA expression was detected using SYBR Premix Ex Taq II Kit (TaKaRa, Dalian, China) on a StepOne qPCR platform (Applied Biosystems, Foster City, CA, USA). β-actin was used as an internal reference gene. The 2−ΔΔCt method was employed to quantify the relative expression levels. Each analysis was performed in triplicate. The primers used were as follows: LINC00342, 5ʹ-CCCAAAGCAGTCCTTCACTACA-3ʹ (forward) and 5ʹ-CTGCAGTTCACTCTGCTGCTT-3ʹ (reverse); β-actin, 5ʹ-CATGTACGTTGCTATCCAGGC-3ʹ (forward) and 5ʹ-CTCCTTAATGTCACGCACGAT-3ʹ (reverse); PD-1 (PDCD1), 5ʹ-CAGTTCCAAACCCTGGTGGT-3ʹ (forward) and 5ʹ-GGCTCCTATTGTCCCTCGTG-3ʹ (reverse); PD-L1 (CD274), 5ʹ-TACTGGCATTTGCTGAACGC-3ʹ (forward) and 5ʹ-TCCAGATGACTTCGGCCTTG-3ʹ (reverse); CTLA4, 5ʹ-AGCCACAGCTGGTGGTATCT-3ʹ (forward) and 5ʹ-GAATTGGGCCCATCGAACTG-3ʹ (reverse); LAG3, 5ʹ-CGACTAGAGGATGTGAGCCAG-3ʹ (forward) and 5ʹ-GATCCAGGTGACCCAAAGGAT-3ʹ (reverse).

Statistical Analyses

The statistical correlation between two variables was estimated via Spearman correlation analysis. Hazard ratios (HRs) with 95% confidence intervals (CIs) were calculated for each variable and independent prognostic factors were identified via univariate and multivariate Cox regression analyses. Statistical significance for the parameters between two groups was assessed by independent t-test. Corrections between the MRLS risk scores and clinicopathological variables were calculated via χ² test. The survival probabilities of patients with ccRCC in the two subgroups were evaluated via the Kaplan–Meier method with Log rank test by utilizing the “survminer” package in R. A nomogram and calibration curves were developed in line with the suggestions proposed by Iasonos et al.^24,30 The correlations of the seven critical m6A-related lncRNAs with the immune checkpoint molecules were validated via Pearson correlation analysis. All statistical analyses were conducted using R software 3.6.1. Statistical significance was defined a two-sided p value of < 0.05.

Results

A total of 14,086 lncRNAs and the expression matrices of 35 m6A-related genes were gathered from the TCGA dataset with the file downloaded from the GENCODE website (Figure 1A). Only when the expression value of a lncRNA was significantly associated with one or more of the 35 m6A-related genes (Pearson R > 0.7 and p value < 0.001) would the lncRNA be regarded as a m6A-related lncRNA. Identified from the TCGA dataset after applying Pearson correlation analysis were 356 m6A-related lncRNAs (Figure 1B), of who 39 were defined as prognostic m6A-related lncRNAs via univariate Cox regression analysis (Figure 1C). The expression levels of these m6A-related lncRNAs in normal and tumor tissues are displayed in Figure 1D.

Figure 1 (A) The components of m6A RNA methylation regulators in writer-, reader-, and eraser-complex; (B) network of m6A-related lncRNAs and m6A-related genes; (C) forest plot of 39 m6A-related prognostic lncRNAs; (D) the expression levels of 39 prognostic m6A-related lncRNAs in renal cell carcinoma tissue and normal tissue.

The impacts of the m6A-related lncRNAs identified herein on ccRCC progression were evaluated by dividing the tumorous tissues into different subgroups via consensus clustering according to the expression of the 39 prognostic m6A-related lncRNAs. The cumulative distribution function (CDF) of the consensus clustering with k = 2–9 and the relative change in area under the curve (AUC) increment are delineated in Figure 2A. Given the maximum AUC increment and the divergence in the expression corrections of the prognostic m6A-related lncRNAs between within groups and between groups, k = 2 was determined as the number of clusters (Figure 2B). The tracking plot of the subgroups for 2 to 9 clusters is presented in Figure 2C. The consensus matrix for optimal k = 2 is delineated in Figure 2D. The difference between these two clusters was discriminated by evaluating the survival probabilities of patients with ccRCC via survival analysis. Results showed that the patients in cluster 2 had better overall survivals than those in cluster 1 (p < 0.001) (HR: 1.860; 95% CIs: 1.375–2.516) (Figure 2E). The median PD-L1 expression in cluster 2 was significantly higher than that in cluster 1 (p < 0.01) (Figure 2F), but no significant difference was found between the two subgroups in terms of the median immunescore (p = 0.31) (Figure 2G). Intriguingly, CD8 T cells were filled in cluster 2 (Figure 2H). Moreover, 16 of the 39 prognostic m6A-related lncRNAs were found to be correlated with PD-L1 expression (Figure 2I). Gene set enrichment analysis (GSEA) revealed several enriched tumor signaling pathways in cluster 2, including the MTOR, Notch, phosphatidylinositol, and VEGF signaling pathways (Figure 2J–M). These results demonstrated that the stratification of the prognostic m6A-related lncRNAs identified herein was a clinically meaningful method, and their expression was strongly associated with the prognosis of patients with ccRCC.

Figure 2 (A) Consensus CDF for k = 2 to 9; (B) area under CDF curve increment for k = 2 to 9; (C) tracking plot for k = 2 to 9; (D) consensus matrix for optimal k = 2; (E) Kaplan–Meier overall survival curves for patients in distinct clusters (p < 0.001); (F) PD-L1 expression levels in distinct clusters (** p < 0.01); (G) immunescore for each patient in distinct clusters (p = 0.31); (H) various types of immune cells in distinct clusters; (I) Correlation matrix of PD-L1 and 39 prognostic m6A-related lncRNAs; (J–M) gene set enrichment analysis indicating that tumor signal pathways were enriched in the cluster 2 subgroup.

Abbreviations: CDF, clustering distribution function; HR, Hazard ratios; CI, confidence interval.

An MRLS model for evaluating the survival probabilities of patients with ccRCC was developed via LASSO Cox analysis on the basis of the 39 prognostic m6A-related lncRNAs identified herein. Seven lncRNAs coupled with their coefficients stuck out (Figure 3A–C). Each patient in the TCGA dataset was valued with a risk score on the basis of the coefficients for the seven lncRNAs. The patients in the low-risk subgroup had better overall survivals than those in the high-risk subgroup in both the training (p < 0.001; HR: 1.644; CIs: 1.430–1.889) and test (p < 0.001; HR: 1.512; CIs: 1.345–1.699) sets (Figure 3D and G, respectively). The risk plots of the training and test sets corroborated the prognostic reliability of the risk scores (Figure 3E and H). The AUCs of time-dependent receiver operating characteristic (ROC) curves revealed that the accuracy of prediction was 0.724 (95% CIs: 0.688–0.760) (Figure 3F) and 0.672 (95% CIs: 0.638–0.706) (Figure 3I) in the training and test sets, respectively. The relative expressions of seven critical m6A-related lncRNAs in the high- and low-risk subgroups of the training and test sets are shown in Figure 3J and K, respectively. All the immune checkpoint molecules, including PD-L1 (Figure 3L), PD-1 (Figure 3M), CTLA4 (Figure 3N), and LAG3 (Figure 3O), were related to more than three of the seven critical m6A-related lncRNAs, suggesting the strong connections between the immune checkpoints with the critical m6A-related lncRNAs. The Kaplan–Meier analyses revealed that the expression levels of all the seven critical m6A-related lncRNAs were significantly correlated with the prognosis of patients with ccRCC (all p values < 0.01) (Figure 3P–V).

Figure 3 (A–B) Least absolute shrinkage and selection operator (LASSO) regression, calculating the minimum criteria; (C) coefficients of selected lncRNAs; (D) Kaplan–Meier curves showed that the high-risk subgroup had worse overall survival than the low-risk subgroup in the training set; (E) distributions of risk scores and survival status of renal cell carcinoma patients in the training set; (F) receiver operating characteristic (ROC) curves of m6A-related lncRNAs significance for predicting the overall survival in the training set; (G) Kaplan–Meier curves showed that the high-risk subgroup had worse overall survival than the low-risk subgroup in the test set; (H) distributions of risk scores and survival status of renal cell carcinoma patients in the test set; (I) receiver operating characteristic (ROC) curves of m6A-related lncRNAs significance for predicting the overall survival in the test set; (J) heatmap of the expression of m6A-related lncRNAs in the high- and low-risk subgroups of the training set; (K) heatmap of the expression of m6A-related lncRNAs in the high- and low-risk subgroups of the test set; (L) correlation matrix of PD-L1 and 7 critical prognostic m6A-related lncRNAs; (M) correlation matrix of PD-1 and 7 critical prognostic m6A-related lncRNAs; (N) correlation matrix of CTLA4 and 7 critical prognostic m6A-related lncRNAs; (O) correlation matrix of LAG3 and 7 critical prognostic m6A-related lncRNAs; (P–V) the expression levels of AC108449.2, AF117829.1, AL008718.3, AL031670.1, CHROMR, COL18A1−AS1, and LINC00342 were significantly correlated with the prognosis of patients with ccRCC (all p values <0.01).

Abbreviations: HR, Hazard ratios; CI, confidence interval.

The connections of risk scores, clinicopathological variables, and expression levels of constituent m6A-related lncRNAs are presented in Figure 4A, which suggests that the patients in the high-risk subgroup were significantly related to a higher mean age, high mean immunescore, and a higher rate of T3-4 or G4 disease than those in the low-risk subgroup. Furthermore, the risk score was correlated with older age, G3-4 disease, and cluster 2 (Figure 4B–I). Patients in the high-risk subgroup had a worse overall survival than those in the low-risk subgroup in different stratifications (Figure 4J–U). These results demonstrated that the prognostic signature combining the seven critical m6A-related lncRNAs could independently forecast the prognosis of patients with ccRCC. In addition, the number of memory B cells, CD8 T cells, follicular helper T cells, and T regulatory cells (Tregs) in the high-risk subgroup was significantly higher than that in the low-risk subgroup (p = 0.008, p = 0.003, p < 0.001, and p < 0.001, respectively) (Figure 4V).

Figure 4 (A) Heatmap of the associations between the expression levels of the 7 critical m6A-related lncRNAs and clinicopathological features in the entire cohort; (B–I) risk scores of renal cell carcinoma patients with different clinicopathological features (including age, gender, grade, cluster, T stage, N stage, M stage, clinical stage); (J–U) Kaplan–Meier survival stratification analyses in multiple subgroups of renal cell carcinoma patients (including patients aged ≤65 or >65 years, male or female patients, patients with G1-2 or G3-4 tumors, patients with T1-2 or T3-4 tumors, patients with N0 or N1 diseases, patients with M0 or M1 diseases); (V) various types of immune cells in the low- and high-risk subgroups.

After all the clinicopathological features were integrated with MRLS, univariate COX regression analysis revealed that the variables, including patient age, tumor grade and stage, and risk score, were significantly predictive for patients with ccRCC in the training (Figure 5A) and test sets (Figure 5C). However, multivariate COX regression analysis indicated that tumor grade was no longer an independent factor for forecasting the prognosis of patients with ccRCC (Figure 5B and D). The survival probability of each patient was evaluated with a total point from each variable by using the nomogram established on the basis of these three independent parameters (ie, patient age, tumor stage, and risk score) (Figure 5E). The higher the total point, the lower the survival possibility of patient bearing ccRCC. The calibration curves regarding the 1-year (Figure 5F), 2-year (Figure 5G), and 3-year (Figure 5H) survival rates revealed that the predictive ability of the nomogram was similar to that of the ideal model, disclosing the high reliability of the nomogram for clinical utility.

Figure 5 (A and B) Univariate and multivariate Cox analyses of risk score along with clinicopathological features (including age, gender, grade, and stage) in the training set; (C and D) univariate and multivariate Cox analyses of risk score along with clinicopathological features (including age, gender, grade, and stage) in the test set; (E) nomogram based on risk score, age, grade, and stage in the entire cohort; (F–H) calibration plots of the nomogram for predicting the probability of overall survival at 1, 2, and 3 years in the TCGA dataset; (I–L) the expression levels of PD-L1, PD-1, CTLA4, and LAG3 in the high- and low-risk subgroup.

In addition, the expression levels of all immune checkpoint molecules, including PD-1 (Figure 5I), PD-L1 (Figure 5J), CTLA4 (Figure 5K), and LAG3 (Figure 5L), in the high-risk subgroup were significantly higher than those in the low-risk subgroup, suggesting the potential high efficacy of immune checkpoint inhibitors for the treatment of patients with a high risk score.

Consequently, the expression levels of LINC00342 in the ccRCC tissues were significantly increased than those in normal kidney tissues from the TCGA cohort (UALCAN) (Figure 6A). Furthermore, the relative expression levels of LINC00342 in the normal populations were significantly lower with that in ccRCC patients of any age at diagnosis, gender, and clinical stages (Figure 6B–D). In the meanwhile, significant correlations were also detected between the expression levels of LINC00342 and pathological grades and subtypes (Figure 6E and F). In addition, LINC00342 expression was positively correlated with immune checkpoints (PDCD1: R = 0.18, p = 2.7e-05; CD274: R = 0.13, p = 0.002; LAG3: R = 0.21, p = 7e-07; CTLA4: R = 0.32, p = 3.1e-14) (Figure 6G–J). In our own cohort, RT-qPCR showed that the expression levels of LINC00342 in the ccRCC tissues were significantly higher than those in the adjacent normal tissues (p < 0.001) (Figure 6K). The positive correlations between LINC00342 and the immune checkpoints were also confirmed in our own cohort (PDCD1: R = 0.363, p = 0.041; CD274: R = 0.423, p = 0.018; LAG3: R = 0.529, p < 0.001; CTLA4: R = 0.591, p < 0.001) (Figure 6L–O).

Figure 6 (A) The relative expressions of LINC00342 in normal individuals and renal cell carcinoma (RCC) patients (UALCAN); (B) the relative expressions of LINC00342 in normal individuals or in RCC patients with different ages (UALCAN); (C) the relative expressions of LINC00342 in normal individuals of either gender and in male or female RCC patients (UALCAN); (D) the relative expressions of LINC00342 in normal individuals or in RCC patients in stage 1, 2, 3, or 4 (UALCAN); (E) the relative expressions of LINC00342 in normal individuals or in RCC patients with Grade 1, Grade 2, Grade 3, and Grade 4 diseases (UALCAN); (F) the relative expressions of LINC00342 in normal individuals or in RCC patients with different tumor subtypes (UALCAN); (G–J) PD-1, PD-L1, LAG3, and CTLA4 genes positively correlated with LINC00342 in RCC patients (GEPIA); (K) the relative expressions of LINC00342 in 20 RCC tissues and paired adjacent normal tissues of our own cohort using qRT-PCR; (L–O) PD-1, PD-L1, LAG3, and CTLA4 genes positively correlated with LINC00342 in 20 RCC tissues and paired adjacent normal tissues of our own cohort using qRT-PCR.

Discussion

Increasing transcriptomic data stored in public databases, such as TCGA and Gene Expression Omnibus, have become available with the technical advancement in sequencing. By analyzing the profiles acquired from these publicly accessible databases or own prospectively maintained institutes, many studies found that m6A modification may serve as a leading factor contributing to cancer pathogenesis,¹⁵ and the expression of ncRNAs, including microRNAs,³¹ lncRNAs,²⁶ and circular RNAs,²⁴ could greatly impact tumorigenesis and cancer progression. However, the mechanism by which m6A modification promotes ccRCC progression in an lncRNA-dependent manner and the connections between m6A-related lncRNAs and immune checkpoints of the immunotherapy for ccRCC remain unclear. In this context, the TCGA dataset was reviewed to extract information regarding 530 ccRCC patients to identify prognostic m6A-related lncRNAs. Confirmed with prognostic significance in the TCGA dataset were 39 m6A-related lncRNAs, 7 of which were utilized to develop an MRLS model for evaluating the overall survival probabilities of patients with ccRCC. The stratification of these patients into high- and low-risk subgroups was done according to the risk score calculated based on the MRLS model. Patients in the low-risk subgroup exhibited an overall survival benefit compared with those in the high-risk subgroup, and GSEA analysis revealed that certain malignant pathways were enriched in the high-risk subgroup. Furthermore, multivariate Cox regression analysis confirmed MRLS as an independent biomarker for predicting the overall survival percentage of patients with ccRCC. A nomogram was constructed to predict the prognosis of these patients on the basis of three independent parameters, namely, patient age, tumor stage, and risk score. All the immune checkpoint molecules, including PD-1, PD-L1, CTLA4, and LAG3, were found to be strongly correlated with the seven critical prognostic lncRNAs and the risk score calculated based on MRLS, indicating the potential efficacy of immune checkpoint inhibitors for the treatment of high-risk patients.

In view of the potential function of lncRNAs in cancer pathogenesis, lncRNAs have been recently become a hot point in cancer research. Several lncRNAs, such as PCAT1,³² LINC00467,³³ and CCAT1,³⁴ play a vital role in prostate cancer progression and invasion. The significantly upregulated lncRNAs CASC11,³⁵ RMRP,³⁶ and UCA1²⁴ enhance cell proliferation and migration in bladder cancer. The overexpression of the lncRNAs UCA1²⁸ and EGFR-AS1³⁷ also promotes cell growth and metastasis in renal cancer.

m6A RNA methylation, the most common modification in RNAs, especially in miRNAs and lncRNAs, is closely related to the pathogenesis of various kinds of cancers.^17–19 An important underlying mechanism accounting for the effects of m6A regulators on cancer development is the modification of specific lncRNAs.¹⁵ The lncRNA GAS5, which is negatively regulated by the m6A reader YTHDF3, can inhibit colorectal cancer cell proliferation by interacting with and triggering YAP phosphorylation and degradation.³⁸ Lan et al³⁹ found that KIAA1429 can facilitate liver cancer progression via m6A-dependent modification of the lncRNA GATA3. The onset and development of cancer disease could be affected by the m6A-dependent modification of lncRNAs.¹⁵ Taking these evidences into accounts, it was reasonable to conclude that m6A modification could be targeted at specific lncRNAs to exert their functions. More attention should be devoted to investigating the interactions and functions of m6A modification and lncRNAs to identify prognostic biomarkers or therapeutic targets of cancers.

Notably, the immune escape mediated by the complex tumor microenvironment could result in the inefficiency of immunotherapy.⁴⁰ Many lncRNAs are reportedly upregulated during the development, differentiation, and activation of immune cells, including macrophages, dendritic cells, neutrophils, T cells, B cells, and bone marrow mesenchymal stem cells.⁴¹ lncRNAs are involved in various processes of the immune response in the tumor microenvironment, as well as in the activation of tumor immunosuppression and immunotherapeutic response in various cancers, such as RCC,^42,43 endometrial cancer,⁴⁴ and liver cancer.⁴⁵ Memory B cells, T follicular helper cells, and Tregs are important factors contributing to worse outcomes in patients suffering from ccRCC.⁴⁶ In this study, the relative expression levels of PD-1, PD-L1, CTLA-4, and LAG3 were significantly upregulated in the high-risk patients, who were more likely to achieve clinically positive responses to immunotherapy. The risk score was in a positive connection with the infiltration of memory B cells, CD8 T cells, T follicular helper cells, and Tregs, indicating that m6A-related lncRNAs participate in regulating the immune microenvironment, a presumption that may help in interpreting the worse outcomes of ccRCC patients in the high-risk subgroup.

The MRLS model was constructed on the basis of seven prognostic m6A-related lncRNAs, namely, AF117829.1, AC108449.2, CHROMR, AL008718.3, COL18A1−AS1, AL031670.1, and LINC00342, some of which had been reported to be in close connection with cancer pathogenesis. LINC00342 is overexpressed in nonsmall cell lung cancer^47,48 and colorectal cancer,^49,50 and it promotes cancer proliferation and progression as a competing endogenous RNA and functions as a potential biomarker of the prognosis of patients bearing these diseases. Wang et al⁵¹ found that COL18A1−AS1 is strongly associated with the overall survival of patients with ccRCC. The lncRNA CHROME is a central component of the noncoding RNA regulatory network maintaining cholesterol homeostasis,^52,53 whose imbalance is the underlying factor of several diseases, such as cancers.⁵⁴ However, few studies have investigated the impacts of these lncRNAs on ccRCC initiation and progression. Furthermore, studies on the interactions of lncRNAs and m6A-related genes are even scarcer. In this study, prognostic m6A-related lncRNAs were identified, and their impacts on the survival probability of patients with ccRCC and cancer development were evaluated.

Of note, several limitations cannot be ignored in this study. Our conclusions were drawn on the basis of the dataset extracted from one public database, namely TCGA, and the prognostic m6A-related lncRNAs selected herein should be further confirmed with other independent ccRCC cohorts. In addition, the interactions and functions of these prognostic lncRNAs and m6A-related genes were evaluated via bioinformatics analysis and thus should be validated by in vitro and in vivo experiments.

Conclusions

This study comprehensively evaluated m6A-related lncRNAs in ccRCC. Seven m6A-related lncRNAs were confirmed as potential biomarkers for forecasting the prognosis of patients with ccRCC and evaluating the efficacy of immunotherapy for these patients. Furthermore, a prognostic and predictive MRLS with a high reliability was constructed to predict the overall survival probability of patients with ccRCC, thus providing insights into the potential roles of m6A-related RNAs in ccRCC tumorigenesis and development.

Abbreviations

m6A, N6-methyladenosine; lncRNA, long non-coding RNAs; RCC, renal cell carcinoma; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristics; AUC, area under the curve; MRLS, m6A-related lncRNAs significance; HR, Hazard ratios; CI, confidence interval; TCGA, The Cancer Genome Atlas; RT-qPCR, Real-Time Quantitative PCR; GEPIA, Gene Expression Profiling Interactive Analysis.

Data Sharing Statement

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Ethics Approval and Consent to Participate

This study was conducted in accordance with the Ethical Standards of the Institutional Ethics Committee of Second Affiliated Hospital of Nanchang University and with the 1964 Helsinki declaration and its later amendments or comparable Ethical Standards. All patients provided informed consent.

Disclosure

The authors declare that they have no competing interests.

Additional information

Funding

This study was supported by the National Natural Science Foundation of P.R. China (Grant Nos. 81560419, 81960512, and 81760457), Jiangxi Provincial “Double Thousand Plan” Fund Project (Grant No. jxsq2019201027), Key Project of Natural Science Foundation of Jiangxi Province (20212ACB206013), and Youth Project of Natural Science Foundation of Jiangxi Province (20212BAB216037).