Exploring the miRNA-mRNA Regulatory Network in Clear Cell Renal Cell Carcinomas by Next-Generation Sequencing Expression Profiles

Altered microRNA (miRNA) expression is a hallmark of many cancer types. The combined analysis of miRNA and messenger RNA (mRNA) expression profiles is crucial to identifying links between deregulated miRNAs and oncogenic pathways. Therefore, we investigated the small non-coding (snc) transcriptomes of nine clear cell renal cell carcinomas (ccRCCs) and adjacent normal tissues for alterations in miRNA expression using a publicly available small RNA-Sequencing (sRNA-Seq) raw-dataset. We constructed a network of deregulated miRNAs and a set of differentially expressed genes publicly available from an independent study to in silico determine miRNAs that contribute to clear cell renal cell carcinogenesis. From a total of 1,672 sncRNAs, 61 were differentially expressed across all ccRCC tissue samples. Several with known implications in ccRCC development, like the upregulated miR-21-5p, miR-142-5p, as well as the downregulated miR-106a-5p, miR-135a-5p, or miR-206. Additionally, novel promising candidates like miR-3065, which i.a. targets NRP2 and FLT1, were detected in this study. Interaction network analysis revealed pivotal roles for miR-106a-5p, whose loss might contribute to the upregulation of 49 target mRNAs, miR-135a-5p (32 targets), miR-206 (28 targets), miR-363-3p (22 targets), and miR-216b (13 targets). Among these targets are the angiogenesis, metastasis, and motility promoting oncogenes c-MET, VEGFA, NRP2, and FLT1, the latter two coding for VEGFA receptors.


Introduction
Approximately 3% of all cancers in adults occur in the kidney; therefore kidney cancer is one of the ten most frequently occurring cancers in western communities [1]. The primary histomorphologic type is clear cell renal cell carcinoma (ccRCC), which accounts for 80-85% of all kidney cancers followed by papillary renal cell carcinoma (pRCC) representing approximately 10% of renal cancers [2]. The gender and age distributions are similar between ccRCC and pRCC; however ccRCC has a worse prognosis with a 5-year survival of 77% [3,4]. This is attributed to an advanced tumor stage (38.1%) and visceral metastasis at diagnosis (14.5%) [5].
MicroRNAs (miRNAs) are a class of small non-coding RNAs (sncRNAs) that can repress gene expression through translational repression or messenger RNA (mRNA) deadenylation and decay by base pairing to partially complementary sites [6]. Deregulated miRNAs have i.a. been associated with formation of metastases, tumor progression, and tumor growth in ccRCC and anti-miRs were suggested as novel therapeutic strategies in the treatment of the disease [7][8][9].
Next-generation small RNA-Sequencing (sRNA-Seq) allows for unbiased quantitative and qualitative sncRNA profiling. When compared to miRNA array platforms, sRNA-Seq additionally enables the discovery of novel miRNAs as well as the detection of other differentially expressed 2 BioMed Research International sncRNAs like small nucleolar RNAs (snoRNAs) and transfer RNA (tRNA)-derived fragments that can mimic miRNA function [10].
Several microarray based studies have demonstrated 21 to 34 differentially expressed miRNAs between ccRCC and normal kidney tissue [11]. SRNA-Seq studies reported more than 100 differentially regulated miRNAs, some of which might serve as diagnostic and prognostic markers [12,13]. Nevertheless, these studies lack detailed information about miRNA targets and bioinformatical analysis is often only focused on miRNAs currently known to miRbase.
Here we used omiRas [14] to analyze a publicly available dataset (GEO: GSE24457) published by Zhou et al. [13], comprising twenty sRNA-Seq libraries of ten ccRCCs and ten adjacent control tissues from the same patient in order to identify sncRNAs with deregualted expression across all cases. After outlier detection with principle component analysis (PCA) samples of nine patients were used for downstream analysis.
We detected 61 sncRNAs as differentially expressed between the groups. Among these were several miRNAs without previous implication in kidney cancer development, like miR-3065-5p. Additionally, we detected seven snoRNAs and two tRNA derived fragments as differentially expressed between ccRCC and control tissues. We connected the deregulated miRNAs to biological pathways composed of differentially expressed genes under potential post-transcriptional control of these miRNAs. To do so, we utilized another publicly available mRNA-Sequencing (RNA-Seq) dataset (see methods). The "interaction network tool" of omiRas allows for the construction of interaction networks of miRNAs and mRNAs, interrogating the information from several miRNA-mRNA interaction databases. Therefore, we in silico assigned functions to significantly deregulated miRNAs and defined miRNAs implicated in the carciogenesis of ccRCC.
Among these is miR-206, which is significantly downregulated in ccRCC. Loss of miR-206 has been associated with hypoxia and under insufficient oxygen supply, angiogenesis is stimulated through upregulation of VEGF [15]. Our analysis revealed that miR-206 can regulate the expression of VEGF and several other genes involved in invasion, metastasis, and angiogenesis (MET, FN1, NRP1, ELMO1, and TAGLN2). This underlines that hypoxia induced loss of miR-206 expression is a critical process in clear cell renal cell carcinogenesis and maintenance.
Overall, our study demonstrates a promising strategy to identify driver miRNAs in cancer development that can afterwards undergo further functional testing.

