Candidate Gene and MicroRNA Biomarkers of Metastatic Renal Cell Carcinoma (MRCC) Response to Sunitinib


 Backgrounds: The incidence of renal cancer is relatively insidious, and some patients have been metastatic renal cancer at the initial visit. Sunitinib is the first-line systemic therapy for patients with metastatic renal cell carcinoma, however, there is scant analysis of its effect on genes and microRNAs.Methods: In this study, 8 differentially expressed microRNAs and 112 differentially expressed genes were designated by analyzing mRNA and microRNA data sets and weighted correlation network analysis (WGCNA).Results: NIPSNAP1 gene showed the most co-expression with other genes. Through the intersection of the microRNA target gene with our differentially expressed genes, we got 26 genes. KEGG and GO analysis showed that these genes were predominantly concentrated in Pathways in cancer, Sphingolipid metabolism and Glycosaminoglycan degradation. After we set the 26 genes and gene of WGCNA do intersection, received six genes, respectively is NIPSNAP1, SDC4, TBC1D9, NEU1, STK40 and PLAUR. Conclusion: Through subsequent cell, molecular and flow cytometry experiments, we found the PLAUR would play a crucial role in renal cell carcinoma (RCC) resistant to sunitinib, which will be available for new ideas to forecast sunitinib resistance and reverse sunitinib resistance.


Introduction
Metastatic Renal Cell Carcinoma (mRCC) is the second only to bladder cancer in the incidence of urinary tract tumors in China, accounting for about 3% of all malignant tumors 1 .Due to the changes in the living environment and the improvement of diagnostic methods, the incidence of exogenous renal cell carcinoma has increased signi cantly in recent years, and it has been one of the most important killers threatening human health 2 . The incidence of renal cancer is relatively occult, with 20% to 30% of patients having metastatic renal cancer at the rst visit 3 , and more than 20% of patients with localized renal cancer will have metastasis after surgery 4 .The prognosis of metastatic renal cancer is poor, with a median survival of only 13 months and a 5-year survival rate of less than 10% 5 . However, due to the effects of occult incidence and limited diagnostic methods of renal cell carcinoma, most of the patients have reached the middle and advanced stage of renal cell carcinoma with varying degrees of metastasis, and are prone to metastasis and recurrence after surgery, resulting in poor prognosis 6 . Cytokine -based nonspeci c immunotherapy or combined tumor -reducing nephrectomy was the de nitive treatment for metastatic renal cancer 7 . With the development of sequencing technology and bioinformatics, diagnostic and therapeutic potential of miRNA has been proved 8 . For more than a decade, microRNA and mRNA target gene therapy has become the most promising treatment for metastatic renal cancer in recent years 9 .
MicroRNAs are a class of endogenous non-coding genes in the cells of animals and plants, with a size of about 21-25 bases, which play an important role in the regulation of gene expression 10 . Now studies have shown that microRNAs play a regulatory role in tumor development and apoptosis 11 . With further research, microRNAs are expected to become targets and molecular markers for tumor diagnosis, treatment and prognosis 12 . Celia Prior et al. identi ed the role of mir-942 in drug resistance from a comparison of patients with signi cant sunitinib e cacy and resistance 13 . Incorvaia Lorena found that a subset of 59 microRNAs involved in mRCC, such as miR-155, miR-22, miR-24, miR-484, miR-335 and miR-492 14 . Madrid et al. showed that INF2, Cdc42 and MAL2 could orderly form a pathway regulating endocytosis and lumen formation 15 . These studies suggest that INF2 may be involved in the development and metastasis of renal cancer through its speci c intracellular functions. Although there are many studies on MRCC related genes, micro and mRNA studies are still lacking.
Sunitinib is a novel drug that selectively targets a variety of receptor tyrosine kinases. Inhibition of receptor tyrosine kinases is thought to "starve" tumors and have the ability to simultaneously kill tumor cells by blocking the blood and nutrient supply needed for tumor growth. Sunitinib combines antiangiogenesis, which stops the blood supply to tumor cells, with anti-tumor mechanisms that directly attack tumor cells. Sunitinib has also been shown to be active against vegfr, platelet-derived growth factor receptor, stem cell factor receptor, glial cell-derived neurotrophic factor, and fms-like tyrosine kinase-3 16 . Sunitinib, as the standard drug for rst-line metastatic renal cell carcinoma, has a better effect than ifn-alfa and improves the production period of patients 17 . First-line observation also found that patients with different metastatic renal cell carcinoma had different tolerance to sunitinib, regardless of age or metastatic site 18 . However, the target genes and micrornas of sunitinib have not been thoroughly studied. The effect, especially in late-stage metastatic renal cell carcinoma, remains mysterious.
In this study, existing microarray expression pro le data of metastatic and non-metastatic renal clear cell carcinoma tumor samples were collected from public databases, and bioinformatics methods were used to screen speci c differentially expressed genes related to tumor metastasis of renal clear cell carcinoma, and functional annotation and biological function analysis of these genes were performed. At the same time, by analyzing the interactions between proteins encoded by differentially expressed genes, the protein-protein interaction network and the correlation analysis with related small-molecule drugs were constructed to provide certain data support for the clinical diagnosis and drug control treatment of renal clear cell carcinoma.

