Identification of a novel immune-related microRNA prognostic model in clear cell renal cell carcinoma.

BACKGROUND: Clear cell renal cell carcinoma (ccRCC) is a type of kidney cancer, and one of the most common malignant tumors. Many studies have shown that certain microRNAs (miRNAs) play an important role in the occurrence and development of ccRCC. Nevertheless, the prognosis of ccRCC patients is very rarely based on these "immuno-miRs". Our aim was thus to determine the relationship between immune-related miRNA signatures and ccRCC.
METHODS: We downloaded the miRNA expression data from 521 KIRC and 71 normal tissues in The Cancer Genome Atlas (TCGA). We used "limma" package and univariate Cox regression analysis to identify differentially expressed miRNAs (DEMs) that related to overall survival (OS). We applied lasso and multivariate Cox regression analyses to construct a prognostic model based on immuno-miRs. We evaluated the performance of model by using the Kaplan-Meier method. Furthermore, Cox regression analysis was used to determine independent prognostic signatures in ccRCC.
RESULTS: A total of 59 significant immuno-miRs were identified. We use univariate Cox regression analysis to acquire 18 immune-related miRNAs which were markedly related to OS of ccRCC patients in the training set. We then constructed the 9-immune-related-miRNA prognostic model (miR-21, miR-342, miR-149, miR-130b, miR-223, miR-365a, miR-9-1, and miR-146b) by using lasso and multivariate Cox regression. Further analysis suggested that the immune-related prognostic model could be an independent prognostic indicator for patients with ccRCC. The prognostic performance of the 9-immune-related-miRNA prognostic model was further validated successfully in the testing set.
CONCLUSIONS: We established a novel immune-based prognostic model of ccRCC based on potential prognostic immune-related miRNAs. Our results indicated that the 9-miRNA signature could be a practical and reliable prognostic tool for ccRCC. 2021 Translational Andrology and Urology. All rights reserved.

Introduction

Kidney cancer is a highly complex disease that may involve different types of kidney tumors, with renal cell carcinoma (RCC) being the most abundant amongst these. There are three types of kidney cancer: clear cell RCC (ccRCC), kidney renal papillary cell carcinoma (KIRP) or papillary RCC (pRCC), and kidney chromophobe (KICH) or chromophobe RCC (chRCC) (1). ccRCC is the most common type of kidney cancer, and is associated with the worst overall survival (OS) rate and high morbidity (2). The relationship between genetic factor alterations and the occurrence and development of ccRCC is not very clear, and thus the treatment options for ccRCC are relatively limited (3,4). Therefore, it is important to understand how certain prognostic factors affect the diagnosis and treatment of ccRCC patients.

MicroRNA (miRNA) is a single-stranded molecule of approximately 22 nucleotides, which are an important part of non-coding RNAs. miRNAs perform a vital biological process in inhibiting messenger RNA (mRNA) translation and expression by binding 3' untranslated regions (UTR) or 5' UTR (5). A great many miRNAs have been discovered by the latest human genome-sequencing technology, and an abundance of evidence suggests that miRNA controls the processes involved in the initiation and development of tumors, including cell proliferation, differentiation, growth, apoptosis, and aging, by regulating oncogenes or tumor suppressor genes (6). A number of studies have reported several miRNAs to be abnormally expressed in KIRC patients and capable of performing many important biological functions. Chen et al. distinguished between normal kidney tissues and ccRCC tissues by constructing a prognostic model of 11 deregulated miRNAs (7). Heinzelmann et al. suggested that specific miRNAs can distinguish between metastatic and non-metastatic ccRCC (8), while Wu et al. showed that a 4-miRNA-expression signature could determine the metastasis status of ccRCC patients (9). Some non-coding RNAs such as lncRNA and circle RNA play an important role in the development of ccRCC. For example, LINC02747 promotes the proliferation of ccRCC by inhibiting miR-608 and activating TFE3 (10). Wang et al. found that hsa_circ_0001451 was significantly down-regulated in RCC tissues and closely related to clinicopathological features and OS (11).

