The MicroRNA Expression Signature of Bladder Cancer by Deep Sequencing: The Functional Significance of the miR-195/497 Cluster

Current genome-wide microRNA (miRNA) expression signature analysis using deep sequencing technologies can drive the discovery of novel cancer pathways regulated by oncogenic and/or tumor suppressive miRNAs. We determined the genome-wide miRNA expression signature in bladder cancer (BC) by deep sequencing technology. A total of ten small RNA libraries were sequenced (five BCs and five samples of histologically normal bladder epithelia (NBE)), and 13,190,619 to 18,559,060 clean small RNA reads were obtained. A total of 933 known miRNAs and 17 new miRNA candidates were detected in this analysis. Among the known miRNAs, a total of 60 miRNAs were significantly downregulated in BC compared with NBE. We also found that several miRNAs, such as miR-1/133a, miR-206/133b, let-7c/miR-99a, miR-143/145 and miR-195/497, were located close together at five distinct loci and constituted clustered miRNAs. Among these clustered miRNAs, we focused on the miR-195/497 cluster because this clustered miRNA had not been analyzed in BC. Transfection of mature miR-195 or miR-497 in two BC cell lines (BOY and T24) significantly inhibited cancer cell proliferation, migration and invasion, suggesting that the miR-195/497 cluster functioned as tumor suppressors in BC. Regarding the genes targeted by the miR-195/497 cluster, the TargetScan algorithm showed that 6,730 genes were putative miR-195/497 targets, and 113 significantly enriched signaling pathways were identified in this analysis. The “Pathways in cancer” category was the most enriched, involving 104 candidate target genes. Gene expression data revealed that 27 of 104 candidate target genes were actually upregulated in BC clinical specimens. Luciferase reporter assays and Western blotting demonstrated that BIRC5 and WNT7A were directly targeted by miR-195/497. In conclusion, aberrant expression of clustered miRNAs was identified by deep sequencing, and downregulation of miR-195/497 contributed to BC progression and metastasis. Tumor suppressive miRNA-mediated cancer pathways provide new insights into the potential mechanisms of BC oncogenesis.