Identi cation of de antly expression gene and microRNA
In order to nd the differential expression of genes and microRNAs associated with sunitinib treatment, the GSE37766 and GSE65615 datasets were downloaded from NCBI GEO database. Because the mRNA dataset samples were repeated too much, the mRNA differential gene screening criteria were set to |logFC|>0.25, adjust. P <0.05.The results (Figure 1a) showed that there were 112 differentially expressed genes in sunitinib-treated patients compared with those who did not use sunitinib , including 84 upregulated genes and 28 down-regulated genes. In the same way, a total of 8 microRNAs (Figure 1b) with different expressions were found in the microRNA data set (|logFC|>1, adjust.p<0.05) consisting of 7 upregulated and 1 down-regulated micoRNA. The 8 differently expressed microRNAs were hsa-miR-30a, hsa-miR-30e, hsa-miR-126, hsa-miR-589, hsa-miR-938, hsa-miR-148b, hsa-miR-708 and hsa-miR-18b.

Weighted correlation network analysis of the DEGs associated with sunitinib
For the identi cation of sunitinib associated modules and genes, WGCNA was performed on the 112 DEGs. First, we constructed the evolutionary tree of the sample based on the expression data. The samples with excessive missing values were deleted and then analyzed. Next, we lter the power values. Finally, we chose 5 as the value of power. Next, we examined the association between the characteristic values of the network module and metastatic renal cell carcinoma. We got two modules, grey and turquoise. The pine green module contained the most differentially expressed genes, which were selected for further analysis. We used the differentially expressed genes contained in the pine green module to construct the co-expression network, as shown in the gure 1c, which contained 54 edges and 26 points, among which the NIPSNAP1 gene showed the most co-expression characteristics with other genes. GWCNA analysis narrowed the range of target mRNAs and removed the less relevant genes.

KEGG analysis and GO annotation
To investigate the biological function of the 26 differentially expressed genes based on the wgcna coexpression network. Gene Ontology (GO) functional and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed. GO results (Figure 2a, Table1) showed that in level 2 these genes mainly play role in binding, catalytic activity and molecular transducer activity. In level 3 the functions of these proteins are mainly concentrated in protein binding (20 genes), transferase activity (10 genes) and ion binding (13 genes). For the molecular function part, most genes were enriched in binding, protein binding and catalytic activity which appears in almost all GO analyses. Some interesting molecular functions were also found such as ion binding, cation binding and metal ion binding. For cellular component, most genes belonged to intracellular membrane-bounded organelle, cytoplasmic part and membrane-bounded vesicle. For the biological process, 16 genes were associated with metabolic process and 15 genes were associated with biological regulation. The results of KEGG (Figure 2b) showed that the main enrichment pathways of the 26 genes are as follows Pathways in cancer (PTGER3, LAMA3 and TXNRD), Sphingolipid metabolism (NEU1 and GLB1) and Glycosaminoglycan degradation (hya and GLB1).