The immune system plays a vital role in the pathogenesis of tumors (12), including KIRC. Furthermore, a great many of immune-related miRNAs have been identified as being associated with KIRC. Some studies have shown that certain miRNAs act as key regulators of the immune response in different tumors. In particular, recent studies have showed that miRNA-mediated mechanisms regulate the activation of special immune cells in the tumor microenvironment. In addition, several miRNAs have been revealed to affect important cancer-related immune pathways that mediate immune cells to secrete immunosuppressive or stimulatory factors (13). Qu et al. identified miR-497-5p to be a potential therapeutic target and biomarker of ccRCC, and elucidated a novel regulatory mechanism of programmed death-ligand 1 (PD-L1) expression (14). Gigante et al. reported that miR-29b and miR-198 are abnormally expressed in T cells of RCC patients, indicating that targeting these miRNAs is beneficial to the diagnosis and treatment of RCC patients (15). Meanwhile, Jasinski-Bergner et al. reported that overexpression of miR-548q and miR-628-5p could cause a downregulation of human leukocyte antigen (HLA-G) and enhance the natural killer (NK) cell­-mediated HLA-G-dependent cytotoxicity, implying that the regulation of miRNAs can control the activity of immune cells. Therefore, miR-548q and miR-628-5p may be future therapeutic targets for RCC patients (16). Increasing evidence has shown that abnormal immuno-miRs are closely related to the occurrence and progression of ccRCC. There are some miRNA targeted therapy for ccRCC, for example, miR-144-3p and miR-21 expression can affect drug sensitivity of metformin, sunitinib and sorafenib. These drugs are common drugs used to treat RCC. Thus, miR-21 may be a potential therapeutic target. These data demonstrate that miR-21 may be the diagnosis, treatment, and prognosis biomarker candidate of ccRCC. Metformin can induce miR-21 expression (17). Gaudelot et al’s study showed that miR-21-silencing can increase drug sensitivity and decrease drug resistance in ccRCC (18). However, few, if any, definitive models currently exist that can systematically predict the tumor immune microenvironment and the overall prognosis of ccRCC patients based on immuno-miRs. Therefore, establishing a reliable and stable prognostic model has important clinical and therapeutic significance.

In this study, we thus obtained differentially expressed immune-related miRNAs (DEIMs) closely related to ccRCC through the “limma” package of The Cancer Genome Atlas (TCGA) and the Immune-miR database. Then, we further detected those immune-related miRNAs that were significantly related to prognosis. At last, we established an immune-related prognostic model by integrating immune-related miRNAs for ccRCC. We aimed to provide novel biomarkers that could effectively predict the prognosis of ccRCC patients.

Finally, we constructed a novel 9-immune-related-miRNA signature model based on TCGA database. The model can effectively predict the OS of ccRCC patients. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis revealed the biological function of miRNA-targeted genes and identified the significant signaling pathways related to ccRCC, which deepened our understanding of the pathogenesis of kidney cancer.