Dataset Collection.
A publicly available sRNA-Seq expression dataset of ten ccRCC and matching normal renal tissue was downloaded from the Gene Expression Omnibus (GEO: GSE24457) database in SRA format. A list of 1,299 significantly upregulated and 1,194 downregulated genes identified in mRNA-Seq data of 65 ccRCC cases from the Cancer Genome Atlas (TCGA) by Wozniak et al. [16] was retrieved in XLS format.

Data Preprocessing.
Raw sequencing files were converted to FASTQ format and the 3 sequencing adapter (TCGTAT-GCCGTCTTCTGCTTGAAA) was removed from the reads with cutadapt [17]. Subsequently, low quality stretches below a SANGER quality score of 20 were additionally trimmed from each end of the reads (-q 20). Only reads with a minimum length of 15 bps after clipping were used for further analysis (-m 15).

Data Analysis
2.3.1. MiRNA Quantification, Outlier Detection, and Differential Expression Analysis. Preprocessed FASTQ-files were submitted to omiRas. Briefly, in omiRas, reads are summarized to UniTags. Singletons are removed from the data set and the remaining tags are mapped to the human genome (hg19) with bowtie [18] allowing at most two mismatches (controlled by option -v 2). Only the alignments in the best stratum are reported (if a read matches to six different genomic loci, two loci with no mismatch, and four loci with one mismatch, only the two alignments are reported) controlled with --strata and --best. Alignments for reads with more than 50 mapping locations are suppressed (-k 50). The mapped tags are annotated with the help of various models of coding and non-coding RNAs retrieved from the UCSC table browser [19]. Tags mapping to exonic regions of coding genes are excluded from further analysis. NcRNAs are quantified for each library independently. For tags mapping to multiple genomic loci the number of reads corresponding to the tag is divided by the number of mapping loci. To account for differences in sequencing depth, tag-counts are normalized (NEV, normalized expression value). Differential expression (corrected value (FDR) < 0.1) is detected with the DESeq bioconductor package [20] that takes biological and technical variance into account.
To reduce noise we introduced an outlier detection prior to differential expression analysis into the omiRas pipeline. The normalized counts are evaluated by PCA with R 3.0.2. The samples identified to be four or more standard deviations away from the mean on the first or second principal component are considered outliers and are removed from analysis.

Identification of miRNA Targets in ccRCC.
MRNA targets (as given in the XLS file of Wozniak et al. [16]) of differentially expressed miRNAs were identified with the "interactive network tool" of omiRas. An interaction between an miRNA and a coding gene is assumed to be valid if the following two criteria apply.
(a) Three or more of seven miRNA-mRNA interaction databases support the interaction.
(b) The expression of the miRNA/mRNA pair is inverse. The miRNA is significantly downregulated and the mRNA is upregulated or vice versa.

Identification of miRNAs Involved in the Deregulation of Genes from the Same Functional Category.
Up-and downregulated genes in ccRCC were mapped to functional Gene Ontology (GO) categories using DAVID Bioinformatics Resources 6.7 [21]. Genes within enriched (FDR < 0.05) categories were committed to the STRING database [22] to determine protein-protein interactions of their gene products. Additionally, miRNAs that might be causative for the deregulation of genes within the category were detected as described above.

Visualization.
PCA and hierarchical clustering of the differentially expressed miRNAs were performed and visualized with R 3.0.2. Networks of genes from the same GO category were visualized with Cytoscape [23]. Visualizations of annotation statistics for each library were taken from omiRas.