MicroRNA target gene prediction
Using miRTarBase, the microRNA-target interacted gene which was experimentally validated were selected based on the DEG microRNAs. Totally, 162 genes were associated with differentially expressed microRNAs. To narrow the range, we cross-referenced the differentially expressed genes with all mRNAtarget genes ( Table 2). The results showed that a total of 26 (0.9%) genes were present simultaneously in both gene sets such as NIPSNAP1, CDK6, PPARA, FGFR1 and CDR1. After that, the 26 genes were intersected with the genes in WGCNA. As shown in the gure 2c, a total of 6 genes were found in the WGCNA and the differential microRNA target gene. These genes are NIPSNAP1, SDC4, TBC1D9, NEU1, STK40 and PLAUR.

PCA and PCoA
To test the representativeness of the differential genes, principal component analysis (PCA) and principal co-ordinates analysis (PCoA) were nished based on the 112 differently expressed genes. PCA is a dimension reduction algorithm commonly used in data mining. PCA uses the idea of dimensionality reduction to reduce the dimensionality of high-dimensional variables in the principle of minimizing data loss, that is, to nd a few comprehensive indexes among numerous variables. PCoA is a visualization method to study the similarity or difference of data. After sorting through a series of eigenvalues and eigenvectors, the principal ordinates are selected. The PCA result (Figure 2d) showed that two groups with (red circle) or without (green circle) sunitinib clearly separated from each other. The rst principal component and the second principal component explained 31.2% and 7.9% of the total variation, respectively. PCoA results ( Figure 2e) show similar results with PCA results that the rst primary coordinate and the second primary coordinate explained the 34.04% and 11.98% of the total variation, respectively. Combining the results of PCA and PCoA, the 112 differently expressed genes fully re ect the therapeutic effect of sunitinib on metastatic renal cell carcinoma.

The upregulated six genes in sunitinib resistant RCC cells
In order to con rm the NIPSNAP1, SDC4, TBC1D9, NEU1, STK40 and PLAUR expression, we grafted sunitinib resistant 786O (SR786O) cells and 786O cells (sunitinib sensitive) into nude mice forming tumor, and the qRT-PCR experiments were performed and the results showed in the gure 3. It indicated that in the six hub genes, only PLAUR were obviously high expressed in sunitinib resistant RCC tissues compared with sensitive RCC tissues, which were paritially in accordance with the previous data acquired from the bioinformatics analysis.

PLAUR regulates sunitinib resistance in RCC cells
To further demonstrate the function of PLAUR in sunitinib-resistant RCC cells, we steady knocked down PLAUR by transfecting lentiviral-mediated short hairpin RNA (shRNA) in SR786O cells ( gure 5a). We cultivated SR786O cells, which were transfected with negative shRNA or PLAUR shRNA-1/-2 A morphological detection of the incubated cells with 2µmol sunitinib for 60 hours, which manifested that PLAUR knocked down cells changed signi cantly from fusiform shape to round shape, and easier to oat compared with negative control ( gure 5c). In order to con rm this phenomenon, the CCK8 experiment was conducted and the result showed that both shPLAUR-1 and shPLAUR-2 remarkably suppressed the proliferation of SR786O cells in contrast to the negative control upon sunitinib administration at indicated concentrations for 60 hours ( gure 5b). To verify the therapeutic action of PLAUR knockdown on sunitinib resistance in SR786O cells, the analysis of SR786O cells by ow cytometry and annexin V-FITC/PI double dyeing experiment were conducted. The proportion of apoptosis in early and late were obviously elevated in PLAUR downregultaed SR786O cells in contrast to the negative control when treated with 2µM sunitinib ( gure 6 a and b). In addition, the cell cycle experiment were performed to detect outcome of the PLAUR knockdown in sunitinib resistant RCC cells. The result suggested the percentage of SR786O cells in G2 phase clearly increased relative to negative control ( gure 6 c and d).
So the apoptosis signaling antibody kit were used in the western blot experiment to test the levels of apoptosis relevant protein ( gure 7a). The result suggested that the expression of cleaved caspase-3, cleaved caspase-9 and cleaved PARP were remarkably upregulated in SR786O cells transfected with shPLAUR-1 and shPLAUR-2 as opposed to negative control ( gure 7 b, c, d and e).