We present the following article in accordance with the TRIPOD reporting checklist (available at http://dx.doi.org/10.21037/tau-20-1495).

Methods

TCGA dataset collection

Level 3 miRNAs from 521 disease samples and 71 normal samples, mRNAs from 535 disease samples and 72 normal samples, and other information including clinical samples were downloaded from TCGA database (https://cancergenome.nih.gov/) (19). The detailed information of miRNA target genes was obtained from MiRTarBase (20). We acquired 245 immune-related miRNAs in online tools that contained specific miRNA-disease associations in the Immunology Database and Analysis Portal (Immune-miR; http://www.biominingbu.org/immunemir/index.html) database (21). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Identification of differentially expressed miRNAs (DEMs) and mRNA in ccRCC

We obtained miRNA and mRNA expression data of ccRCC samples from TCGA dataset. We then used the R package “limma” to screen DEMs and mRNAs (DEmRNAs) between tumor samples and normal samples (22). DEMs and DEmRNAs were considered using log2 |fold change (FC)| >1 and adjusted P value <0.01 as the screening criteria.

The construction of the prognostic immune-related miRNA signature

First, we used the R package “caret” to randomly divide the sample into two categories (training set and testing set). We applied univariate analysis to select the independent risk miNRAs in the training samples. In addition, we used the multivariate Cox regression method (survival package) to identify corresponding coefficients of the ccRCC prognostic signature (23). The risk score of every tumor sample from TCGA cohorts was derived from the immune-related miRNA signature. All tumor samples were randomly divided into high- and low-risk groups with the median of the risk score. We used a Kaplan­-Meier curve and log-rank test to evaluate survival analysis in the training, testing and whole set. The predictive performance of the immune-related miRNA signature was evaluated by the area under the curve (AUC) of 1-, 3-, and 5-year overall survival dependent on the receiver operating characteristic (ROC) curve according to the R software package, “survivalROC” (24).

Functional enrichment analysis of miRNA target genes

The relationship between miRNAs and target genes was obtained from miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/). We applied the intersection between the differentially expressed genes (DEGs) and the miRNA target genes to clarify the molecular mechanism of the important genes in kidney cancer. Finally, we used the R package, “clusterProfiler”, (25) for functional enrichment analysis based on the intersection genes.

Statistical analysis

All statistical analysis were performed on R software (version 3.6.1). DEMs and mRNAs (DEmRNAs) were obtained by R package “limma”. Training set and testing set of tumor samples were divided by R package “caret”. Univariate Cox regression analysis were used to assess the relationship between miRNA and overall survival. The 9-miNRA model was established based on multivariate Cox proportional hazards regression model. Both the univariate and the multivariate Cox regression analysis were performed with the ‘survival’ package. The ROC curves and corresponding area under the curve (AUC) were drawn by the package of ‘survivalROC’ in R.

Results

Identification of differentially expressed immune-related miRNAs and mRNAs in ccRCC

We selected 59 DEIMs (|log FC| >1, P<0.01 after P value adjustment) between ccRCC tissues and normal tissues. Among the miNRAs, 26 miRNAs were upregulated and 33 were downregulated (). The cluster heatmap of 59 DEIMs is displayed in . In addition, there were a total of 6,204 DEGs, including 2,634 upregulated and 3,579 downregulated genes based on 535 tumor samples and 72 normal samples.

Figure 1

Identification of DEIMs using “limma” package in R software. (A) Volcano plot showing the immune-related miRNA expression change between ccRCC and normal tissue. The cutoff was log|FC| >1, adjusted P<0.01. (B) The cluster heatmap of the DEIMs between ccRCC tissues and normal tissues. DEMs, differentially expressed immune-related miRNAs; FC, fold change; ccRCC, clear cell renal cell carcinoma.

Screening of immune-related miRNAs with the prognostic value in ccRCC

All tumor samples were randomly divided (N=521) based on miRNA expression profiles into training samples (N=261) and testing samples (N=260). To identify DEIMs with prognostic characteristics, the expression of 59 miRNAs in the training group was evaluated by univariate Cox analysis in. Ultimately, we found that 18 immune-related miRNAs were closely related to OS. The details of immune-related miRNAs are shown in .

Table 1

General characteristics of prognostic immune-related miRNAs

miRNA symbol	HR	Z	95% CI	P value
hsa-miR-21	2.54	4.0	1.61–4.00	6.5E-05
hsa-miR-106b	1.99	2.9	1.26–3.15	3.4E-03
hsa-miR-155	1.68	2.3	1.08–2.63	0.022
hsa-miR-142	1.58	2.0	1.01–2.46	0.044
hsa-miR-532	0.49	−3.0	0.31–0.78	0.0024
hsa-miR-193a	1.86	2.7	1.19–2.92	0.0068
hsa-miR-342	1.72	2.4	1.10–2.68	0.017
hsa-miR-10a	0.63	−2.0	0.40–0.99	0.046
hsa-miR-149	2.03	3.0	1.28–3.21	0.026
hsa-miR-130b	2.18	3.3	1.38–3.45	9.1E-04
hsa-miR-223	1.77	2.5	1.13–2.78	0.013
hsa-miR-365a	2.10	3.2	1.32–3.32	0.0016
hsa-miR-9-1	3.08	4.4	1.87–5.07	9.1E-06
hsa-miR-204	0.44	−3.6	0.28–0.69	3.6E-04
hsa-miR-146b	2.06	3.1	1.31–3.22	0.0017
hsa-miR-31	1.93	2.8	1.23–3.04	0.0046
hsa-miR-451a	0.54	−2.6	0.34–0.86	0.0086
hsa-miR-183	1.99	3.0	1.27–3.14	0.0029

HR, hazard ratio; CI, confidence interval.

Construction of an immune-related prognostic model for ccRCC

We finally selected 9 miRNAs to construct a prognostic model in training samples based on the multivariate Cox analysis. The detailed formula of the model was as follows: risk score = (0.41 × expression level of hsa-miR-21) + (–0.06 × expression level of hsa-miR-342) + (0.02 × expression level of hsa-miR-149) + (0.47 × expression level of hsa-miR-130b) + (0.27 × expression level of hsa-miR-223) + (0.04 × expression level of hsa-miR-365a) + (–0.03 × expression level of hsa-miR-204) + (–0.12 × expression level of hsa-miR-14b) + (0.08 × expression level of hsa-miR-91). The risk scores were obtained by combining the expression level of miRNA with its corresponding regression coefficients. The AUC values of the ROC curve for 1-, 3-, and 5-year OS were respectively 0.714, 0.751, and 0.72 in the training set, 0.715, 0.667, and 0.680 in the testing set, and 0.715, 0.707, and 0.697 in the whole set. The results suggested that this prognostic model based on immune-related miRNAs demonstrated a good sensitivity and specificity (). We applied multivariate Cox regression to obtain regression coefficients. The risk scores were related to miRNA expression and regression coefficients. Based on the median risk score, we divided the KIRC patients into a high-risk group and low-risk group. Survival analysis revealed that the survival rate of patients in the high-risk group was lower than that in the low-risk group in each of the training, testing, and whole sets (P<0.0001; ). In the training set, we applied a univariate and multivariate Cox regression model to characterize miRNAs associated with OS and tumorigenesis (). The Kaplan-Meier method identified hsa-miR-21, hsa-miR-342, hsa-miR-149, hsa-miR-130b, hsa-miR-223, hsa-miR-365a, hsa-miR-9-1, hsa-miR-204, and hsa-miR-146b as the 9 miRNAs that were closely related to patients’ OS (P value <0.05; ). Among these 9 miRNAs, 8 miRNAs [hsa-miR-21, hsa-miR-342, hsa-miR-149, hsa-miR-130b, hsa-miR-223, hsa-miR-365a, hsa-miR-9-1, and hsa-miR-146b; hazard ratio (HR) >1] were negatively associated with survival and 1 miRNA (hsa-miR-204, HR <1) was positively associated. This indicated that these 8 miRNAs had high-risk characteristics, because their high expression levels showed that patients had a shorter OS, but low expression levels of 1 miRNA indicated a shorter patient OS ().

Figure 2

ROC analysis based on time and survival analysis for the 9-miRNA signature in ccRCC. (A) ROC analysis of 1-, 3-, and 5-year OS in the training set, (B) testing set, and (C) the whole set. (D) The risk score and Kaplan-Meier analysis for the 9-miRNA model in the training set, (E) testing set, (F) and the whole set. ROC, receiver operating characteristic; ccRCC, clear cell renal cell carcinoma; OS, overall survival.

Figure 3

Survival analysis identified 9 immune-related miRNAs related to OS in ccRCC patients. According to the median expression of miRNAs, patients were divided into high-expression groups and low-expression groups. (A) OS for patients with miR-21; (B) OS for patients with miR-130b; (C) OS for patients with miR-9-1; (D) OS for patients with miR-342; (E) OS for patients with miR-223; (F) OS for patients with miR-146b; (G) OS for patients with miR-149; (H) OS for patients with miR-365a; (I) OS for patients with miR-130b; ccRCC, clear cell renal cell carcinoma; OS, overall survival.

Independence between the immune-related prognostic signature and other clinical factors

The clinical information of ccRCC patients including gender, age, clinical stage and TNM stage was analyzed in depth in order to evaluate the independent predictive ability of the immune-related miRNA model based on univariate and multivariate Cox regression analyses. Univariate analysis showed that clinical stage (P<0.001), TNM classification (P<0.001), and the immune-related miRNA model (P<0.001) were significantly related to OS ().

Table 2

Univariate and multivariate analyses of overall survival in ccRCC patients of TCGA

Variable	Univariate analysis			Multivariate analysis
	HR (95% CI)	P value		HR (95% CI)	P value
Age	0.7 (0.45–1.1)	0.11		0.59 (0.36–0.96)	0.034
Gender	1.1 (0.7–1.7)	0.67		1.0 (0.6–1.66)	1.000
Clinical stage	6.8 (3.98–11.6)	2.5E-12		15.16 (1.08–211.42)	0.043
T classification	18.9 (7.28–48.9)	1.5E-09		2.92 (0.74–11.46)	0.125
M classification	4.4 (2.71–7.1)	2.3E-09		0.23 (0.02–2.77)	0.250
N classification	4.8 (2.13–11.0)	1.6E-04		2.49 (0.84–7.41)	0.101
Prognostic model	2.27 (1.45–3.57)	4.2E-04		1.79 (1.08–2.95)	0.024

ccRCC, clear cell renal cell carcinoma; TCGA, The Cancer Genome Atlas; HR, hazard ratio; CI, confidence interval.

Functional enrichment analysis of miRNA-target genes associated with ccRCC

A total of 9 miRNA-target genes were validated by the miRTarBase dataset. According to the miRTarBase, there were 3,878 target genes for the 9 miRNAs. In addition, a total of 6,204 DEGs were obtained. There were 830 intersection genes between the 3,878 DEGs and the 6,204 validated target genes of the 9 miRNAs (). The molecular function (MF, ), cellular component (CC, ), and biological process in GO analysis (BP, ) were depicted i n “dotplot” from “clusterProfile” packages. BP analysis mainly enriched kidney development, renal system development, and kidney epithelium development. CC functional analysis mainly included collagen-containing extracellular matrix, external side of plasma membrane, and cell-cell junction. MF analysis mainly contained coenzyme binding, growth factor binding, and virus receptor activity. The results of KEGG enrichment analysis showed that it identified about 830 intersection genes associated with ccRCC, of which counts >10 were mainly enriched in Epstein-Barr virus infection, cell adhesion molecules (CAMs), Hippo signaling pathway, AGE–RAGE signaling pathway in diabetic complications, antigen processing and presentation pathway, and intestinal immune network for immunoglobin A (IgA) production (P<0.01 after P value adjustment; ). Additionally, in order to provide a detailed description of the complex relationship between KEGG pathway and target genes, the “pathway-gene” network was constructed, as shown in .

Figure 4

Functional enrichment analysis of miRNA target genes associated with ccRCC. (A) The 9-miRNA signature predicted 830 DEGs. (B) MF, (C) CC, (D) BP, (E) KEGG pathway analysis. (F) The “pathway-gene” network based on the KEGG signaling pathway. ccRCC, clear cell renal cell carcinoma; DEGs, differentially expressed genes; MF, molecular function; CC, cellular component; BP, biological process; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Discussion

KIRC is a highly malignant solid tumor that easily metastasizes to the rectum or lungs, and thus can seriously impinge on the survival of patients. Consequently, it has become increasing important to find novel prognostic markers that can aid in the treatment of patients. Therefore, understanding the pathogenesis of ccRCC is essential for the diagnosis and treatment of ccRCC. Assuming miRNA is the stabilizers of human tissues, miRNAs are considered to be novel biomarkers (26,27). Recently, a plethora of studies have indicated that specific miRNAs can play critical roles in the development of ccRCC (28-31). In this study, we systematically analyzed a miRNA dataset from TCGA database and screened DEMs using the “limma” package. We then identified those immune-related DEMs that are closely related to prognosis by using univariate and multivariate Cox regression analyses. Finally, a total of 59 immune-related DEMs were identified between 521 tumor samples and 71 normal samples. Some studies have shown that models based on multiple miRNAs have more stable performance than single miRNA models. Therefore, we established a 9-immune-related-miRNA signature (miR-21, miR-342, miR-149, miR-130b, miR-223, miR-365a, miR-9-1, miR-204, and miR-146b) that possessed a strong predictive diagnostic ability and prognostic significance for patients.

Several studies have shown that certain miRNAs are abnormally expressed during the occurrence and development of RCC (32-34). However, owing to the molecular heterogeneity of cancer, the reproducibility of the methods for identification of these miRNAs has limitations (35). In addition, the number of patient volunteers studied in real life is relatively small. TCGA database can provide gene expression data of solid tumor tissue and normal tissue. Making full use of these expression data is conducive to forming a better understanding of the underlying biological mechanism, which is necessary for improving the diagnosis of ccRCC and for creating a prognostic signature (36).

To identify an effective biomarker for cancer diagnosis and prognosis, we constructed a 9-immune-related miRNA prognostic signature from the ccRCC dataset. The 9-miRNA signature was an independent prognostic factor of ccRCC, with the survival rate of patients in the high-risk group being significantly lower than that of the low-risk group (). Out of the 9 miRNAs, miR-21 was regarded as an oncogene. Ample research indicates that miR-21 can promote cancer development by inhibiting tumor suppressor genes (32), and recently, Liang et al. reported that miRNA-21 could promote cell proliferation and differentiation, and prevent apoptosis by activating mTOR-STAT3 signaling pathway (37). Additionally, miR-21 was also found to be closely related to various tumors, including those of the breast (38), colon (39), bladder (40), liver (41), and kidney (42). Okato et al. revealed that the ectopic expression of miR-149 significantly inhibited cancer cell migration and invasion in ccRCC cells (43); meanwhile, Li et al. demonstrated that the downregulation of miR-130b inhibited cell migration, cell proliferation, and induced cell apoptosis of RCC, suggesting that miR-130b may be an effective biomarker for diagnosis and a therapeutic target for the treatment of RCC (44). Additionally, Xiao et al. discovered that the overexpression of miR-223-3p could enhance cell proliferation and metastasis in renal cancer cells (45), and Hildebrandt et al. showed that miR-9-1 may not only be involved in the occurrence and development of ccRCC but may also participate in the metastatic recurrence of KIRC (46). Moreover, Xiong et al. revealed that miR-204 can directly regulate the target gene, RAB22A, in inhibiting RCC proliferation and invasion (47). Finally, Yang’s findings suggested that HOXA11-AS could promote RCC growth and invasion by modulating the miR-146b-5p-MMP16 axis (48). Although miR-342 and miR-365a have not been reported in other studies, our results indicated that they showed significant differential expression in ccRCC samples, which suggests that these miRNAs may serve as significant biomarkers and therapeutic targets for treating the occurrence and development of ccRCC.

Although this study developed 9-miRNA immune-related prognostic model in ccRCC based on miRNA expression profiles in TCGA, there were still some limitations. First, the difference of the races, gender, age of ccRCC patient, tumor stages would lead to heterogeneity. Second, although our results were validated in testing set in TCGA dataset, the results were not verified by experiments such qPT-PCR in vitro. Further evidences are required to confirm our results by experiments on a great number of samples and clinical patients.

Conclusions

This study demonstrated that our 9-immune-related-miRNA signature could be a practical and reliable prognostic model for ccRCC. Furthermore, through evaluation using a testing set and a whole set of ccRCC patients, we showed that this prognostic model has the potential to be a valuable biomarker in ccRCC.

The article’s supplementary files as