MicroRNA Expression Profile in Penile Cancer Revealed by Next-Generation Small RNA Sequencing

Penile cancer (PeCa) is a relatively rare tumor entity but possesses higher morbidity and mortality rates especially in developing countries. To date, the concrete pathogenic signaling pathways and core machineries involved in tumorigenesis and progression of PeCa remain to be elucidated. Several studies suggested miRNAs, which modulate gene expression at posttranscriptional level, were frequently mis-regulated and aberrantly expressed in human cancers. However, the miRNA profile in human PeCa has not been reported before. In this present study, the miRNA profile was obtained from 10 fresh penile cancerous tissues and matched adjacent non-cancerous tissues via next-generation sequencing. As a result, a total of 751 and 806 annotated miRNAs were identified in normal and cancerous penile tissues, respectively. Among which, 56 miRNAs with significantly different expression levels between paired tissues were identified. Subsequently, several annotated miRNAs were selected randomly and validated using quantitative real-time PCR. Compared with the previous publications regarding to the altered miRNAs expression in various cancers and especially genitourinary (prostate, bladder, kidney, testis) cancers, the most majority of deregulated miRNAs showed the similar expression pattern in penile cancer. Moreover, the bioinformatics analyses suggested that the putative target genes of differentially expressed miRNAs between cancerous and matched normal penile tissues were tightly associated with cell junction, proliferation, growth as well as genomic instability and so on, by modulating Wnt, MAPK, p53, PI3K-Akt, Notch and TGF-β signaling pathways, which were all well-established to participate in cancer initiation and progression. Our work presents a global view of the differentially expressed miRNAs and potentially regulatory networks of their target genes for clarifying the pathogenic transformation of normal penis to PeCa, which research resource also provides new insights into future investigations aimed to explore the in-depth mechanisms of miRNAs and other small RNAs including piRNAs in penile carcinogenesis regulation and effective target-specific theragnosis.


Introduction
Molecular research has emerged as a promising tool for in-depth understanding into carcinogenesis, with its cellular components, signaling pathways and pharmaceutical targets in cancer therapeutics. Unfortunately, the remarkable progress has not encompassed all cancer types. Penile cancer (PeCa) is such a relatively rare tumor entity within the developed countries of the world, only about 1,640 penile cancer cases were newly diagnosed and 320 cancer-related deaths were estimated in 2014 in the United States nationally [1]. Whereas, the morbidity and mortality rates of PeCa in most developing countries were considerably higher, accounting for a 3-10 fold increase at the rates of developed countries over the past decades [2]. Penile cancer is associated with several established risk factors such as human papillomavirus (HPV) infection, poor hygiene, phimosis with chronic inflammation, smoking and some epigenetic alternations including histone methylation modifications [3,4]. Besides, a few genes during carcinogenesis, proliferation, invasion and metastasis of PeCa have been identified. For instance, p53, cyclin-dependent kinase inhibitor 2A (CDKN2A) and matrix metallopeptidase 9 (MMP-9) were proven to play crucial roles in PeCa carcinogenic routes and epithelial-mesenchymal transition (EMT), respectively [5]. Nonetheless, as mentioned above, additional comprehensive research to further elucidate the exact mechanisms of molecular changes underlying the initiation, development and progression of PeCa is imperative for us to explore novel and valuable preventive, early detection, targeted therapies and prognosis prediction approaches for this genitourinary disorder.
Mature microRNAs (miRNAs) are a group of endogenously expressed, evolutionally conserved, single-stranded and non-coding small RNAs (smRNAs, approximately 19-23 nucleotides in length) that act as posttranscriptional regulators of gene expression in a wide range of organisms [6]. To date, more than 1000 human miRNAs have been identified and tremendous observations have been achieved in connecting the alternative and aberrant expression levels of miRNAs with the initiation and progression of human diseases, especially for various kinds of cancer types [7][8][9], as no less than half miRNA genes are located in cancer-associated regions or fragile sites in the whole genome [10]. miRNA deregulation in tumor biology was first observed in miRNA-15a and miRNA-16-1 within locus 13q14, both of the two miRNAs which targeted B-cell lymphoma 2 (Bcl-2) mRNA were deleted or downregulated in most majority of the chronic lymphocytic leukemia (CLL) cases, resulting in tumorigenesis by reducing apoptotic activities [11]. Presently, it has been widely accepted that distinct cancer types, stages, regression grades or differentiation states might possess individual miRNA expression profiles, which hold great promise in acting as reliable molecular biomarkers for cancer diagnosis and may also reveal new pathogenic signaling pathways correlated with carcinogenesis and progression for targeted molecular therapeutics [12]. Ever since, a variety of techniques have been introduced to conduct miRNA profiling. For instance, quantitative reverse-transcription polymerase chain reaction (qRT-PCR) assays, Northern blotting analyses and microarrays have ever been extensively used for miRNA profiling studies [13]. Recently, a newly developed technology named next-generation sequencing (NGS) has attracted much attention in small RNA (smRNA) profiling due to its unique advantages in test specificity, sensitivity and novel smRNAs identification [14][15][16]. Through which, we can detect almost all smRNAs such as known and novel miRNAs even with extremely low abundance, small cytoplasmic RNAs (scRNAs), nuclear RNAs (snRNAs), nucleolar RNAs (snoRNAs) and piwi-interacting RNAs (piRNAs) [17].
In the present study, we aimed at filling the gap in information regarding the miRNA profile in human penile cancer and evaluated a possible correlation between differentiated miRNA expression levels with tumorigenesis and progression. By applying next-generation sequencing (NGS) technology, we performed smRNA-seq for ten paired fresh-frozen specimens of penile squamous cell carcinoma and matched histologically normal penile tissues which were adjacent to the cancer. We have addressed the general information of the miRNA expression profiles in both the paired penile tissues and especially found that a number of miRNAs were aberrantly expressed in penile cancers compared with matched normal tissues. Notably, the gene ontology (GO) analysis of potential miRNA target genes indicated that deregulated miRNAs in enriched biological processes, molecular functions and cellular components were involved in cell growth, cell shape, axonogenesis, protein activity regulation and angiogenesis, which together participated in the transformation of normal cells to malignant lesions. Besides, the kyoto encyclopedia of genes and genomes (KEGG) pathway analysis successfully enriched several tightly cancer-associated pathways which were all well-documented to participate in cancer initiation, development, invasion, metastasis and also promising therapy. Our results provide a valuable research resource to further elucidate the regulatory roles of miRNAs in penile tumorigenesis and progression, which also holds great promise in facilitating the development of novel diagnostic biomarkers and effective therapeutic strategies for PeCa.