Discussion
A growing number of researchers are concerned about the role of microRNAs in cancer development because they regulate gene expression and ultimately alter cell survival cycles. For instance miR-18b, hsa-miR-30a hsa-miR-30e which were included in mir-17-92 cluster were dysregulated in ccRCC in a study by Chow, Tsz-fung F. et al 33 . However, information on microRNAs as markers of sunitinib tolerance to MRCC is limited. We found eight differentially expressed microRNAs by comparing sunitinib-sensitive patients with drug resistance. They were hsa-mir-30a, hsa-mir-30e, hsa-mir-126, hsa-mir-589, hsa-mir-938, hsa-mir-148b, hsa-mir-708 and hsa-mir-18b, respectively. This is partly consistent with a previous study which found that hsa-miR-18 involved in resistance to sunitinib in Renal Cell Carcinoma Cells 34 . Gámez-Pozo et al also reported that twenty-eight microRNAs of the 287 include the has-miR-126 were related to resistance to rst-line sunitinib in advanced renal cell carcinoma patients 35 . However, in another study, candidate MicroRNA biomarkers in response to sunitinib therapy in metastatic renal cell carcinoma were not the same as in this study 4 (miR-155, miR-484, miR-221, miR-222, miR-425, miR-133, miR-410, miR-141, miR-628 and miR-942).To sum up, the differential expression of microRNA in MRCC patients to some extent re ects the difference in response to sunitinib drug treatment. In any case, more detailed experiments are needed to demonstrate the role of mRNAs in sunitinib resistance.
Sunitinib resistance is a challenge in advanced renal cell carcinoma. Therefore, it is necessary to understand the underlying mechanisms of sunitinib resistance. We compared patients with or without different responses to sunitinib and mined 26 potential hub genes through WGCNA analysis. These 26 hub genes were then intersected with microRNA target genes to obtain 6 genes in total. These six genes are NIPSNAP1, SDC4, TBC1D9, NEU1, STK40, PLAUR. Based on previous studies NIPSNAP1 belong to a highly conservative family of proteins of unknown function which is mainly associated with cognitive de cits 36 , pain transmission 37 , and neurological dysfunction 38 . Studies have shown that the NIPSNAP1 gene is inactivated during malignant metastasis of prostatic epithelial cells 39 .NIPSNAP1 was identi ed as a tumor suppressor in prostate cancer in a previous study identi ed as a tumor suppressor in prostate cancer 40 . We hypothesized that sunitinib might activate the NIPSNAP1 gene in some people to activate anti-tumor immune protein activity. Much work remains to be done. Syndecan-4 (SDC4) SDC4 plays an important role in the development and metastasis of renal cell carcinoma 11 . The decrease of tumor metastasis and necrosis will lead to the upregulation of SDC4 41,43 . In ros1 positive lung cancer, an activating KIT mutation of SDC4 induces crizotinib resistance 42 . We suspect that sunitinib may also cause mutations in SDC4 that lead to changes in resistance.
MicroRNA-125b was reported to promote invasion and metastasis of gastric cancer by targeting NEU1 44 . No signi cant changes in microrna-125b were found in our results, and sunitinib may not affect microrna-125b. Plasminogen activator receptor (PLAUR) has long been associated with human colon cancer, and abnormalities in the foxm1-PLAUR signaling pathway play a key role in the progression and metastasis of human colon cancer 45 . The overexpression of PLAUR and SERPINE1 has been correlated with poor survival outcomes and a high-grade tumor in RCC 46 . Our results also implied the association between PLAUR and metastatic renal cell carcinoma. Only in our results, the microRNAs related to PLAUR were hsa-mir-622 and hsa-mir-335 which had no signi cant difference in expression.
We found PLAUR were obviously upregulated in sunitinib resistant RCC cells in contrast to sunitinib sensitive RCC cell tissues. While, downregulation of PLAUR make the sunitinib resistant RCC cells acquired better response to sunitinib compared to the other ve genes in the loss-of-function analysis. Previous studies have reported the PLAUR induces non-small-cell lung cancer acquired ge tinibresistance through EGFR/p-AKT/survivin signaling pathway 47 . In our study, knockdown of PLAUR in sunitinib resistant RCC cells promoted apoptosis in the sunitinib treatment and upregulated the expression of cleaved caspase-3, caspase-9 and PARP in contrast to the negative control.