Introduction
In developed countries, bladder cancer (BC) is the fifth most commonly diagnosed tumor and the second most common cause of death in patients with genitourinary tract malignancies. In the United States, there are more than 73,000 new cases and 14,000 deaths annually [1]. In Japan, the number of new BC patients was estimated at 17,461 in 2007 and the number of deaths was estimated at 7,008 in 2011 [2].
BCs can be classified into two categories, non-muscle-invasive tumors and muscle-invasive tumors. Although 70%-80% of patients are diagnosed with non-muscle-invasive tumors, high recurrence rates (50%-70%) are observed in these patients. Moreover, among recurrent cases, 15% of BCs progress to muscle-invasive disease [3]. The five-year survival rate for patients with non-muscle-invasive BC is close to 90%, whereas that of patients with muscle-invasive BC is approximately 60% [4]. Furthermore, nearly 80% of patients with lymph node metastases die within the first five years [5]. Since most clinical trials of chemotherapeutics for advanced BC have shown limited benefits, new prognostic markers and effective treatment strategies based on current cancer-genome analyses are necessary.
The discovery of non-coding RNA in the human genome was an important conceptual breakthrough in the post-genomic sequencing era [6]. Improved understanding of non-coding RNA is necessary for continued progress in cancer research. miRNAs constitute a class of small, non-coding RNA molecules, 19-22 nucleotides in length, that modulate gene expression. Regulation is achieved through imperfect pairing with target messenger RNAs (mRNAs) of protein-coding genes and transcriptional or post-transcriptional regulation of their expression [7]. Currently, 2,042 human mature miRNAs are registered at miRBase release 20.0 [http://microrna. sanger.ac.uk/]. A growing body of evidence indicates that miRNAs also contribute to the initiation, development and metastasis of various types of cancers. Many human cancers show aberrant expression of miRNAs that can function either as tumor suppressors or oncogenes [8]. Therefore, identification of aberrantly expressed miRNAs is the first step toward elucidating miRNA-mediated oncogenic pathways in human cancers.
The development of high-throughput, deep sequencing technology has rapidly uncovered novel information about miRNAs. Deep sequencing analysis seems to be superior to microarray-or PCR-based methods that are limited to known miRNAs and usually do not contain the full list of known miRNAs sequences. Deep sequencing analysis will become the gold standard method for comprehensive miRNA analysis in cancer genomics. miRNA expression signatures of BC obtained by deep sequencing technology have led to four recent publications [9][10][11][12] and this study constitutes the fifth report. A total of ten small RNA libraries were sequenced (five bladder carcinomas and five matched, histologically normal samples of urothelia), leading to 13,190,619 to 18,559,060 clean small RNA reads in this analysis. We detected 933 known miRNAs and 17 new miRNA candidates in our series of samples.
In this study, we hypothesized that the miR-195/497 cluster functioned as a tumor suppressor through targeting of several oncogenic genes involved in specific cancer-related pathways in BC. Based on this hypothesis, we performed gain-of-function studies in microRNAs transfectants and also attempted to identify novel miR-195/497 cluster-mediated molecular targets and pathways by in silico analysis. Understanding the role of miRNAs will provide an effective and promising strategy for miRNAassociated evidence-based therapeutics for the treatment of BC.

Sequencing and annotation of small RNAs
Ten small RNA libraries (from five BCs and five NBEs) were sequenced by using deep sequencing technology. The patients' characteristics are shown in Table 1. Initially, 13,905,124 to 18,938,856 raw sequence reads were produced for the libraries. After filtering and trimming of the reads for low quality and adaptors, from 13,190,619 to 18,559,060 clean small RNA reads were obtained ( Table 2). The distribution of nucleotide lengths of clean small RNA reads varied from 16 to 38 nucleotides in each library. The most abundant length was 22 nucleotides ( Figure S1). All of the high-quality clean reads larger than 18 nucleotides were mapped to the human genome and these genome-matched reads were divided into different categories of small RNAs according to their biogenesis and annotation (Table 3). Both BCs and NBEs contained multiple and heterogeneous small RNAs species that included miRNA, degradation fragments of non-coding RNAs (tRNA, rRNA, snRNA, snoRNA, etc.), genomic repeats, mRNAs, or unclassified RNAs. The most abundant RNA category from each library was miRNA. The sequence data from this study were deposited in DDBJ Sequence Read Archive (accession number; DRA001043).

Identification of downregulated miRNAs based on miRNA expression signatures of BC
We separated the high quality clean read sequences in each library into two categories, miRNAs or non-coding RNAs using NCBI GenBank data, Rfam data and miRBase. In this study, a total of 933 known miRNAs and 17 candidate new miRNAs were detected from high quality clean read sequences (Tables S1 and  S2). The expression levels of the new miRNA candidates were low and the functional roles were not assessed enough (Table S2). For these reasons, we excluded them from the analysis. However, the functional significance of these new miRNA candidates is important and they will be studied in the future. Here, we focused on the 933 known miRNAs and established the miRNA expression signature of BC by deep sequencing. Differentially expressed miRNAs were identified by using the edgeR program with a general linear model. A total of 60 downregulated miRNAs were selected for this analysis (Table 4).
Next, we mapped the downregulated miRNAs in the human genome and found several miRNAs that were close together and were termed clustered miRNAs, such as miR-1/133a, miR-206/ 133b, let-7c/miR-99a, miR-143/145 and miR-195/497 (Table 5). Because the role of the miR-195/497 cluster had not been reported in BC, we focused on it for further studies. As shown in Figure 1 (Figure 2A Predicted target genes and cancer pathways of the miR-195/497 cluster To explore the biological functions of the miR-195/497 cluster, we performed in silico analyses using two independent algorithms, TargetScan, and GeneCodis3, and public microarray expression data approved by GEO. The workflow for in silico analysis is shown in Figure 4. At first, the TargetScan database indicated that 6,730 genes were predicted targets of the miR-195/497 cluster, which shared the same seed sequence 'GCUGCU' in their 39-untranslated regions (39UTR).
Based upon the KEGG pathway classification, these genes were implicated in 113 pathways, and ''Pathways in cancer'' was the most significantly enriched (Table S3). A total of 104 genes were involved in this pathway. We searched the oncogenes targeted by the miR-195/497 cluster based on the hypothesis that putative oncogenes would be upregulated in clinical specimens due to downregulation of tumor suppressive miR-195/497 miRNAs. By using GEO database analysis, 27 genes were selected as miR-195/ 497-regulated oncogenic genes in BC ( Figure 5 and Table S4). To confirm this analysis, we selected the top two upregulated genes (BIRC5 and WNT7A), and validated whether these miRNAs were regulated by miR-195/497 (below).  The mRNA and protein expression levels of BIRC5 and WNT7A were downregulated in miR-195and miR-497-transfected BOY and T24 cells in comparison with control transfectants (Figure 6A, B).
We next performed luciferase reporter assays to determine whether BIRC5 and WNT7A mRNAs had functional target sites for miR-195 and miR-497 in BC. We used a vector encoding the 3-9UTRs of BIRC5 and WNT7A mRNAs, including putative target sites for miR-195 and miR-497 (positions 198-204 and 197-203, respectively; Table S5). We found that the luminescence intensity was significantly decreased in the presence of the target site of BRIC5 by miR-195and miR-497-transfectants ( Figure 7A). We obtained the same result with a luciferase reporter assay for the WNT7A gene ( Figure 7B). These data suggested that miR-195 and miR-497 bound directly to specific sites in the 39-UTR of both BRIC5 and WNT7A mRNAs.

Discussion
It is believed that miRNAs regulate the expression of approximately 30% of all genes in the human genome [29]. In normal cells, the interaction of miRNAs and target messenger RNAs (mRNAs) is tightly regulated, whereas this regulation is often missing in cancer cells. A growing body of evidence suggests that miRNAs are aberrantly expressed in many human cancers and that they play significant roles in the initiation, development, and metastasis of these cancers [30]. Therefore, identification of aberrantly expressed miRNAs is an important step in the analysis of human oncogenesis. The latest analytical method in genomics, deep sequencing technology, has the potential to identify novel tissue-specific miRNAs including those in human cancer tissues [31]. Deep sequencing technology is currently the gold standard for comprehensive analysis of human miRNAs. In this study, deep sequencing identified a large set of miRNAs that are differentially expressed between BC and normal epithelia.
Determining the expression signatures of miRNAs can rapidly and reliably reveal miRNAs that are aberrantly expressed in  cancers. Therefore, we previously analyzed miRNA expression signatures through the use of PCR-based arrays and searched for tumor suppressive miRNAs in BC. In this way, we successfully identified tumor suppressive miRNAs such as miR-1, miR-133a, miR-145 and miR-218 [18][19][20]25,[32][33][34][35][36]. We compared our previous results with the deep sequencing signatures described herein, allowing us to include the original data set with the new signature. To date, there have been four studies of BC miRNA expression signatures by deep sequencing analyses [9][10][11][12]. In published data, miR-1, miR-145, miR-143, miR-125b, and miR-100 have been listed as the most frequently downregulated miRNAs [9][10][11][12], data confirmed by our own signature. This fact supports the effectiveness of the deep sequencing signature, and it suggests that there are novel tumor suppressive miRNAs in these new data. The advantage of this strategy is to be able to identify new miRNA candidates in BC that have not been previously reported. We found 17 such sequences that might constitute new miRNA candidates. However, they were detected at very low levels of expression in comparison with known miRNAs in both cancer and benign tissues. The functional significance of these new miRNA candidates is unknown, and it will be important in the future to examine the relationship between these candidate miRNAs and BC oncogenesis. In the current study, we did not use tissue microdissection, therefore the precise proportion of epithelial cells in the normal samples and the precise proportion of tumor cells in the tumor samples are unknown. However, the thin-flap of normal bladder mucosa was consisted mostly of NBE cells, whereas tumor samples were constituted primarily of cancerous epithelial cells. Therefore, we believe that minor contamination of other cells was not significant and did not affect expression levels of the microRNAs.
MicroRNAs can be categorized into families based on their sequences or into clusters based upon their genomic loci. Clustered miRNAs are transcribed coordinately and have similar functions regulating the same targets [37]. For example, expression of the miR-1/133a cluster frequently reduced several types of cancers including BC. This cluster functions as a tumor suppressor targeting the same oncogenic genes such as TAGLN2 and PNP [18][19][20][21][22][23]. Similar examples have been reported with other clustered miRNAs such as miR-15a/16 and miR-143/145 [14][15][16][25][26][27][28]. Thus, these studies showed that clustered miRNAs might work in combination to accomplish their function throughout the same biological processes and the same targeted genes. With regard to the miR-195/497 cluster, it has been shown that miR-195 was downregulated in several cancers [38][39][40] and inhibited glucose uptake, proliferation, and cell cycle by targeting GLUT3 and CDK4 in BC [41,42]. It has also been demonstrated that miR-497 was downregulated and functioned as a tumor suppressor by targeting BCL2 in several cancers [43,44]. The miR-195/497 cluster also functioned as a tumor suppressor by targeting RAF1/CCND1 in breast cancer [45]. To validate these past reports, we evaluated the expression levels of miR-195 and miR-497 in BC clinical specimens   miRNAs are unique in their ability to regulate many protein coding genes. Bioinformatic predictions indicate that miRNAs regulate more than 30% of protein coding genes [29]. In cancer cells, normal miRNA -mRNA networks are disrupted by aberrantly expressed miRNAs. We have taken the position that identification of novel cancer pathways and responsible target genes regulated by the tumor suppressive miR-195/497 cluster is an important first step in understanding BC oncogenesis. Based on this view, we identified molecular pathways and target oncogenes that were regulated by the miR-195/497 cluster using genomewide gene expression data and in silico analysis. The miR-195/497 cluster belongs to the miR-15/16/195/424/497 family that shares the same 39-UTR binding seed sequence (AGCAGCA). The TargetScan database identified 6,730 genes that were candidate targets of the miR-195/497 cluster. KEGG pathway analysis indicated that those genes were implicated in 113 pathways, including ''Pathways in cancer'', ''MAPK signaling'', ''Endocytosis'', ''Insulin signaling'' and ''Wnt signaling''.
Among these pathways, we focused on ''Pathways in Cancer'' because of their direct relevance. According to our hypothesis, target genes of miR-195/497 should be upregulated in cancer cells. When we examined the expression levels of 104 genes in the ''Pathways in cancer'' group, 27 genes were upregulated in BC clinical specimens, suggesting that these genes were candidate miR-195/497 targets and functioned as oncogenes in BC. To confirm the reliability of our in silico analysis of potential target genes, we checked whether BIRC5 and WNT7A genes were in fact regulated by the miR-195/497 cluster in BC. Our data demonstrated that BIRC5 and WNT7A were directly regulated by the miR-195/497 cluster in BC. BIRC5 is a member of the inhibitor of apoptosis (IAP) family preferentially expressed by many cancers including BC [46]. It can be measured in urine at mRNA and protein levels, and urine BIRC5 was reported as a diagnostic biomarker for BC [47]. Recent studies from other groups have shown that several microRNAs such as miR-542-3p, miR-218, miR-708 and miR-203 directly inhibited BIRC5 expression in different types of cancer [48]. Our study is the first to suggest that BIRC5 might be regulated by the miR-195/497 cluster in BC cells. Therefore, our strategy of identifying novel tumor-suppressive, miRNA-regulated molecular targets/pathways in BC is an effective approach for miRNA-cancer research.

Conclusions
In summary, we constructed miRNA expression signatures of clinical BC specimens using deep sequencing as a gold standard method. A total of 60 miRNAs were significantly reduced in BC specimens. We identified several downregulated miRNAs that were located within the same genomic region constituting a miRNA cluster, including miR-1/133a, miR-143/145, miR-99a/let-7c, miR-195/497 and miR-144/451a. Downregulation of the miR-195/497 cluster was a frequent event in BC clinical specimens. Restoration of these miRNAs significantly inhibited cancer cell proliferation, migration and invasion, suggesting that miR-195/ 497 function as tumor suppressors in BC. Our analysis of the miR-195/497 cluster suggested that it regulated putative cancer pathways and oncogenic genes. These findings will contribute to the analysis of BC oncogenesis. Identification of the tumorsuppressive miR-195/497 cluster in human BC could provide new information on potential therapeutic targets in the treatment of BC.

Clinical specimens and cell culture
Tissue specimens for miRNA screening using deep sequencing were obtained from five BC patients and five NBEs. BC patients had undergone cystectomy or transurethral resection of bladder tumors (TUR-BT) at Kagoshima University Hospital between 2003 and 2010. Five NBEs were derived from patients with non-BC diseases. No contamination of the cancer cells was detected by pathologists. Another series of 29 BCs and 20 NBEs were obtained from Kagoshima University Hospital between 2003 and 2010. This group was subjected to real-time RT-PCR to evaluate miRNA expression levels. The samples were staged in accordance with the tumor-node-metastasis classification system of the American Joint Committee on Cancer/Union Internationale Contre le Cancer (UICC) and histologically graded [49]. Written informed consent was obtained from all patients and our study was approved by the Bioethics Committee of Kagoshima University. The patients' backgrounds and clinicopathological characteristics are summarized in Tables 1 and 6.
We used two human BC cell lines: BOY, which was established in our laboratory from an Asian male patient, 66-years old, who was diagnosed with stage III BC with lung metastasis; and, T24, which was invasive and obtained from the American Type Culture Collection. These cell lines were maintained in minimum essential medium (MEM) supplemented with 10% fetal bovine serum in a humidified atmosphere of 5% CO 2 and 95% air at 37uC.

Tissue collection and RNA extraction
Tissue samples were immersed in RNAlater (QIAGEN, Valencia, CA, USA) and stored at 220uC until RNA extraction was conducted. Total RNA, including miRNA, was extracted from frozen fresh tissues using the mirVana TM miRNA isolation kit (Ambion, Austin, TX, USA) following the manufacturer's protocol. The integrity of the RNA was checked with the RNA 6000 Nano Assay Kit and a 2100 Bioanalyzer TM (Agilent Technologies, Santa Clara, CA, USA).

Small RNA library preparation and profiling by deep sequencing
Total RNAs from five BCs and five NBEs were subjected to deep sequencing. The experimental process included the following steps: small RNA isolation, cDNA library preparation, and sequencing. Total RNA of each sample was size-fractionated on a 15% PAGE gel, and an 18-30 nucleotide fraction was collected. The 59-RNA adapter was ligated to the RNA with T4 RNA ligase. Ligated RNA was size-fractionated, and the 36-50 nucleotide fraction was excised. The 39-RNA adapter was subsequently ligated to precipitated RNA using T4 RNA ligase. Ligated RNA was size-fractionated, and the 62-75 nucleotide fraction (small RNA+adaptors) was excised. Small RNAs ligated with adaptors were subjected to RT-PCR to produce sequencing libraries. PCR products were purified and cDNA libraries, 106-118 bp in length, were collected. Then, cDNA libraries were sequenced by Genome Analyzer IIx (GAIIx) (Illumina Inc., CA, USA) at Beijing Genomics Institute at Shenzhen.

Construction of the miRNA expression signature in BC
The 50 nucleotide sequence tags were produced using HiSeq sequencing. Then, the low quality reads were filtered out according to base quality value, the adaptor sequence was trimmed at the 39-primer terminus, and the 59-adaptor contaminants formed by ligation were removed. The high-quality clean reads that were larger than 18 nucleotides were mapped to the human genome using the SOAP program [50]. Small RNA tags were aligned to the miRNA precursor/mature miRNA of corresponding species in miRBase release 20.0. If there were no miRNA information for a particular species in miRBase release 20.0, small RNA tags were aligned to the miRNA precursor/ mature miRNA of all plants/animals in miRBase release 20.0. To identify small RNA tags originating from coding exons, repeats, rRNA, tRNA, snRNA, and snoRNA, NCBI GenBank data and Rfam data were used. To identify new miRNA candidate genes, all hairpin-like RNA structures encompassing small RNA tags were identified using Mireap (http://sourceforge.net/projects/ mireap/). Assessment of differentially expressed miRNAs was determined with the edgeR program, version 3.0.8 [51] with the general linear model (GML) method. P-values are adjusted for multiple testing using the Benjamini and Hochberg method [52] and P-values with false discovery rate (FDR),0.01 are considered as significant. Predicting target genes and cancer pathways of miR-195 and miR-497 We applied in silico analysis to investigate the biological significance of differentially expressed miRNAs obtained by deep sequencing. The workflow for in silico analysis of target genes is shown in Figure 4. First, we identified predicted target genes of miR-195 and miR-497 using the TargetScan database [http:// www.targetscan.org/]. We performed GeneCodis3 [http:// genecodis.cnb.csic.es/] analysis using all candidate genes to identify the biological processes or pathways potentially regulated by miR-195 and miR-497 [53]. The GeneCodis3 software assigned a number of the putative target genes to known pathways in KEGG [http://www.genome.jp/kegg/pathway.html]. These data facilitated understanding of miRNA-regulated molecular networks in human cells. We performed gene expression analyses of all candidate genes involved in each of the pathways identified by GeneCodis3 software analysis using microarray expression data, which were approved by the Gene Expression Omnibus (GEO) and were assigned GEO accession numbers (GSE11783 and GSE31684). In the microarray expression data, we examined 90 BCs and six NBEs collected from patients, none of whom had been exposed to chemotherapy before surgery. Statistical analyses were conducted using the Mann Whitney U-test.

Mature miRNA transfection
As described previously [32], the BC cell lines were transfected with Lipofectamine RNAiMAX transfection reagent (Invitrogen) containing Opti-MEM (Invitrogen) with ten nM Pre-miR TM or mirVana TM miRNA Mimic Negative Control (Applied Biosystems) that are mature miRNA molecules for gain-of function experiments. Cell were seeded in six-well plates for wound healing assays (20610 4 per well), in 24-well plates for Matrigel invasion assays (5610 4 per well), and in 96-well plates for XTT assays (3,000 per well).

Cell proliferation, Matrigel invasion, and wound healing assays
Cell proliferation was determined by using an XTT assay (Roche Applied Sciences, Tokyo, Japan) performed according to the manufacturer's instructions. A cell invasion assay was carried out using modified Boyden chambers consisting of Transwell-precoated Matrigel membrane filter inserts with eight-mm pores in 24-well tissue culture plates (BD Biosciences, Bedford, MA, USA). MEM containing 10% fetal bovine serum in the lower chamber served as the chemoattractant, as described previously [32]. Cell migration activity was evaluated by wound healing assay. Cells were plated in six-well dishes, and the cell monolayer was scraped using a P-20 micropipette tip. The initial gap length (zero h) and the residual gap length 24 h after wounding were calculated from photomicrographs. All experiments were performed in triplicate.
Plasmid construction and dual-luciferase reporter assay miRNA target sequences were inserted between the XhoI-PmeI restriction sites in the 39-UTR of the hRluc gene in the psiCHECK-2 vector (C8021; Promega, Madison, WI, USA). BOY cells were transfected with 15 ng vector, ten nM miRNAs (Pre-miR TM or mirVana TM miRNA Mimic Negative Control), and one mL Lipofectamine 2000 (Invitrogen) in 100 mL Opti-MEM (Invitrogen). The activities of firefly and Renilla luciferases in cell lysates were determined with a dual-luciferase assay system (E1910; Promega). Normalized data were calculated as the ratio of Renilla/firefly luciferase activities. Detailed vector sequences are shown in Table S5.

Statistical analysis
The relationships between two variables and numerical values were analyzed using the Mann-Whitney U test, and the relationship between three variables and numerical values was analyzed using the Bonferroni-adjusted Mann-Whitney U test. The relationship between miR-195 and miR-497 expression was analyzed by the Spearman rank correlation. Expert Stat View analysis software (ver. 4; SAS institute Inc., Cary NC) was used in these analyses. In the comparison of three variables, a nonadjusted statistical level of significance of P,0.05 corresponded to the Bonferroni-adjusted level of P,0.0167.

Supporting Information
File S1 Tables S1-S5. Table S1. The number of total reads and expressed reads of 'known miRNAs' in deep sequencing.