Ethics statement
Informed consents were signed by all the thirteen patients themselves with penile squamous cell carcinoma who had received partial penectomy and the investigation was evaluated thoroughly and then approved strictly according to the regulations by the Ethics Committees on Human Research of the First Affiliated Hospital of Anhui Medical University (Approval No. PJ20140906).

Clinical specimens collection
Penile cancer and matched non-malignant specimens were collected from a total of 13 patients who underwent partial penectomy from 2013 to 2015 at Department of Urology, the First Affiliated Hospital of Anhui Medical University. Immediately after surgery, cancerous and matched normal penile tissues were separately stored in RNAlater solution (Ambion, USA) and frozen at -80°C according to a standard operation procedure (SOP) guided by the uropathologist. Briefly, tumor tissues containing more than 80% tumorous lesions were included for further analyses. Correspondingly, the matched normal adjacent tissues were obtained from where at least 1.0 cm apart from the visible tumorous tissues in penis. Each cancerous and non-malignant tissues were diagnosed histopathologically by Hematoxylin-eosin (HE)-staining.

SmRNA library construction and deep sequencing
Total RNA was collected from ten pairs of fresh-frozen penile cancers and matched histologically normal penile tissues by using TRIzol reagent (Life Technologies, USA). Qualified RNA with A260/A280 ratio more than 1.8 was subsequently assessed for RNA integrity by conducting the DNA chip assay. Two smRNA libraries were constructed from pooled and integrated RNAs of ten individuals for each group (cancerous and matched normal penile tissues) by using a smRNA sample prep kit (Illumina, USA) following the standard protocols with specific modifications [15]. Briefly, the appropriate fractions ranged from 18 nt to 28 nt were separated, purified via 15% (w/v) denaturing polyacrylamide gel electrophoresis and then linked to RNA adaptors (5'-GTTCAGAGTTCTACAGTCCGACGATC and 3'-TCGTATGCCGTCTTCT GCTTG) followed by RT-PCR amplification. Next, the PCR products were further purified on agarose gels to establish the libraries. Finally, the two libraries which constructed from ten cancerous penile tissues (pooled and homogenized to a library) and ten matched adjacent normal penile tissues (pooled and homogenized to another library) were sequenced by applying Illumina Hiseq 2000 (Illumina, USA) at Beijing Genomics Institute (Shenzhen, China) strictly following the standard protocols [18].