Conclusion
A total of 112 differentially expressed genes and 8 microRNAs were found in two-group (with or without sunitinib) comparisons of patients with metastatic renal cell carcinoma. Weighted correlation network analysis (WGCNA) results showed that NIPSNAP1 gene showed the most co-expression with other genes. Through the intersection of the microRNA target gene with our differentially expressed genes, we got 26 genes. KEGG and GO analysis showed that these genes were mainly concentrated in Pathways in cancer, Sphingolipid metabolism and Glycosaminoglycan degradation. After we put the 26 genes and gene of WGCNA do intersection, received six genes, respectively is NIPSNAP1, SDC4, TBC1D9, NEU1, STK40 and PLAUR. The PLAUR may induce renal cell carcinoma resistance to sunitinib. Postoperative specimens and serum PLAUR levels may serve as a biomarker to forecast the RCC patients' response to sunitinib, and inhibiting PLAUR may help sunitinib resistant RCC patients acquired better response to sunitinib.
The mRNA and miRNA combined datasets for myocardial infarction analysis were extracted from the Series Matrix File, respectively. The probe ID was converted into gene symbol or miRNA ID through the platform annotation information table. The same gene symbol or miRNA ID was incorporated. Expression matrixes of mRNA and miRNA samples were screened out according to ID. Other data sets are preprocessed in the same way.

Identi cation of DEGs
The identi cation of differentially expressed genes (DEGs) were performed with R software 3.4.3 20, ,edgeR (version 3.5.2) 21 http://bioconductor.org/packages/release/bioc/html/edgeR.html . The t-test method was used to test the difference between two groups. After t-test, all the genes obtained corresponding p-values. In order to make the results more reliable, we conducted multiple tests to correct the p-values through the Benjamin and Hochberg method. Finally, the threshold of DEGs was set as |logFC|>0.25, adjust.p<0.05. The obtained differential gene counts values were used for volcano plot and heat map.

