The mRNA-seq data, lncRNA-seq data, and clinical information of ccRCC in The Cancer Genome Atlas Kidney Clear Cell Carcinoma database (TCGA-KIRC) were downloaded from the Xena Functional Genomics Explorer (Xiao et al., 2020). All 1793 immune-related genes were collected from ImmPort (https://www.immport.org) (Bhattacharya et al., 2018). There were 2483 rows of data that were included in the file we downloaded from ImmPort; then, the duplicate genes were eliminated and 1793 genes were left. We screened the differentially expressed genes (DEGs) of KIRC using the “limma” package with parameters of adjusted p-value < 0.05 and log₂|FC| > 1. Then, the intersection of DEGs and 1793 immune-related genes were obtained using a Venn diagram.

The Database for Annotation, Visualization, and Integrated Discovery (DAVID) was used to obtain the functional enrichment results of 586 DEGs. We selected the top15 enrichment results of the biological process (BP), cellular component (CC), and molecular function (MF) to present. The top 20 signaling ways of Kyoto Encyclopedia of Genes and Genomes (KEGG) results were shown. p-value < 0.05 was considered significant. The protein-protein interaction (PPI) network of all DEGs was pictured by Cytoscape (v.3.6.1) using the “string” plugin. The confidence score was set as 0.40. Then, the PPI network was further analyzed by Cytoscape using the “MCODE” plugin to figure out the important modules. The top five modules were finally selected. Next, the “CytoHubba” plugin was conducted to screen the hub genes from all 586 DEGs with the criteria of degree >10.

Patients of ccRCC were divided into high-expression groups and low-expression groups according to the median of each gene expression. The “survival” package and R software were used to evaluate the prognostic value of the hub genes. GEPIA was utilized to verify the overall survival (OS) prognosis of the above hub genes (Tang et al., 2017). Then, multivariate cox proportional hazard regression analysis was conducted by SPSS (version 23.0). The value of p < 0.05 was significant. There were nine genes (AGER, HAMP, LAT, LTB4R, NR3C2, SEMA3D, SEMA3G, SLC11A1, and VAV3) finally screened after the above three steps. The Metascape database was used for GO enrichment analysis for hub genes (Zhou et al., 2019).

Starbase (Yang et al., 2011) and miRTarbase (Hsu et al., 2011) databases were used to predict miRNA reversely based on nine genes. The prediction results were merged and LAT was removed because there was no miRNA predicted. Cytoscape and plugin “CytoHubba” were conducted to obtain the top 15 miRNA according to the score of degree. The expression of miRNA was explored by UALCAN (Chandrashekar et al., 2017). The survival analysis of miRNA was identified by OncoLnc (http://www.oncolnc.org/). Next, LncBase (Paraskevopoulou et al., 2013) and Starbase databases were used to predict lncRNA based on miRNA, both of which had different expressions and were related to patients’ survival. There were 29 lncRNA in the intersection of prediction results, which were screened. Graphpad prism (v 8.2.1) was used to analyze the differential expression of lncRNA through an unpaired t-test, OncoLnc was used for survival analysis, and SPSS was used to conduct multivariate cox proportional hazard regression analysis. Finally, we constructed a new immune-related ceRNA network including NNT-AS1 (lncRNA), hsa-miR-186-5p (miRNA), NR3C2, and VAV3.

LinkedOmics (Vasaikar et al., 2018) was utilized to conduct the Gene set enrichment analysis (GSEA) of two genes in the ceRNA network (NR3C2 and VAV3). Connectivity Map (Cmap) (Lamb, 2007) was a common small-molecular drug-prediction database, so we predicted potential therapeutic chemicals for ccRCC through three genes (NR3C2, VAV3, and HAMP) associated with miR-186-5p. PubChem (Kim et al., 2018) was used to find their structure and other information.

In order to explore the relationship between immune-infiltration level and the two immune genes, we used CIBERSORT (Chen et al., 2018; Kong et al., 2020; Meng et al., 2021a) (https://cibersort.stanford.edu/) to obtain the immune infiltration level of TCGA-KIRC data. The algorithm was run using LM22 signature and 1,000 permutations. The 530 samples of TCGA-KIRC were divided into two groups according to the median of VAV3 or NR3C2. Then, R software and the “ggplot2” package were used to further present the immune cell infiltration level difference between two groups. The Wilcoxon test was used to calculate the p value. p value <0.05 was considered significantly.

Paired normal and cancer tissues of ccRCC were obtained from Department of Urology, Union Hospital, Tongji Medical College, Wuhan, China, with the approval of the Ethics Committee of Huazhong University of Science and Technology. According to the instructions, we performed RNA Isolation and Real-time PCR analysis with SYBR Green mix (Thermo, Massachusetts, USA). Primers of NNT-AS1 and hsa-miR-186-5p were obtained from RiboBio (Guangzhou, China). Primers of VAV3, NR3C2, and normalized gene GAPDH were obtained from Sangon Biotech (Shanghai):

GAPDH

Forward 5′‐GAG​TCA​ACG​GAT​TTG​GTC​GT‐3′

Reverse 5′‐GAC​AAG​CTT​CCC​GTT​CTC​AG‐3′

VAV3

Forward 5′‐AGA​GAA​ACG​GAC​CAA​TGG​ACT‐3′

Reverse 5′‐GGT​GGT​GTT​CCA​GAA​TAG​TTC​C‐3′

NR3C2

Forward 5′‐GAA​AGA​CGG​TGG​GGT​CAA​GTT‐3′

Reverse 5′‐ACC​GGA​AAC​ACA​GCT​TAC​GTT‐3′

We screened 2858 upregulated DEGs and 2354 downregulated DEGs in the TCGA-KIRC database, setting the parameters of p < 0.05 and log₂|FC| >1 (Figure 1A). Taking the intersection of DEGs and 1793 immune-related genes, 396 upregulated immune genes and 190 downregulated immune genes were finally determined and used for the following analysis (Figure 1B).

To determine the biological functions and potential signaling ways of DEGs, we used the DAVID database to perform GO and KEGG analysis of them. We selected the top 15 terms of GO and the top 20 terms of KEGG since there too many enrichment results. The results indicated that upregulated DEGs were enriched in biological processes (BP) of T-cell costimulation and interferon-gamma mediated signaling ways, while downregulated DEGs participated in cell proliferation and migration (Figures 2A,B). Cellular component (CC) enrichment results showed that both upregulated and downregulated DEGs were involved in the formation of various membrane and receptor complexes (Figures 2C,D). Molecular function (MF) results suggested that upregulated DEGs were associated with the tumor necrosis factor receptor, tumor necrosis factor-activated receptor activity, and so on (Figure 2E). Downregulated DEGs were related to transforming growth factor-beta receptor binding, fibroblast growth factor receptor binding, and so on (Figure 2F). KEGG enrichment indicated that both groups participate in various important signaling ways in ccRCC development (Figures 2G,H).

We used Cytoscape and the “string” plugin to construct the PPI network of all 586 DEGs. Then, “MCODE” plugin was utilized to divide the DEGs into several modules. We selected the top five modules for the following enrichment analysis. Genes in module 1 were associated with neutrophils, monocytes, lymphocyte chemotaxis, and positive regulation of angiogenesis (Figures 3A,B). Genes in module 2 were involved in cell adhesion molecules, positive regulation of B-cell proliferation, and natural killer cell-mediated cytotoxicity (Figures 3C,D). The TNF-signaling pathway and IL8 production were enriched in module 3 (Figures 3E,F). The NF-kappa B signaling pathway, Rap1 signaling pathway, and JAK-STAT signaling pathway were enriched in module 4 (Figures 3G,H). It was worth mentioning that the above three pathways were also enriched in the other four modules. DEGs in module 5 were related to the PI3-Akt signaling pathway and MAPK signaling pathway (Figures 3I,J).

First, we used Cytoscape and “CytoHubba” plugin to initially screen hub genes. There were 449 genes with degree >10 selected. Then, the OS rate of these genes was evaluated through the “survival” package. GEPIA was an online website by which we verified the prognostic value of the above genes. Multivariate cox proportional hazard regression analyses were further determined by the hub genes (Table 1). In addition, there were nine genes (AGER, HAMP, LAT, LTB4R, NR3C2, SEMA3D, SEMA3G, SLC11A1, and VAV3) finally selected. However, since later we found that there were no miRNA prediction results for LAT, we decided to remove it (Figures 4A–H). AGER, HAMP, LTB4R, and SLC11A1 were risk factors for ccRCC, while NR3C2, SEMA3D, SEMA3G, and VAV3 were protective factors for ccRCC. We obtained the GO enrichment results of the eight hub genes using the meta scale database. Interestingly, there were two clusters in which one (AGER, HAMP, SLC11A1, and VAV3) participated in macrophage activation (Figure 4I), and the eight hub genes were associated with cell proliferation and some immune processes (Figure 4J).

We predicted 673 miRNAs through Starbase and miRTarbase for all the eight hub genes (Figure 5). Then, Cytoscape and “cytoHubbe” plugin were used again to obtain the top 15 miRNAs associated with the hub genes (Figure 6A). Next, we explored the expression of miRNA using the UALCAN database and the prognostic value of them using the OncoLnc website. Then, hsa-miR-186-5p (alias, hsa-miR-186) was identified (Figures 6B,C). The related genes NR3C2 and VAV3 were negatively associated with hsa-miR-186-5p, which was consistent with the known mechanism of miRNA.

As mentioned before, we predicted 1417 lncRNAs using Lncbase and 157 lncRNAs using Starbase for hsa-miR-186-5p. There were 29 lncRNAs in the intersection (Figure 7A). We checked their expression and conducted survival analysis and multivariate cox proportional hazard regression analyses to determine the final lncRNA. In addition, according to the ceRNA hypothesis, lncRNA should be negatively related to miRNA. Finally, NNT-AS1 was the only one suitable. NNT-AS1 had a much lower expression level in ccRCC tissues than in normal tissues (Figure 7B). The low expression group showed poorer survival (Figure 7C). The results of cox regression were presented in the forest plot (Figure 7D).

As a result, a new ceRNA network (NNT-AS1—hsa-miR-186-5p—NR3C2/VAV3) was constructed by us. What is more, GSEA results of NR3C2 and VAV3 indicated that both genes were associated with the cell cycle, cytokine–cytokine receptor interactions, the fatty acid metabolism, and the citrate cycle (TCA cycle), which made great significance in the development of ccRCC (Figure 8A). To further verify the results from bioinformatics analysis, we performed qRT-PCR analysis on 30 pairs of clinical samples. As expected, the results showed that NNT-AS1, VAV3, and NR3C2 were remarkably downregulated in ccRCC compared to normal kidney tissues, while hsa-miR-186-5p was significantly highly expressed in ccRCC (Figures 8B–E). Moreover, the correlation analysis results further suggested that NNT-AS1 was negatively related to miR-186-5p (Figure 8F) but positively related to VAV3 and NR3C2 (Figures 8G,H). In addition, MiR-186-5p was negatively associated to VAV3 and NR3C2 (Figures 8I,J). In conclusion, the overall results were consistent with the hypothesis of the ceRNA network (NNT-AS1—hsa-miR-186-5p—NR3C2/VAV3).

Since VAV3 and NR3C2 are both immune-related genes, we decided to further explore how these two genes influenced the immune infiltration of the ccRCC microenvironment using the CIBERSORT database. The results showed that the proportions of CD4⁺ memory resting T-cells, monocytes, M1, M2, resting DCs, and resting Mast cells were higher in VAV3high than in VAV3low groups. On the contrary, plasma cells, CD8+T cells, and Tregs were much more in VAV3low than in VAV3high groups (Supplementary Figure S1). In the NR3C2high group, the proportions of naïve B cells, memory resting CD4+T cells, resting NK cells, monocytes, M2, resting DCs, and resting Mast cells were higher, while the fractions of plasma cells, CD8+T cells, memory activated CD4+T cells, Th cells, γδT cells, and M0 macrophages were lower (Supplementary Figure S2).

To figure out the association between clinical characteristics and the two immune-related genes, we analyzed the mRNA expression levels of VAV3 and NR3C2 in different clinical subgroups. The classification standard was based on the previous literature (Meng et al., 2021b; Lou et al., 2021). The results showed that the expression of VAV3 had no significance between the elderly group (age ≤ 60 years) and young group (age > 60 years) (Supplementary Figures S3A,B). Besides, VAV3 expressed higher in females than in males. However, the expression level of VAV3 decreased with the increase of T stage, N stage, M stage, TNM stage, and G grades (Supplementary Figures S3C–G). It may meant that VAV3 was not only downregulated in ccRCC but also decreased as the malignancy of tumors increased. Similarly, the expression level of NR3C2 was not correlated with age (Supplementary Figure S4A) and gender (Supplementary Figure S4B) but importantly negatively related to T stage, N stage, M stage, TNM stage, and tumor G grade (Supplementary Figures S4C–G).

Using the Cmap database and three hub genes associated with hsa-miR-186-5p, we predicted potential drugs for ccRCC. NR3C2 and VAV3 were downregulated in ccRCC, while HAMP was upregulated. Drugs with negative connectivity scores were considered therapeutic. We finally selected five small molecular drugs in which enrichment < −0.7 and p-value < 0.01 (Figure 9; Table 2). Their structure and other information were presented through PubChem, in which Prestwick-691 could not be found.