NGS data analyses and novel miRNAs identification
The smRNA sequencing data were further analyzed referring to the methods described in the previous publications [15][16][17][18]. Concretely, after deleting 5' adaptor and trimming 3' acceptor sequences, removing contaminated and low-quality reads (Q<10, the quality value was calculated by Q = ASCII character code-64) such as reads without insert fragment or containing poly(A) stretches, the qualified unique tags were mapped to GRCh37.p5 by introducing SOAP2.0 ultrafast tool (http://soap.genomics.org.cn) [19]. Tags with no more than one mismatch were picked out for subsequent analyses. The conserved known miRNAs were identified by aligning the mapped tags to miRBase tool (version 18.0) in Homo sapiens [20], and the tags not matched in miRBase were further aligned against the sequences of other non-coding RNAs (rRNAs, tRNAs, snRNAs, snoRNAs) presented on Rfam database (http://www.sanger.ac.uk/ resources/databases/rfam.html) [21], repeats database as well as piRNA database [22]. In consideration of a few smRNA tags might map to more than one category, we followed the priority rules described in the previous publication to guarantee that every unique smRNA was mapped to unique annotation as follows: miRNA, Rfam, repeats, piRNA and mRNA (122.228.158.106/ mr2_dev) [23].
After identifying the known miRNAs, the remaining sequences from all cancerous and matched normal penile specimens not annotated by the above-mentioned databases were picked out of the two libraries to be aligned with the integrated human transcriptome for predicting novel miRNAs. All hairpin-like structures containing unclassified smRNA tags consisting of no less than 45 reads were predicted using classical pre-microRNAs prediction tool Mireap (http://sourceforge.net/projects/mireap/) [24] strictly according to the following rules: 1) Length of miRNA sequences range from 18 to 26 nt; 2) Free energy of a miRNA precursor is less than -18 kcal/mol; 3) Bulge allowed for miRNA and the paired precursor miRNA Ã is less than 4; 4) Space between miRNA and the paired precursor miRNA Ã is no more than 35 nt. Meanwhile, RNA Fold (http://rna.tbi.univie.ac.at/RNA/RNAfold.html) [25], a program computes the minimum free energy (MFE) and backtraces an optimal secondary structure, was lent to predict all the essential secondary structures of potential miRNA precursors. All the raw data generated by next generation sequencing has been deposited in the designated database ArrayExpress (Accession number: E-MTAB-3087).

Bioinformatics analyses for differentially expressed miRNAs
The approaches described in our previous publication for bioinformatics analyses were introduced to perform the target genes prediction in the present study [17][18]. Briefly, the targeted genes of differentially expressed miRNAs from cancerous and matched normal penile tissues were predicted using MicroCosm Targets, formerly called miRBase targets (http://www.ebi.ac. uk/enright-srv/microcosm/) [26][27], TargetSpy (http://targetspy.org/) [28] and miRanda (http://www.microrna.org) [29]. The concrete predictions were strictly according to the rules described as follows: 1) Any mismatch should not be found at the seed region; 2) The free energy of the miRNA/target duplex should be less than -20 Kcal/mol; 3) The total score for a miRNA-mRNA pair should exceed 140 [17][18].
Subsequently, the enriched GO terms (http://www.geneontology.org) [30] and KEGG pathways (http://www.genome.jp/kegg/) [31] analyses were conducted for the predicted target genes of differentially expressed miRNAs from normal and cancerous penile specimens. Concretely, after mapping the potential target genes to the dataset of GO annotations and KEGG pathways, sorting out the enriched biological processes and signaling pathways through Hypergeometric and Fisher's tests, the enrichment ratios and p-values for GO and KEGG were calculated via introducing the identical formula mentioned in our previous report [17], then the False Discovery Rate (FDR) adjustment was performed to judge the significance of differences in multiple testing [32]. Eventually, the key GO terms and pathways were identified only under the circumstance that enrichment ratio!2.0 and FDR p-value<0.05 simultaneously.

qRT-PCR verification for altered miRNA expression
Mature miRNAs quantification was performed by routine quantitative real-time PCR using an Applied Biosystems StepOne Real-Time PCR System (Applied Biosystems, USA) and a SYBR premix Ex-Taq II kit (Takara, Japan) with optimized reaction conditions and using the specific primers listed in S1 Table. Briefly, total RNA was extracted from 13 paired cancerous and matched normal penile tissues. Among which, the extracted RNA from the identical ten patients applied to NGS was pooled and validated, while each total RNA from another three patients was separately used for further analyses. After the cDNA synthesis, triplicate reactions were performed at 95°C/10 min, the subsequent 40 amplified cycles were conducted at 95°C/15 s and 60°C/60 s. Meanwhile, U6 snRNA was lent as an internal control to normalize the miRNA expression. For data analysis, the threshold cycle (Ct) was defined as the fractional cycles when the fluorescence passed a fixed threshold and further applied to calculate relative miRNA abundances (44Ct = Ct miRNA -Ctsn RNAU6 ). Relative expression fold changes were calculated using 2 -44 Ct formula [16]. The miRNA concentration differences between cancerous and normal penile tissues were analyzed using unpaired t-tests [18]. When a p-value was found to be less than 0.05, it would be considered statistically significant.

Results
Overview of whole genome smRNA-sequencing data from the paired penile specimens All smRNAs of 18-32 nucleotides from the paired fresh-frozen specimens of penile squamous cell carcinoma (hereafter referred to as "cancerous tissues") and matched histologically normal penile tissues adjacent to the cancer (hereafter referred to as "normal tissues"), which have been diagnosed histopathologically by HE-staining (Fig 1), were deep-sequenced in order to obtain a comprehensive view of smRNAs profile.
For each sample, 13,796,889 (out of 15,702,009 reads) and 14,051,355 sequence reads (out of 15,698,197 reads) aligned to the human genome sequence dataset were obtained, representing 704,514 (out of 966,914) and 787,447 (out of 1,090,353) unique tags, respectively (S2 Table). Among which, 6,898 unique tags corresponding to 4,795,955 reads and 7,173 unique tags corresponding to 3,815,482 reads were matched to known miRNAs for normal and cancerous tissues, respectively. Meanwhile, 5,877 corresponding to 172,216 reads and 6,272 corresponding to 119,724 reads were matched to known piRNAs for normal and cancerous tissues, respectively. The rest of the sequences were matched into other types of RNAs, including mRNAs, rRNAs, sRNAs, snRNAs, snoRNAs, tRNAs, repeats and so on (Fig 2 and S2  Table). Among which, the repeats mainly consisted of rRNA either based on total reads or unique tags (S3 Table).
The majority of these clean reads were 22 nt in size, consistently followed by 23-nt and 21-nt RNA fragments (S1 Fig). Furthermore, to understand the distribution pattern of all the unique tags on chromosomes, the detailed location of each tag on chromosome was evaluated. As a result, the chromosome distributions between the cancerous and normal tissues were
To identify potential novel miRNAs in penile tissues, the unclassified tags were further processed using Mireap tool. Only those tags with reads count more than 50 and strictly matching the default parameters were classified as candidates of novel mature miRNAs. On the basis of this analysis, each 10 miRNA precursors (whose potential stem-loop structures were shown in S3 Fig) encoding the most abundant predicted novel mature miRNA candidates with lengths ranged from 22 to 24 nt were separately identified in the two smRNA libraries from penile tissues. Notably, 7 abundantly expressed novel miRNAs were overlapped between cancerous and matched adjacent normal penile tissues (Table 2), which indicated the identically histological origin and chromosomal location of the miRNAs in both tissues once again. Meanwhile, these novel miRNAs originated from same precursors have also shown heterogeneity and preference, by generating relatively more mature products from the 3' ends compared with the 5' ends and among the top 10 novel mature miRNAs in either penile tissues, 9 miRNAs initiated with adenine or uridine (Table 2), which patterns were consistent with the previous report [17].

Differentially expressed profiles of miRNAs in penile cancers
When comparing the distinct levels of miRNA abundance between normal and cancerous tissues, two quantities are frequent of great importance: one is the p-value corresponding to the ttest responses to the issue whether the statistically significant difference of expression levels between the paired groups can be achieved. The other one is the fold-change (FC) value, which could illustrate what is the magnitude of the miRNA abundance difference between the paired specimens. Notably, even very large FC values may have insignificant p-values according to the t-test caused by considerable variabilities. Similarly, relatively small magnitude of the difference may also mean a tiny p-value due to the highly consistent technical replicates within each samples. Thus, a useful visualization method called "volcano plot", which can simultaneously exhibit fold-change and p-values on the X-axis and Y-axis separately, has been introduced to conduct our analysis. Generally, points located in the upper-left or upper-right regions of the plot and indicated with the green color should draw more attention, as they had both small p-values (p<0.05) and high FC values (FC!2). In this comparative study, the overview of the volcano plot generated by miRNAs profiles in penile cancer and matched normal tissues were shown (S4 Fig). As a result, significantly differential expression levels of 98 miRNAs were found between PeCa and adjacent normal penile specimens after adjustment for multiple testing. By applying additional filtering with a stringent cut-off of 50 reads count in both samples, we identified 56 miRNAs consisted of 30 (53.6%) miRNAs that were downregulated (Table 3) and 26 (46.4%) that were upregulated (Table 4) in the patients and could discriminate penile carcinoma from matched normal tissues.

Altered miRNA expression in cancerous penile tissues validated by qRT-PCR
To validate the altered miRNA expression in PeCa profiled by smRNA-seq, qRT-PCR was performed on the penile cancer tissues and matched histologically normal tissues. Firstly, to exclude the possible individual differences between the samples subjected to the sequencing, the 10 paired penile specimens from the patients were also mixed as a pool for the validated assay. Subsequently, we randomly picked 5 upregulated and 5 downregulated miRNAs among the differentially expressed miRNAs to conduct the qRT-PCR analysis and the expression levels were presented using the normalized fold change of cancerous tissues versus normal tissues. The results revealed that the expression levels and altered expression tendency of the 10 deregulated miRNAs were quite consistent with the results obtained from NGS technology. From both techniques, miR-509-3p showed the largest fold-change of downregulated expression levels in cancerous penile tissues, and miR-107 possessed the largest fold-change of upregulated expression levels in cancerous penile tissues, which together indicated that the expression levels of miRNAs detected by smRNA-seq sequencing were reliable (S5 Fig). Moreover, to validate the altered miRNA expression pattern revealed by the NGS data in individual patients, we picked 4 deregulated miRNAs to conduct the qRT-PCR analysis and the expression levels were presented using the normalized fold change of cancerous tissues versus normal tissues in each patient, respectively. According to the result, we found that although the population-based effects and inherent differences inevitably emerged as the previous report [33], the tendency of the expression patterns was consistent with the NGS findings, as miR-107 and miR-223-3p showed the upregulated expression levels in cancerous penile tissues, while miR-1247-5p and miR-509-3p showed the downregulated expression levels in cancerous penile tissues, which again verified the smRNA-seq sequencing data (S6 Fig). GO enrichment analysis of genes differentially expressed in cancerous and normal penile tissues The essential biological functions of the putative target genes were classified via the Gene Ontology system. Since individual genes were related to distinct GO terms, gene sets which showed the similar deregulated expression patterns in multiple cancers were thought to be functionally crucial to the tumorigenic processes [30]. In order to better understand the biological behaviors of deregulated known miRNAs, potential target genes and biological functions targeted by these miRNAs were predicted by using MicroCosm, TargetSpy and miRanda softwares [26][27][28][29] (S1 File). After adjusting the number of total genes count annotated by a specific GO term!10, enriched ratio!2.0 and FDR p-value<0.05 to obtain the most convincing result for enrichment analysis, the top 10 most enriched GO terms including biological processes, cellular components and molecular functions were presented in Table 5. Among which, the putative target genes of deregulated miRNAs between the cancerous and normal penile specimens appeared to be mostly associated with the enriched cellular components consisted of cytoskeletons such as adherens junction (which junction anchors transmembrane proteins and facilitates cell cytoskeletons attachment), stress fiber (which consists of short actin filaments), cell-cell junction (which communicating junction commonly connects adjacent cells). Meanwhile, the most enriched molecular function among the top 10 GO terms was beta-catenin binding (which functions as the interactive behaviour with beta-catenin in a selective and noncovalent manner). Besides, these putative target genes were tightly related to a wide range of biological processes including growth regulation (which extensively modulates the organism's development), cell shape regulation (which functions as the cell surface configuration commander), axonogenesis (which process generally modulates the axon's morphogenesis and shape determination), cyclin-dependent protein kinase activity regulation (which procedure regulates serine/threonine kinase activity comprehensively), peptidyl-serine phosphorylation (which process mediates the conversion of peptidyl-serine to peptidyl-O-phospho-L-serine) and positive angiogenesis regulation (which procedure facilitates the formation of blood vessels).

KEGG pathway analysis of genes differentially expressed in cancerous and normal penile tissues
After GO analysis, the putative target genes of these miRNAs differentially expressed in cancerous and normal penile specimens were introduced to KEGG to conduct a pathway enrichment analysis, which eventually showed 22 target genes that were annotated for differentially expressed miRNAs, and 5 KEGG pathways were enriched. Based on the statistically significant definition with FDR p-values<0.05, the annotated most enriched pathways were found to be tightly linked to "dorso-ventral axis formation (GRK-EGFR)" (which modulates the key effectors in axis formation such as Rho and Notch) and "colorectal cancer" (which is involved in the carcinogenic signaling transductions of genomic alternations resulted from chromosomal and microsatellite instabilities) ( Table 6). All together, these enriched pathways participated in cell cycle, cell survival, proliferation, growth, differentiation, genomic instability and also apoptotic activity. Especially, we have summarized the specific signaling pathways participated in the comprehensive modulation of the putative target genes of the differentially expressed miRNAs between the cancerous and adjacent normal penile tissues, as follows: TGF-beta signaling pathway, Wnt signaling pathway, MAPK signaling pathway, p53 signaling pathway, PI3K-Akt signaling pathway and Notch signaling pathway (Detailed regulatory networks were available at KEGG pathway database on http://www.genome.jp/kegg/pathway.html [31]).

Discussion
Penile cancer is associated with a number of potential risk factors including HPV infection, phimosis with chronic inflammation, poor hygiene and smoking, which disease predominantly affects aged males between 50 and 70 years old in most developing countries by high morbidity and mortality [34][35]. Although some oncogenes, chromosomal aberrations, epigenetic alterations, regulatory pathways and several core machineries involved in the penile cancer carcinogenesis, progression and invasion have been identified over the past decades [36][37], fine- tuning of the definite molecular concepts of penile cancer and especially the exact mechanisms involved in the development and progression of this highly mutilating disease remain largely pending. In-depth researches to further elucidate the detailed molecular events underlying the initiation and progression of penile cancer and its regulatory signaling pathways are eagerly necessary to direct prevention, early detection, targeted therapy and prognosis prediction.
Notably, an increasing number of studies indicated that tumorigenesis, progression and invasion were regulated transcriptionally and post-transcriptionally in a strict and also cooperative manner by small non-protein-coding RNAs [38][39]. Among the variety kinds of smRNAs, microRNA (miRNA) molecules are a class of evolutionally conserved and endogenously expressed RNAs, which have been emerging as critical regulators that function in gene expression at posttranscriptional level [40][41]. Accumulating evidence has indicated that particular miRNAs either act as tumor oncogenes or suppressors, whose loss or overexpression aberrantly, could not only differentiate cancers arising from distinct physiological sites but also possess considerable prognostic significances [42][43]. Accordingly, as malignant tumors have shown additional traits beyond the acquisition of enhanced growth potency, decreased cell death, sustained angiogenesis as well as the invasion and metastasis ability, it is convincing that miRNAs could also act as master regulators of cancer therapy [44]. Thus, to generate the differentially expressed profiles of smRNAs in human penile cancer tissues and their matched adjacent normal tissues is prerequisite for thorough understanding of the definite regulatory roles of non-coding RNAs in penile cancer initiation, progression and invasion processes.
Recently, the newly developed next-generation sequencing (NGS) technology has been widely used for miRNA quantitative studies, by overriding the limitations in test specificity and inability to detect novel or low-abundance miRNAs such as rno-miR-509-5p/3p and rno-miR-1306-5p which can not be detected by the traditional quantitative assay [45]. Based on NGS, bioinformatic analyses of distinct miRNAs expression signatures and putative miRNAs binding sites have indicated several novel potential gene targets and signaling pathways of deregulated miRNAs in various cancers such as colorectal cancer [46] and breast cancer [47]. Particularly, for the genitourinary cancers, the comprehensively comparative miRNA expression profiles of prostate cancer [14,48], bladder cancer [16], renal cancer [33], testicular germ cell tumors [49] and adrenocortical tumors [50] have been well-documented. Whereas, no any systematic studies have been conducted on miRNAs profiling in human penile cancer by using deep sequencing or even microarray.
In this present study, via utilizing state-of-the-art next-generation sequencing technology on smRNA libraries prepared from ten human penile cancer tissues and their matched adjacent tissues confirmed by the histopathologic diagnosis and then sequenced by employing the same protocols on the identical platform to guarantee the negligible inter-laboratory or interplatform technical variations [18,49], we have quantified 14,051,355 and 13,796,889 total sequencing reads aligned to the human genome sequence dataset, which represented 787,447 and 704,514 unique tags for the cancerous and adjacent normal tissues. Among these unique tags, 806 and 751 genome-widely expressed known miRNAs accompanied with 211 and 189 novel miRNAs were identified in PeCa and normal penis, respectively. Therein, the majority of highly expressed miRNAs were let-7-5p miRNA family members, which expression pattern was extremely similar as the miRNA profile in human testis [17,51]. The possible explanation of the consistence was that mature form of let-7 family were highly conserved as their critical roles in developmental regulation of human reproductive organs including penis and testis were decided [52].
Furthermore, based on the "volcano plot" method which defined small p-values (p<0.05) and high fold change values (FC!2) simultaneously and FDR adjustment, we have promoted this comparative study by focusing on exploring the miRNAs which showed the significantly aberrant expression pattern in penile cancer compared with matched histologically normal tissues, as differentially-expressed miRNAs were widely considered to be attractively novel prognostic biomarkers candidates and available targets for the early theragnosis of penile cancer. As a result, quantitatively differential expression analysis allowed the identification of a 56-miRNA signatures distinguishing penile cancer from normal penis. To make these signatures more concrete, the most commonly-used assay qRT-PCR was applied to validate the results obtained from NGS [14][15][16][17][18]53]. Consequently, both the randomly-selected downregulated and upregulated miRNAs showed the almost unanimous altered expression pattern with the NGS data, indicating the miRNA profiles detected by the high-throughput sequencing technology were reliable.
Among the deregulated miRNAs, the most abundantly downregulated and upregulated miRNAs were miR-320a and miR-107, respectively. In consideration of the transformation of normal cells into malignant cancer cells may be governed by several common rules and undergo the identical biological processes [54], we have compared the deregulated miRNA expression patterns with other cancers with great interest. To begin with the most downregulated miRNA in penile cancer miRNA-320a, we found it often played as a tumor-suppressive miRNA in bladder cancer [55] and colorectal cancer [56] by targeting different oncogenes. Referring to the most abundantly and also significantly upregulated miRNA in penile cancer miR-107, it was consistently revealed the tumorigenic and metastatic potency to human colorectal cancer [57], pancreatic cancer [58] and gastric cancer [59]. Meanwhile, miR-508-3p, the most significantly downregulated miRNA in penile cancer, has been revealed a lower expression level in renal cancer tissues, in the same study, another downregulated miRNA named miR-509-3p which can suppress the proliferation of renal cancer cells was also found to be downregulated in penile cancer, together indicating its tumor suppressive role in cancer development [60].
To lend more concrete proofs to identify the potential modulatory roles of the miRNAs aberrantly expressed in penile cancer, we compared their expression patterns from NGS with previous reports in additional kinds of cancers. As expected, most majority of the miRNAs, either upregulated in the cancer tissues, such as miR-223-3p [61], miR-421 [62], miR-424-5p [63], miR-1260b [64] etc, or downregulated in the penile cancer, such as miR-205-5p [65], miR-211-5p [65][66], miR-365-3p [67] and miR-1247 [68] etc, were entitled the similar roles in carcinogenesis of bladder cancer, retinoblastoma, breast cancer, nasopharyngeal cancer, pancreatic cancer and several other cancer types [61][62][63][64][65][66][67][68]. Moreover, as miRNAs are highly tissue specific and can act as both tumor promoters or suppressors [69], we have then illuminated the definite roles of some "Janus miRNAs" in penile carcinogenesis regulation. For instance, expression of miR-23b-3p in cancer is somewhat controversial because it was either upregulated and oncogenic in renal cancer [70] or downregulated and execute suppressive effect in prostate cancer [71] and bladder cancer [72]. Similarly, recent studies have shown that miR-199a-5p played opposite roles in cancer initiation and development of various cancer types [73]. Our NGS data have revealed that miR-23b-3p was downregulated and miR-199a-5p was upregulated in penile cancer tissues, indicating their tumor-suppressive and oncogenic effect, respectively. More importantly, we have determined the regulatory function of several fresh miRNAs, including miR-887-3p, miR-1277-3p, miR-3120-5p etc, whose roles in carcinogenesis or cancer progression were largely unknown before.
miRNA databases enabled us to obtain the various currently validated targets of the differentially expressed miRNAs identified by our NGS data. Predicted target analyses suggested participation of the deregulated miRNAs in a number of biological processes involved in cell growth, axonogenesis, cell shape, protein activity regulation and angiogenesis. Accordingly, these top enriched GO terms were closely related to transformation of normal cells to malignant lesions. For instance, CDKN2A and MMP-9 among the putative target genes were proven to play crucial roles in penile carcinogenesis [5], notably, these two genes were successfully enriched in the regulation of cell growth, cyclin-dependent protein kinase activity and angiogenesis, respectively. Besides, IL-1A and IL-1B, which have been identified as the key machineries in inflammation modulation to cause penile cancer [3][4], were also successfully enriched in the current GO terms. Meanwhile, a remarkable GO term named axonogenesis, which was not usually associated with carcinogenesis but always played vital roles in the generation and outgrowth of axons during neuron development, was surprisingly enriched. Intriguingly, axonogenesis has been a recently described phenomenon of paramount importance in prostate cancer by participating in the proliferation and hormonal regulation of prostate cancer [74]. Therefore, our exploratory result has listed the predicted target genes and hinted a similar role of axonogenesis involved in penile cancer progression. Simultaneously, we have addressed that the enriched cellular components and molecular functions mainly consisted of cytoskeletal structures such as cell junctions, stress fibers and relevant beta-catenin binding regulation, respectively. This was reasonable since cellular morphogenesis, which determined cell shape and polarity, cell-cell junctions as well as cell-matrix adhesions through their interactions with cadherins and integrins, has been well established to accompany progression and invasion after carcinogenesis [75].
According to the KEGG pathway analysis, the tightly cancer-associated pathway "colorectal cancer" was successfully enriched. From which, the containing TGF-β [76], Wnt [77], MAPK [78] PI3K-Akt [79], p53 [80] and Notch signaling pathways [81] were previously documented to be involved in cancer development, progression, invasion, metastasis and also promising therapeutic targets screening. Besides, another pathway was also enriched. Although it was not usually considered as the direct carcinogenesis-modulated networks, the core molecular machinery Rho employed in the pathway was also definitely related with cancer biology [82]. Thus, our results have provided additional in-depth insights into the system-level molecular mechanisms of the tumorigenesis of penile cancer.
In addition, we have identified a distinct subset of 23-28-nt length smRNAs named piRNAs, which also differentially expressed in human penile cancers versus their matched adjacent normal penile tissues. These piRNAs are suggested to bind PIWI proteins to mediate posttranslational control of transposable elements expression and modulatory a crucial role in epigenetic changes via interaction with crucial proteins such as MAEL [83]. piRNAs have previously identified as smRNAs expressed predominantly in the germline of mammals [17,84]. However, there were mounting evidences revealed piRNAs pathway activities in various cancers including breast cancer, gastrointestinal tract cancer, endometrioma and ovarian cancer recently [85][86]. For instance, piR-651 expression was upregulated in gastric cancerous tissues compared with paired non-cancerous tissues and the piR-651 level was further associated with cancer TNM stage, indicating piRNAs may also act as potential biomarkers for various cancer diagnosis [86]. Thus, in the present study, we have summarized the top 10 most abundant piR-NAs in cancerous and adjacent normal penile tissues (S4 Table). From which, we found 6 out of 10 top abundant piRNAs expressed in paired penile tissues, again indicated the tissues with identically histological origin possessed the similar expression profiles of the smRNAs (Table 1). Simultaneously, we have compared the differentially expressed piRNAs in the paired penile tissues via a stringent cut-off of 50 reads count in either sample under the visualization "volcano plot" method (S5 Table and S7 Fig). Future studies will help to elucidate the detailed biogenesis and regulatory functions of these piRNAs in penile cancer initiation and progression.
Another potential advantage of our sequencing effort was the capability to evaluate possible miRNA modification events. Several studies have described miRNAs were modified to a significant degree in adult mouse brain and human testis, which modifications may enhance the miRNA stability, reinforce the specific miRNA-mRNA interactions and also execute some regulatory behaviors [17,87]. In our current study, we also found some miRNAs from human penile tissues with similar modification profiles, among which, 6 modifications were overlapped in the top ten "5' end added" modifications in both of the paired specimens, indicating the same tissue origin not only possessed the similar expression profiles of the smRNAs, but also showed the highly similar modification pattern [18]. Intriguingly, some miRNAs such as miR-320a, which was differentially expressed between the cancerous and normal penile tissues (Table 3), possessed slightly different modification behaviors between the paired specimens as well (S6 and S7 Tables). The distinct modifications including adding or trimming in each 5' or 3' terminus of the deregulated miRNAs expressed in penile cancers and normal adjacent tissues might be related to the carcinogenesis, which findings also provided promising avenues for the design of innovative strategies in fighting against penile cancer.

Conclusions
In summary, to the best of our knowledge, we have obtained the first data for the comprehensive characterizations of miRNAs in penile cancer by next-generation deep sequencing, which technology enabled to precisely quantify expression levels of deregulated miRNAs that may serve as attractively novel diagnostic markers in PeCa. Several differentially expressed miRNAs potentially target networks of genes and signaling pathways probably be involved in the malignant transformation of normal penis to PeCa. Besides, we also explored the novel miRNAs and deregulated piRNAs in cancerous and matched adjacent normal penile tissues. Our data provide an important platform for future investigations aimed to characterize the concrete roles of specific miRNAs and other regulatory smRNAs including piRNAs in penile cancer initiation, progression, invasion and also effective theragnosis. Notably, further studies to conduct the deep sequencing in separated rather than pooled clinical specimens may offer considerable benefits to such efforts.