Functional analysis of DEGs
In order to further clarify the functional annotation of DEGs, The common enrichment analysis tools AGRIgo 22 (http://systemsbiology.cau.edu.cn/agriGOv2/index.php) and DAVID 23 version 6.8 https://david.ncifcrf.gov/ were used to GO enrichment analysis of Gene Ontolog BP function and explore KEGG pathways analysis, respectively. The cutoff criterion was P-Value<0.05 and count>1. Protein sequences were extracted from the human genome protein bank based on the gene symbol using our own Perl scripts. We performed separate notes of DEGs using kaas 24 (https://www.genome.jp/tools/kaas/).
Construction of the WGCNA co-expression network R 3.5.0 package WGCNA 25 (1.6.6) was used to build the co-expression network. First of all, we need to process our data, the screening criteria were as follows: (a) retain 75% of the genes with higher absolute deviation, (b) the MAD value was greater than 0.01. Construct a sample tree to check whether or not outliers exist and determine the soft threshold (power value). Construct a hierarchical clustering tree. to show each module, and draw the correlation graph of each module. Visual gene network (TOM plot) was constructed, the nodes of the top100 were screened according to the wight value between nodes, and the network was exported and plotted by Cytoscape 26 . In detail, no outlier samples were found in this study. The soft threshold (power value) was 6 for the following analysis.

MicroRNA target gene prediction
MicroRNA target genes that were supported by strong experimental evidences (Reporter assay or Western blot) were downloaded from mirTarbase 27 (http://mirtarbase.mbc.nctu.edu.tw/php/download.php). Inhouse perl scripts were used to extract the data information and remove the irrelevant data based on the differently expressed microRNAs.

Principal component analysis (PCA) and principal co-ordinates analysis (PCoA)
An in-house perl scripts were used to obtain the expression data of the differently expressed gene and the grouping information. Then R (version 3.5.0) packages ggplot2 28 , plyr 29 , scales 30 , and ggbiplot 31 were applied together to do the principal component analysis. The X-axis represents the percentage of the rst principal component to the total variation. The Y-axis represents the percentage of the second principal component in the total variation. Principal co-ordinates analysis was completed by R (version 3.5.0) ggplot2 and vegan 32 with the hellinger method based on the expression data and grouping information.

Cells and cell culture
The human renal cell carcinoma cell line 786O was purchased from the American Type Culture Collection (ATCC, Rockville, MD, USA). In order to cultivate sunitinib-resistant 786O cells, we successively incubated sunitinib-sensitive 786O cells to an initial dosage of sunitinib (2 µM) and consecutively incremental concentrations up to 5 µM. The sunitinib-sensitive 786O and -resistant 786O cells were maintained in RPMI 1640 medium, which were supplemented with 10% FBS, 100IU/ml penicillin and 100 µg/ml (Thermo Fisher Scienti c, MA, USA).
Cell counting kit-8 (CCK-8) experiment SR786O cells were seeded into 96-well plates at the approximately 1000 cells per well with 200 µl of culture medium. 24 hours late, 10 µl CCK8 reagents (Dojindo laboratories, Japan) were added into each well, and the plates were incubated for 2 hours at 37℃. Finally the optical density were detected at the wavelength of 450 nm.
Flow cytometry of cell apoptosis Cell apoptosis was evaluated using the Annexin V-FITC/PI Apoptosis Detection Kit according to the explanatory memorandom, and all steps were carried out on ice to keep a 4 ℃ temperature. The transfected SR786O were treated with 5 µM sunitinib, washed with PBS and re-suspended in FITC Annexin V binding buffer, followed by incubated with Annexin V-FITC and stained with PI (Thermo Fisher Scienti c). In the end the proportion of apoptotic cells was assessed through ow cytometry.

Flow cytometry cell cycle experiment
The transfected SR786O cells were seeded in 6-well plates (5×10 5 cells/well) and were incubated with 5µM sunitinib. Subsequently, the SR786O cells were washed with PBS and suspended in DNA staining solution added with 5µl of permeabilization solution for half an hour at normal temperature. In the end, the proportion of cells in the cell cycle phase was evaluated on the SONY SH800 Cell Sorter (Tokyo, Japan).

Quantitative reverse transcription PCR
Total RNA was isolated using TRIzol Reagent (Thermo Fisher Scienti c, MA, USA), and cDNA was synthesized using reverse transcription kit (Tadara Biotechnology Co., Dalian, China). Afterwards, qRT-PCR was carried out using the ABI StepOnePlus system (Thermo Fisher Scienti c, Inc.). Relative expression levels was computed through -2 △△Ct method, and the β-actin was served as internal control.
The mRNA and miRNA combined datasets for myocardial infarction analysis were extracted from the Series Matrix File, respectively. The probe ID was converted into gene symbol or miRNA ID through the platform annotation information table. The same gene symbol or miRNA ID was incorporated. Expression matrixes of mRNA and miRNA samples were screened out according to ID. Other data sets are preprocessed in the same way.

Identi cation of DEGs
The identi cation of differentially expressed genes (DEGs) were performed with R software 3.4.3 20, ,edgeR (version 3.5.2) 21 http://bioconductor.org/packages/release/bioc/html/edgeR.html . The t-test method was used to test the difference between two groups. After t-test, all the genes obtained corresponding p-values. In order to make the results more reliable, we conducted multiple tests to correct the p-values through the Benjamin and Hochberg method. Finally, the threshold of DEGs was set as |logFC|>0.25, adjust.p<0.05. The obtained differential gene counts values were used for volcano plot and heat map.

Functional analysis of DEGs
In order to further clarify the functional annotation of DEGs, The common enrichment analysis tools AGRIgo 22 (http://systemsbiology.cau.edu.cn/agriGOv2/index.php) and DAVID 23 version 6.8 https://david.ncifcrf.gov/ were used to GO enrichment analysis of Gene Ontolog BP function and explore KEGG pathways analysis, respectively. The cutoff criterion was P-Value<0.05 and count>1. Protein sequences were extracted from the human genome protein bank based on the gene symbol using our own Perl scripts. We performed separate notes of DEGs using kaas 24 (https://www.genome.jp/tools/kaas/). Construction of the WGCNA co-expression network R 3.5.0 package WGCNA 25 (1.6.6) was used to build the co-expression network. First of all, we need to process our data, the screening criteria were as follows: (a) retain 75% of the genes with higher absolute deviation, (b) the MAD value was greater than 0.01. Construct a sample tree to check whether or not outliers exist and determine the soft threshold (power value). Construct a hierarchical clustering tree. to show each module, and draw the correlation graph of each module. Visual gene network (TOM plot) was constructed, the nodes of the top100 were screened according to the wight value between nodes, and the network was exported and plotted by Cytoscape 26 . In detail, no outlier samples were found in this study. The soft threshold (power value) was 6 for the following analysis.

MicroRNA target gene prediction
MicroRNA target genes that were supported by strong experimental evidences (Reporter assay or Western blot) were downloaded from mirTarbase 27 (http://mirtarbase.mbc.nctu.edu.tw/php/download.php). In-house perl scripts were used to extract the data information and remove the irrelevant data based on the differently expressed microRNAs.
Principal component analysis (PCA) and principal co-ordinates analysis (PCoA) An in-house perl scripts were used to obtain the expression data of the differently expressed gene and the grouping information. Then R (version 3.5.0) packages ggplot2 28 , plyr 29 , scales 30 , and ggbiplot 31 were applied together to do the principal component analysis. The X-axis represents the percentage of the rst principal component to the total variation. The Y-axis represents the percentage of the second principal component in the total variation. Principal co-ordinates analysis was completed by R (version 3.5.0) ggplot2 and vegan 32 with the hellinger method based on the expression data and grouping information. The transfected SR786O cells were seeded in 6-well plates (5×10 5 cells/well) and were incubated with 5µM sunitinib. Subsequently, the SR786O cells were washed with PBS and suspended in DNA staining solution added with 5µl of permeabilization solution for half an hour at normal temperature. In the end, the proportion of cells in the cell cycle phase was evaluated on the SONY SH800 Cell Sorter (Tokyo, Japan).
Quantitative reverse transcription PCR Total RNA was isolated using TRIzol Reagent (Thermo Fisher Scienti c, MA, USA), and cDNA was synthesized using reverse transcription kit (Tadara Biotechnology Co., Dalian, China). Afterwards, qRT-PCR was carried out using the ABI StepOnePlus system (Thermo Fisher Scienti c, Inc.). Relative expression levels was computed through -2 △△Ct method, and the β-actin was served as internal control.
Western blot experiment RIPA kit (Beyotime Biotechnology, China) and phosphorylation protease inhibitor were used to abstract the protein of cells. The proteins were separated by gel electrophoresis (SDS-PAGE) and then transferred to PVDF membranes (Millipore, Bedford, MA, USA). The membranes were incubated with anti-PLAUR (cat.   Figure 1 Volcano gure of the differently expressed mRNA (a) and miRNA (b) (|logFC|>1,adjust. P <0.05). The red dots represent patients with sunitinib treatment, meanwhile, the green or blue dots represent patients without sunitinib treatment. Weighted correlation network of the differently expressed genes in turquoise module. X-axis represented the weighting power and y-axis represented the quadratic correlation index derived from log (k) and log (P(k)) of the corresponding network (c). differently expressed genes. X-axis and Y-axis of principal co-ordinates analysis represent the rst principal co-ordinates and the second co-ordinates component to the total variation, respectively (e).

Figure 3
The six hub genes were veri ed by qRT-PCR experiments in sunitinib resistant and sensitive mRCC tissues. The PLAUR expression levels were consistent with the previous data obtained from the bioinformatics analysis