Results
The samples of patient P 10 were identified as outliers (see methods) and consequently tissues from the patient were not considered for downstream analysis. Therefore, differences in the sncRNAs between 9 patients' kidney tumors and adjacent normal renal tissue have been assessed.
Overall, 306,958,418 sequences were processed. On average, 7.5% of reads were discarded during adapter clipping and quality trimming. The proportion of sncRNAs in the libraries ranged from 63% to 75%, whereas miRNAs accounted for 37-74% of sncRNAs (see Table 1). The length distribution after adapter clipping revealed a clear peak at 22 base pairs (bps) across all samples (Figure 1(b) (see Figure S1 in Supplementary Material available online at http://dx.doi.org/10.1155/2014/948408)), which is typical of mammalian miRNAs [24]. Besides ncRNA encoding loci reads were annotated to intergenic regions ( Figure 1(c), S1) that may harbor novel miRNAs. The difference in miRNA proportions between libraries was attributed to varying degrees of snoRNAs, as well as tRNA derived fragments ( Figure 1(d), S1). 1,672 sncRNAs were expressed in at least one library (Supplementary Table T1); 41 of these were significantly downregulated and 20 were upregulated in ccRCC tissues. The most differentially up-/downregulated miRNAs are listed in Table 2. The PCA of the normalized expression values of these miRNAs indicated a clear separation of normal and cancer tissue samples via the first two principle components (Figure 2(a)).

Detection of miRNA Targets and Functional Enrichment.
To assert the influence of miRNAs on the gene expression pattern in ccRCC, we detected differentially expressed mRNAs with conserved seed sequences of inversely expressed miRNAs in their 3 UTR, supported by at least  sequences in ∼1/3 of all mRNAs (49 targets). Other miRNAs with more than ten targets are miR-135a-5p (32 targets), miR-206 (28 targets), miR-363-3p (22 targets), and miR-216b (13 targets). Some of the predicted interactions in the network have recently been validated in different cell types, like the upregulation of the c-Met oncogene (MET) [29], fibronectin (FN1) [30], or vascular endothelial growth factor A (VEGFA) [31] due to loss of miR-206 or the upregulation of E2F1 [32], CCND1 [33], and CDKN1A [34] due to loss of miR-106a-5p. A similar analysis was performed for upregulated miRNAs and downregulated mRNAs and the interactions are given in Supplementary Table T2.
In order to identify the influence of miRNAs on the metastatic potential of ccRCCs, we examined a subset of upregulated genes, enriched in GO category "cell motion" (FDR = 0.0002), and detected interactions of their gene products as well as the influence of miRNAs on the expression of the corresponding mRNAs. The interaction network is given in Figure 5. It is comprised of 36 mRNAs/gene products and six miRNAs. The network is centered around 12 highly connected gene/protein nodes, among these are transforming growth factor beta 1 (TGFB1), the vascular cell adhesion molecule-1 (VCAM1), the cell surface receptor CD44, and the intercellular adhesion molecule (ICAM1). The NRP2 and FLT1 mRNAs, both coding for VEGFA receptors, are potential targets of three and two miRNAs, respectively. Other targeted mRNAs from the network include cadherin 13 (CDH13), integrin A4 and A5 (ITGA4, ITGA5), chemokine CCL5, fibronectin (FN1), neuropilin-1 (NRP1), and the MET gene coding for the hepatocyte growth factor receptor.

Discussion
Our study links coding-and non-coding transcriptome data of normal and ccRCC tissue from two distinct studies. By the use of several miRNA-mRNA interaction databases available in omiRas we are able to provide new insights into the influence of aberrant miRNA expression on hundreds of deregulated genes. Hypoxic regions of tumors are often resistant to both radiotherapy and chemotherapy [35] and under insufficient oxygen supply angiogenesis is stimulated through upregulation of VEGF. VEGF-mediated angiogenesis is thought to play a critical role in tumor growth and metastasis [36]. We show that miR-206 is among the most significantly downregulated miRNAs in ccRCC and miR-206 loss has likewise been reported under hypoxic conditions [15]. As given in Figure 3, VEGFA is one of the predicted targets of miR-206. This prediction has gained support by the study of Zhang et al., who showed that miR-206 downregulation promotes proliferation and invasion of laryngeal cancer by regulating VEGF expression [37].
In line with downregulation of miR-206 and upregualtion of VEGFA under hypoxia, MET overexpression can be caused by hypoxia [38]. The overexpression of the cell surface  receptor Met tyrosine kinase MET is a hallmark of various types of cancers including pRCC and ccRCC [39]. In the normal kidney, ligand binding to MET mediates the activation of the MAPK, STAT, and a variety of other signaling pathways and promotes increased cell growth, scattering and motility, invasion, protection from apoptosis, branching morphogenesis, and angiogenesis [40]. MET overexpression in transgenic mice led to spontaneous development of hepatocellular carcinoma and and inactivation of the MET gene to tumor regression [41]. We identified MET as a proposed target of miR-206 and this prediction is supported by Yan and colleagues [29], who showed that miR-206 targets MET and inhibits rhabdomyosarcoma development.
Fibronectin (encoded by FN1) is a multifunctional extracellular glycoprotein and its increased expression is significantly associated with a higher probability of metastasis, poorer overall survival, and distant metastasis development [42]. FN1 expression was significantly increased in hypoxic mouse embryonic stem cells [43] and downregulation of miR-206 has been shown to induce FN1 upregulation in mouse lung tissues [30].
NRP1 encodes a receptor for VEGF and a block to NRP1 suppresses tumor growth due to decreased angiogenesis and cell proliferation [44]. In two different lung cancer cell lines, hypoxia regulated NRP1 expression differently. In the A549 AC cell line the expression increased, whereas a decreased expression was reported in SKMES-1 SCC [45]. Regulation of VEGF as well as VEGF receptor gene expression has been ascribed to transcription factor ETS1 [46]. Oikawa et al. showed that hypoxia induced ETS1 expression in human bladder cancer cell lines [47]. In several tumors, high ETS1 expression has been associated with decreased survival, angiogenesis, and poor prognosis [46].
Myristoylated alanine-rich C kinase substrate (MAR-CKS) controls mucus granule secretion by airway epithelial cells and directed migration of leukocytes, stem cells, Protein-protein interaction miRNA-mRNA interaction and fibroblasts. SiRNA knockdown of MARCKS expression in invasive lung cancer cell lines reduced migration [48]. In bladder cancer, MARCKS expression is induced under hypoxia [49]. Currently, an interaction between miR-206 and NRP1, ETS1, or MARCKS has not been experimentally verified but is predicted by our analysis. Therefore, these interactions represent an interesting target for experimental validation.
Other suggested targets of miR-106a by our analysis in ccRCC are i.a. E2F1, CDKN1A, and PREX. Upregulation of the transcription factor E2F1, a key regulator of proliferation and apoptosis, might be a driving force in the local and vascular infiltration of ccRCC [51]. In glioma, cell growth is inhibited by miR-106a due to post-transcriptional downregulation of E2F1 and accordingly downregulation of p53 [32]. Campbell and colleagues [52] showed that a common predicted target of miR-206 and miR-106a, PREX, is essential for metastasis formation in several cancer types by influencing physical migration processes through effects upon Rac1-driven motility.
One miRNA without previous implications in ccRCC is miR-3065-5p, an antisense miRNA to miR-338 [53]. Due to its novelty, currently no experimentally validated interactions are known for this miRNA. Based on our predictions it might target the previously described NRP2, FLT1, ETS1, and MARCKS and its loss can thereby contribute to angiogenesis in ccRCC.
Taken together the analysis underlines the role of hypoxia as a key factor in kidney tumor angiogenesis [54], which might i.a. be regulated by the loss of miR-206. We show that miR-206 has several targets that are upregulated under hypoxic conditions. A couple of these predicted interactions have already been experimentally validated, which highlights the validity of our bioinformatical in silico approach.

Conclusion
In this study we demonstrate how the combined analysis of miRNA and mRNA data with omiRas can explain differential gene expression signatures via loss/gain of posttranscriptional control by deregulated miRNAs. Without the integration of coding gene expression, differentially expressed miRNAs can only serve as biomarkers with nonspecific function. We assign roles to miRNAs without laborious functional testing by concentrating on the "low hanging fruit, " namely, interactions between miRNAs and mRNAs with increased likelihood of interaction due to the support of several databases. These promising candidates can afterwards undergo further functional testing. The validity of this strategy is underlined by the fact that several predictions from our study have recently been validated in cell lines. This makes omiRas a beneficial tool in cancer research. Furthermore, omiRas is not only limited to cancer studies, but is a useful tool to integrate coding gene expression profiles into miRNA analysis for any dataset that compares two different biological conditions.