Integrative miRNA and mRNA analysis in penile carcinomas reveals markers and pathways with potential clinical impact

INTRODUCTION

Penile carcinoma (PeCa) is an aggressive and mutilating disease with high incidence in developing countries, with few therapeutic options available and high morbidity [1–3]. Several risk factors have been identified in PeCa etiology, including poor penile hygiene, phimosis and human papillomavirus (HPV) infection [4, 5]. A recent analysis of cancer registries showed that PeCa patient survival has not improved in Europe nor United States in the last 20 years [6].

The surgical approaches used for the management of PeCa patients are total or partial penectomy, which are associated with high morbidity [7]. The main prognostic factor in PeCa is the presence of lymph node metastasis, which is reported in 20 to 65% of cases [8, 9]. The overall five-year survival in patients with lymph node metastasis is 27% [10]. Novel molecular markers in association with lymph node metastasis [11], disease free survival or overall survival [12–14], recurrence and perineural invasion [15] have been described in PeCa. However, to date, none of them are used in the clinical practice to predict lymph node metastasis and prognosis.

Messenger RNA (mRNA) and micro-RNA (miRNA) expression signatures have diagnostic, prognostic and predictive values in a number of diseases [16]. To our knowledge, two reports described mRNA profiles in PeCa. In 56 PeCa samples, Kroon et al. [17] reported a 44-probe classifier that predicted lymph node metastasis. However, this classifier was not able to predict lymph node metastasis in an independent set of cases. In a previous study by our group, integration of the transcriptome of 33 PeCa samples with DNA methylation profiles [18] revealed several dysregulated oncogenic pathways associated with PeCa development and progression.

To date, five reports have described the differential expression of miRNAs in PeCa. Our group has previously reported that down-expression of SLC8A1, mediated by hsa-miR-223, promotes lower intracellular calcium concentrations, reduced apoptosis and increased cell proliferation in penile tumors [19]. Recently, hsa-miR-218 and miRNA-146a down-expression and EGFR overexpression were reported as associated with high-risk HPV penile tumors [20, 21]. Using next generation sequencing, Zhang et al. [22] identified 56 differentially expressed miRNAs when comparing 10 matched PeCa with adjacent non-cancerous tissues. Hartz et al. [11] reported a miRNA-based signature associated with unfavorable prognosis in 24 PeCa and concomitant inguinal lymph node metastasis. Although these studies have contributed to our understanding of the disease, none have described an integrated analysis using miRNA and mRNA data in matched samples.

Approaches that integrate multiple omic profiles (i.e. mRNA and miRNA) provide a meaningful and comprehensive understanding of the biological processes involved in cancer development and progression. In addition, they may identify the genes that drive tumorigenesis, which have the potential to be applied in translational oncology [23–25].

In this study, integrated miRNA and mRNA profiles from the same set of PeCa samples were investigated, in order to gain insight into the mechanisms of penile carcinogenesis. In addition, miRNA and mRNA data were evaluated according to clinical and pathological characteristics, including lymph node metastasis and HPV infection status.

RESULTS

Distinct miRNA and mRNA expression profiles in PeCa

The unsupervised clustering analysis using miRNA expression data revealed two clusters separating tumors (N = 23) from non-neoplastic penile tissues (NPT = 12) (Supplementary Figure 1). The miRNA profile was composed by 81 differentially expressed miRNAs in PeCa (17 down-expressed and 64 overexpressed, Supplementary Table 1) (P-value < 0.01 and FDR < 5%). hsa-miR-31-5p showed the highest expression levels (FC = 352.4) in tumors compared to normal tissues, while hsa-miR-891a-5p presented the lowest expression levels (FC = –149.3).

The mRNA expression analysis comparing 23 PeCa samples with NPT revealed 2,697 differentially expressed transcripts (947 overexpressed and 1,750 down-expressed, CI = 99.9%, fold-change > 2).

Potential molecular signatures related with prognosis in PeCa were investigated for tumor samples using unsupervised hierarchical clustering analyses for both, miRNA and mRNA profiles. Three main clusters were detected for the PeCa samples according to miRNA profiles (Figure 1A). Cluster 3 (eight cases) was mainly enriched by the tumors with aggressive features (five cases with T3-T4 tumor size and five with lymph node metastasis). A similar analysis was performed with mRNA profiles (Figure 1B), which also revealed three clusters, where cluster 3 (seven cases) was comprised of patients with poor prognosis (six cases with T3-T4 tumor size, five with lymph node metastasis, and four with perineural invasion). Although not statistically significant, cluster 3 for both miRNA (P = 0.39) and mRNA (P = 0.49) profiles were enriched with patients with shorter overall survival (Supplementary Figure 2).

Integration of miRNA and mRNA expression profiles reveals potential disrupted pathways in PeCa

Integrative analysis was performed using 81 differentially expressed miRNAs and 2,697 mRNAs. Based on predicted (mirWalk 2.0) and/or experimentally validated interactions (miRTarBase), 68 miRNAs that potentially regulate 255 mRNAs were identified (Supplementary Table 2), representing 598 miRNA/mRNA interactions with negative correlation (r Spearman < 0) and inverted fold change (Supplementary Figure 3).

The main canonical pathways detected by IPA and confirmed by KOBAS 2.0 revealed enrichment of the Human Embryonic Stem Cell Pluripotency, VEGF signaling, Molecular Mechanisms of Cancer, B Cell Receptor, PDGF, ERBB, Matrix Metalloproteases and PI3K/AKT signaling pathways involving transcripts detected in the integrative analysis (Figure 1C, Table 1).

miRNA and mRNA validation

Eight miRNAs and 10 transcripts were evaluated by RT-qPCR analysis in the same set of samples used for microarray and in an additional group of tumors (validation set) (Figure 2; Supplementary Table 3). The candidates for validation were selected according to the following criteria: (a) negative correlation between the mRNA and its miRNA regulator detected in the integrative analysis, (b) prognostic association, (c) high FC, (d) low FDR, and (e) low P-value in the miRNA and mRNA microarray analysis. Overexpression of hsa-miR-20a-5p, hsa-miR-29b-3p, hsa-miR-31-5p, hsa-miR-224-5p, hsa-miR-106a-5p, hsa-miR-17-5p, hsa-miR-223-3p and down-expression of hsa-miR-145-5p were confirmed by RT-qPCR analysis in the validation samples set (Figure 2A). All, with the exception of hsa-miR-20a-5p and hsa-miR-29b-3p, were also confirmed as significantly dysregulated in the same set of samples used for microarray (Supplementary Table 3).

Significant down-expression of AR, ERBB4, FGFR1, NRN1 and PPARG and overexpression of MMP1, MMP12, NRAS and SPP1 were confirmed by RT-qPCR assays in the validation set of samples (Figure 2B; Supplementary Table 3). Although not significant, DNMT3A presented decreased expression in the tumor samples.

In order to optimize the accuracy to distinguish tumors from non-neoplastic tissues, a higher discriminating power for the assessed markers was prioritized over the miRNA/mRNA interaction results. Three miRNAs and three mRNAs, presenting higher AUC (AUC for hsa-miR-31-5p = 0.861, hsa-miR-224-5p = 0.739 and hsa-miR-223-3p = 0.733, MMP1 = 0.923, MMP12 = 0.865 and PPARG = 0.851) were selected to construct a molecular classifier. miRNA and mRNA classifier accuracy was 79% and 89%, respectively. Sensitivity and specificity were 82% and 74% for miRNA classifier and 92% and 83% for mRNA classifier, respectively (Supplementary Figure 4).

MMP1 is a potential prognostic marker for lymph node metastasis in PeCa

The integrative analysis results were compared with clinical and pathological data including lymph node metastasis, tumor size and perineural invasion (Supplementary Table 4). Statistically significant associations were identified for lymph node metastasis (12 genes and three miRNAs), tumor size (eight genes and one miRNA) and perineural invasion (UHMK1 gene) (P < 0.01, FDR < 20%). Among these genes, increased expression of MMP1 was observed in cases with lymph node metastasis, further confirmed by RT-qPCR in the validation set of samples (Figure 3A).

MMP1 expression was a better predictor factor of lymph node metastasis than well established pathological parameters (AUC: MMP1 = 0.751; histological grade = 0.672; perineural invasion = 0.596; tumor size = 0.376) (Figure 3B, Supplementary Table 5). Higher MMP1 expression levels (in both microarray and RT-qPCR analysis) were associated with shorter survival, although not significant (Figure 3C). The hsa-miR-145, which regulates MMP1, was down-expressed compared with NPT. However, no association with this miRNA and pathological characteristics was observed.

Differentially expressed miRNA and mRNA according to HPV status

Seven of 23 PeCa samples (30.4%) presented high-risk HPV genotype (6 HPV16 cases and 1 HPV18). The mRNA hierarchical clustering analysis revealed four of seven HPV-positive cases grouped together in cluster 1. In addition, the miRNA analysis presented the cluster 2 enriched with five of seven HPV-positive tumors. The cluster 3, which comprised patients with poor prognosis, showed a limited number of HPV-positive cases (2 for mRNA and 1 for miRNA cluster analysis) (Figure 1A and 1B).

The integrative analysis revealed five miRNAs (hsa-let-7b-5p, hsa-miR-146a-5p, hsa-miR-185-5p, hsa-miR-29b-3p and hsa-miR-505-3p) down-expressed in HPV positive cases compared with negative cases. Thirteen transcripts (down-expression of CSF1 and PKD2 and overexpression of PPM1B, INPP5A, LONRF1, WASF3, PRKG1, NTF3, NBEA, EGR1, RGS5, NTRK2 and OLFM1) were also detected as dysregulated in HPV positive cases (FDR > 20%) (Supplementary Table 4).

DISCUSSION

Dysregulation of the transcriptome and miRNA machinery is a common process in cancer development and progression. In the present study, gene expression and miRNA profiles were able to distinguish one cluster (cluster 3), which contained the majority of cases presenting aggressive clinicopathological features (lymph node metastasis and T3-T4 tumor size).

Integrative analysis of miRNA and mRNA is an important tool to identify potential diagnostic and prognostic markers in a variety of tumors [23, 24]. Recently, Zhang et al. [22] described 384 differentially expressed miRNAs with potential involvement in penile malignant transformation. Comparison of the differentially expressed miRNAs found in the present study with those described by Zhang et al. [22] revealed a 15% overlap. A plausible explanation for this is the use of different methods in the two studies (next generation sequencing of a pool of 10 samples versus the Taqman Low density array with a larger number of cases used in our study).

Using in silico molecular analysis, we described 598 interactions involving 68 miRNAs and 255 mRNAs. hsa-miRNA-31-5p was identified as having the highest fold change (FC = 352.4) in tumor samples, as previously described in several epithelial tumors (head and neck squamous cell, colorectal, prostate and lung cancer) [26, 27], as well as in premalignant colonic lesions [27]. In addition, this miRNA in combination with hsa-miR-224-5p and hsa-miR-223-3p were able to distinguish tumors from NPT, indicating its importance in PeCa development. A variety of experimentally validated target genes regulated by hsa-miR-31-5p has been described in different tumor types [26], including AR and DNMT3A genes. The androgen receptor (AR) is a critical transcription factor with prognostic value in breast cancer and prostate [28, 29] as well as in bladder carcinomas [30].

The expression of AR was significantly decreased in penile tumors, both in the microarray and in the set of cases used for validation (p < 0.001). Interestingly, the ESR1 mRNA was also down-expressed in PeCa, with several miRNAs that potentially regulate this gene being overexpressed (Supplementary Table 2). In addition, AR and ESR1 presented the highest number of predicted functional interactions with other genes (Figure 1C), suggesting their potential as molecular drivers in PeCa. Functional loss of the AR gene by promoter hypermethylation has been described as involved in PeCa development [31]. A previous study by our group on penile tumors did not reveal AR hypermethylation [18]. In addition to hypermethylation, dysregulation of the miRNAs that target AR may also contribute to PeCa development and progression. We found four overexpressed miRNAs predicted as involved in the regulation of the AR gene (hsa-miR-31-5p, hsa-miR-34a-5p, hsa-miR-205-5p and hsa-miR-185-5p) and possibly associated with its down-expression.

Potential dysregulation of the VEGF signaling pathway was also observed by in silico analysis. Genes involved in this pathway have already been shown to be regulated by miRNAs [32]. The hsa-miR-424-5p overexpression, which has previously been associated with angiogenesis [33], was described in our study. This miRNA regulates the FGFR1 gene, confirmed by RT-qPCR as down-expressed. A clinical trial to test anti-angiogenic therapy in PeCa patients (NCT02279576) started in 2014. Overall, our findings give additional support to the relevance of genes and miRNAs associated with angiogenesis in PeCa.

While the role of VEGF and angiogenesis remains unexplored in PeCa, the involvement of the EGF pathway has been reported. Recently, EGFR gene mutation, amplification and overexpression have been described in PeCa, suggesting the potential of using therapeutic strategies targeting EGFR in a selected group of patients [34, 35]. Currently, EGFR-targeted therapies are found to be clinically useful in PeCa patients (approximately 30 cases) [34–39]. Although the EGFR expression levels were not identified in this study, four ligands of EGFR (AREG, EREG, TGFA and EPGN) were overexpressed, as well as the downstream effector NRAS. These data suggest that this pathway is dysregulated and, hence contribute to penile carcinogenesis.

Previously, we have shown the association of PPARG loss with poor prognostic features in PeCa (advanced clinical and T stage and lymph node metastasis) [12]. We have also demonstrated that PPARG loss leads to its down-expression as a result of a gene dosage effect. Increasing PPARG activation using agonists has been shown as a powerful strategy to either inhibit cell proliferation or induce apoptosis [40]. The use of PPARG gene as a therapeutic target is already being applied in a considerable number of tumors, including breast, prostate, colorectal and thyroid [40].

The metalloprotease genes family was shown to have the highest fold-changes, including MMP1 (FC = 51.1), MMP10 (FC = 52.5) and MMP12 (FC = 89.9). In the validation set of PeCa samples, overexpression of MMP1 and MMP12 was confirmed by RT-qPCR. In addition, a classifier with MMP1 and MMP12 in combination with PPARG was able to discriminate PeCa from NPT. Molecular alterations have been reported as preceding morphological and pathological changes [41, 42]. Therefore, markers that could indicate regions with genetic alteration may be of great value to define margin assessment. Along with the fact that matrix metalloproteinase was one of the most important pathways dysregulated in PeCa, these results suggest that MMPs genes might act as oncogenic drivers for PeCa development.

Integrative analysis revealed that MMP1 and MMP12 may be regulated by hsa-miR-145-5p, which was down-expressed in PeCa. Similarly, Zhang et al. [22] reported down-expression of this miRNA in 10 PeCa. Down-expression of hsa-miR-145 has been described in several cancer types, mediating suppression of cell growth, invasion and metastasis [43]. In prostate cancer, hsa-miR-145 down-expression has been reported as a predictor of poor prognosis [44].

Although hsa-miR-145 down-expression was unable to predict poor prognosis in our cohort of cases, its target MMP1 showed increased expression levels in patients with lymph node metastasis. MMP1 gene expression levels were a better predictor of lymph node metastasis than tumor size, histological grade and perineural invasion. In agreement with this finding, Zheng et al. [45] showed that hsa-miR-145 overexpression altered MMP1 and MMP9 mRNA and protein levels, with subsequent inhibition of invasion, metastasis and angiogenesis in gastric cancer cells. Although some prognostic molecular makers have been reported in PeCa, the stratification of the patients is currently based on clinical and histopathological features [11, 31, 39, 46, 47]. In the present study, MMP1 showed superior performance in discriminate lymph node metastasis in PeCa compared with established clinical-pathological parameters. The Supplementary Table 4 shows a list of transcripts and miRNA potentially associated with prognosis and clinical features in PeCa (P < 0.01 and FDR < 20%). However, none of them were confirmed as predictors of lymph node metastasis (data not shown).

Although controversial, HPV infection has been reported as a risk factor involved in the PeCa etiology influencing the prognosis [5]. Gregoire et al. [48] reported an association between HPV positive infection and poor prognosis in penile cancer patients. However, HPV infection was also reported as conferring a survival advantage in these patients [49, 50]. Lont et al. [50] reported a 5-year cancer-specific survival rate of 93% for HPV-positive and 78% for HPV-negative cases. The cluster 3 (enriched with cases showing poor prognosis features) of our miRNA and mRNA profiles revealed a limited number of HPV positive cases, which suggests a less aggressive phenotype. In addition, a set of five miRNAs was down-expressed in HPV positive cases, including the hsa-let-7b-5p, hsa-miR-146a-5p, hsa-miR-185-5p, hsa-miR-29b-3p and hsa-miR-505-3p. Recently, Peta et al. [21] reported lower expression levels of hsa-miR-146a in high-risk HPV positive than in negative tumors.

In conclusion, our integrative analysis was able to identify miRNA and mRNA related with cancer development and progression. Furthermore, MMP1 is a predictive marker of lymph node metastasis in PeCa. We also pinpointed PPARG, VEGF, EGFR and matrix metalloproteinase pathways as dysregulated in PeCa samples, endorsing their involvement as potential targets for PeCa treatment.

MATERIALS AND METHODS

Patients and sample collection

A total of 101 samples was included in this study: 59 penile carcinomas (PeCa), 26 surrounding normal tissues (SNT) and 16 normal glans (NG) (Supplementary Figure 5). Fresh frozen PeCa tissue and SNT were obtained from untreated patients who underwent tumor resection at A.C. Camargo Cancer Center, São Paulo, Barretos Cancer Hospital, Barretos, and Faculty of Medicine, Botucatu, SP, Brazil. Normal glans samples were obtained from autopsies. The majority of the penile cancer tissue samples were confirmed histologically as usual penile carcinomas. All samples used in this study, including SNT, were submitted to cellular macrodissection and histology confirmation. Tumor samples presented at least 80% of tumor cells and SNT were composed by normal epithelial cells. Written informed consent was obtained from all patients or relatives. The Human Research Ethic Committees of the aforementioned Institutions approved the study (Protocols #1230/09: A.C. Camargo Cancer Center; #363–2010: Barretos Cancer Hospital, and #501.229/2013: Faculty of Medicine, Botucatu, SP, Brazil).

Previously, we reported a transcriptomic analysis (4 × 44K, Agilent Technologies, Santa Clara, CA, USA) of 33 PeCa and a pool of four NG [18]. Twenty-three PeCa samples had tissue available for miRNA expression analysis. Seven SNT (paired with 7 PeCa) and five NG samples were also included in the miRNA expression analysis as control. Integrative analysis of mRNA and miRNA expression data was performed for 23 PeCa samples.

Quantitative RT-PCR was applied in the same group of samples used in the array experiments (21 PeCa, 6 SNT and 5 NG for miRNA microarray analysis and 20 PeCa and 3 NG for mRNA microarray analysis) and in the validation set of samples (33 PeCa, 20 SNT and 7 NG for miRNA and 36 PeCa, 9 SNT and 10 NG for mRNA analysis). Non-neoplastic penile tissues (NPT) composed by SNT and NG samples were compared with tumor tissue. Clinical and histopathological data for the PeCa samples are shown in Table 2.

Human papilloma virus (HPV) status was investigated for all PeCa samples (Linear Array HPV Test Genotyping, Roche Molecular Diagnostics, Branchburg, NJ, USA). HPV-positive cases were detected in seven of the 23 cases (30.4%) used for microarray analysis and in nine of the 36 PeCa (25%) in the validation set (Table 2).

miRNA and mRNA profiles in PeCa

Total RNA was obtained from macrodissected fresh frozen tissues using miRNAeasy Kit (Qiagen, Venlo, Limburg, Netherlands). TaqMan Human MicroRNA Assay System Set v2.0 (Applied Biosystems, Foster City, CA, USA) was used for miRNA expression analysis. A set of pre-defined primers (Megaplex RT primers™, Pool A, Applied Biosystems, Foster City, CA, USA) was used for cDNA synthesis as recommended by the manufacturer. Data was normalized using the Pfaffl model [51] with MammU6, RNU44 and RNU48 as references. Low abundant miRNAs (undetermined quantification cycle) that were observed in more than 20% of samples from each comparison were excluded from the analysis. Biological groups were compared using the two-sample t-test with BRB ArrayTools software (v. 4.4.0) [52], establishing a two-tailed P-value of < 0.01 with a low false discovery rate ratio (FDR< 5%) and at least two fold changes (FC > 2) as significant.

mRNA expression analysis was performed using the Whole Human Genome 4x44K microarray platform (Agilent Technologies) as described by Kuasne et al. [18] Transcriptomic data are available on the Gene Expression Omnibus (GEO) database (GSE57955). Unsupervised hierarchical clustering analysis of miRNA and mRNA expression was accomplished using one minus correlation metric and complete linkage. The Kaplan Meier curve and log rank test were undertaken to estimate the overall survival [53].

Integrative analysis

Integrative analysis using miRNA and mRNA expression data from 23 tumor samples was based on predicted and experimentally validated miRNA/mRNA interactions. Predicted miRNA/mRNA interactions were performed using miRWalk 2.0 [54], in at least 10 of 12 target prediction tools (Last access: October 2015). In addition, experimentally validated miRNA/mRNA interactions by reporter assays were obtained from miRTarBase database [55] (Last access: October 2015). Subsequently, negative miRNA/target mRNA correlation (r < 0) presenting inverted FC was considered.

Transcripts regulated by miRNAs were submitted to in silico analysis using Ingenuity Pathways Analysis (IPA, Ingenuity^® System) and KEGG Orthology Based Annotation System (KOBAS 2.0) software [56]. Protein-protein interactions were assessed with known and predicted physical interactions using I2D version 2.0 (http://ophid.utoronto.ca/i2d). The resulting networks were visualized using NAViGaTOR version 2.3 [57]. http://ophid.utoronto.ca/navigator. Clinical associations were performed using two-sample t-test (P < 0.01, FDR < 20%).

Evaluation of differentially expressed miRNAs and mRNA by reverse transcription quantitative PCR (RT-qPCR)

Eight miRNAs (hsa-miR-20a-5p, hsa-miR-29b-3p, hsa-miR-31-5p, hsa-miR-224-5p, hsa-miR-106a-5p, hsa-miR-17-5p, hsa-miR-223-3p, and hsa-miR-145-5p) and ten mRNAs (AR, DNMT3A, ERBB4, FGFR1, MMP1, MMP12, NRAS, NRN1, PPARG, and SPP1) were evaluated by RT-qPCR. The candidates for validation were selected according to the following criteria (in descending order of importance): (a) miRNA-mRNA pair detected in the integrative analysis with negative correlation, (b) prognostic association, (c) high FC, (d) low FDR, and (e) low P-value in the individual miRNA and mRNA analysis.

The RT-qPCR experiments followed the MIQE guideline recommendations [58]. RNU48 and GUSB were used as references for miRNAs and mRNA, respectively, as previously reported [19]. cDNA synthesis was performed using total RNA and the TaqMan miRNA Reverse Transcription Kit (Applied Biosystems, Foster City, CA, USA), according to the manufacturer’s recommendations. MicroRNA expression was assessed using the TaqMan^® MicroRNA Assay (Applied Biosystems, Foster City, CA, USA). For mRNA expression analysis, cDNA synthesis was conducted as previously described [59]. The reactions were carried out by automated pipetting (QIAgility, Qiagen, Courtaboeuf, France) in duplicate using TaqMan^® Universal PCR Master Mix, No AmpErase^® UNG (Applied Biosystems, Foster City, CA, USA) (miRNA) or Syber Green Master Mix (mRNA) using the 7900 Real time PCR System (Applied Biosystems, Foster City, CA).

Relative quantification of the expression levels was calculated according to Pfaffl method [51]. Parametric t test was applied to compare tumors with non-neoplastic penile tissue (NPT) (comprising SNT and NG samples) and according to clinicopathological features. Statistical analysis was performed using GraphPad Prism5 (San Diego, CA, USA) and SPSS version 21.0 (SPSS, Chicago, IL, USA).

Classifiers were designed to distinguish PeCa from NPT, using the markers presenting the higher area under the ROC (Receiver Operating Characteristic) curve (AUC). The support vector machine (SVM) method was applied and performance was assessed by leave-one-out-cross-validation (LOOCV) using BRB ArrayTools software (v. 4.4.0). A similar approach was carried out to predict lymph node status. Overall survival analysis was performed using Kaplan-Meier and log rank test. In this analysis, gene expression was dichotomized as low (bellow median) and high (above median) expression.

ACKNOWLEDGMENTS

The authors would like to acknowledge the A.C. Camargo Cancer Center biobank, SP and Barretos Cancer Hospital, SP, for providing human specimens. We are also grateful to the Pathology Departments of A. C. Camargo Cancer Center, SP, Barretos Cancer Hospital, SP and Botucatu Medical School, SP, Brazil.

CONFLICTS OF INTEREST

The authors declare no conflicts of interest.

GRANT SUPPORT

The study was supported by São Paulo Research Foundation (FAPESP 2009/52088-3, 2010/51601-6 and fellowship 2013/03667-6) and National Council for Scientific and Technological Development (CNPq).

Authors’ contributions

Conceived and designed the experiments: HK, MCBF, ABL, MP, JJMM, FAM, SRR. Clinical care and patients’ treatment: EFF, GCG, AL, JCSTF, Histopathological analysis: CSN. Supervised the study: SRR. Wrote the manuscript: HK, MCBF, SAD, SRR. All authors read and approved the final manuscript.

REFERENCES

1. Zukiwskyj M, Daly P, Chung E. Penile cancer and phallus preservation strategies: a review of current literature. BJU Int. 2013; 112:21–6.

2. Bleeker MC, Heideman DA, Snijders PJ, Horenblas S, Dillner J, Meijer CJ. Penile cancer: epidemiology, pathogenesis and prevention. World J Urol. 2009; 27:141–50.

3. Guimarães GC, Rocha RM, Zequi SC, Cunha IW, Soares FA. Penile cancer: epidemiology and treatment. Curr Oncol Rep. 2011; 13:231–9.

4. Steinestel J, Al Ghazal A, Arndt A, Schnoeller TJ, Schrader AJ, Moeller P, Steinestel K. The role of histologic subtype, p16(INK4a) expression, and presence of human papillomavirus DNA in penile squamous cell carcinoma. BMC Cancer. 2015; 15:220.

5. Chaux A, Cubilla AL. The role of human papillomavirus infection in the pathogenesis of penile squamous cell carcinomas. Semin Diagn Pathol. 2012; 29:67–71.

6. Verhoeven RH, Janssen-Heijnen ML, Saum KU, Zanetti R, Caldarella A, Holleczek B, Brewster DH, Hakulinen T, Horenblas S, Brenner H, Gondos A, Group ESW. Population-based survival of penile cancer patients in Europe and the United States of America: no improvement since 1990. Eur J Cancer. 2013; 49:1414–21.

7. Mossanen M; Holt S, Gore JL, Lin DW, Wright MD. 15 Years of penile cancer management in the United States: An analysis of the use of partial penectomy for localized disease and chemotherapy in the mestastatic setting. Urologic Seminars: Seminars and Original Investigations. 2016; pii: S1078–1439 (16) 30162–4.

8. Ficarra V, Akduman B, Bouchot O, Palou J, Tobias-Machado M. Prognostic factors in penile cancer. Urology. 2010; 76:S66–73.

9. Chalya PL, Rambau PF, Masalu N, Simbila S. Ten-year surgical experiences with penile cancer at a tertiary care hospital in northwestern Tanzania: a retrospective study of 236 patients. World J Surg Oncol. 2015; 13:71.

10. Barski D, Georgas E, Gerullis H, Ecke T. Metastatic penile carcinoma - an update on the current diagnosis and treatment options. Cent European J Urol. 2014; 67:126–32.

11. Hartz JM, Engelmann D, Fürst K, Marquardt S, Spitschak A, Goody D, Protzel C, Hakenberg OW, Pützer BM. Integrated loss of miR-1/-101/-204 discriminates metastatic from non-metastatic penile carcinomas and can predict patient outcome. J Urol. 2016; 196:570–8.

12. Busso-Lopes AF, Marchi FA, Kuasne H, Scapulatempo-Neto C, Trindade-Filho JC, de Jesus CM, Lopes A, Guimarães GC, Rogatto SR. Genomic profiling of human penile carcinoma predicts worse prognosis and survival. Cancer Prev Res (Phila). 2015; 8:149–56.

13. Zargar-Shoshtari K, Spiess PE, Berglund AE, Sharma P, Powsang JM, Giuliano A, Magliocco AM, Dhillon J. Clinical Significance of p53 and p16(ink4a) Status in a Contemporary North American Penile Carcinoma Cohort. Clin Genitourin Cancer. 2016; 14:346–51.

14. Rocha RM, Ignácio JA, Jordán J, Carraro DM, Lisboa B, Lopes A, Carvalho KC, da Cunha IW, Cubilla A, Guimarães GC, Vassallo J, Soares FA. A clinical, pathologic, and molecular study of p53 and murine double minute 2 in penile carcinogenesis and its relation to prognosis. Hum Pathol. 2012; 43:481–8.

15. da Silva Amancio AM, da Cunha IW, Neves JI, da Silva Quetz J, Carraro DM, Rocha RM, Zequi SC, Cubilla AL, da Fonseca FP, Lopes A, do Perpétuo Socorro Saldanha da Cunha M, Lima MV, Vassallo J, et al. Epidermal growth factor receptor as an adverse survival predictor in squamous cell carcinoma of the penis. Hum Pathol. 2016

16. Iorio MV, Croce CM. MicroRNA dysregulation in cancer: diagnostics, monitoring and therapeutics. A comprehensive review. EMBO Mol Med. 2012; 4:143–59.

17. Kroon BK, Leijte JA, van Boven H, Wessels LF, Velds A, Horenblas S, van’t Veer LJ. Microarray gene-expression profiling to predict lymph node metastasis in penile carcinoma. BJU Int. 2008; 102:510–5.

18. Kuasne H, Cólus IM, Busso AF, Hernandez-Vargas H, Barros-Filho MC, Marchi FA, Scapulatempo-Neto C, Faria EF, Lopes A, Guimarães GC, Herceg Z, Rogatto SR. Genome-wide methylation and transcriptome analysis in penile carcinoma: uncovering new molecular markers. Clin Epigenetics. 2015; 7:46.

19. Muñoz JJ, Drigo SA, Barros-Filho MC, Marchi FA, Scapulatempo-Neto C, Pessoa GS, Guimarães GC, Trindade Filho JC, Lopes A, Arruda MA, Rogatto SR. Down-Regulation of SLC8A1 as a Putative Apoptosis Evasion Mechanism by Modulation of Calcium Levels in Penile Carcinoma. J Urol. 2015; 194:245–51.

20. Barzon L, Cappellesso R, Peta E, Militello V, Sinigaglia A, Fassan M, Simonato F, Guzzardo V, Ventura L, Blandamura S, Gardiman M, Palù G, Fassina A. Profiling of expression of human papillomavirus-related cancer miRNAs in penile squamous cell carcinomas. Am J Pathol. 2014; 184:3376–83.

21. Peta E, Cappellesso R, Masi G, Sinigaglia A, Trevisan M, Grassi A, Di Camillo B, Vassarotto E, Fassina A, Palù G, Barzon L. Down-regulation of microRNA-146a is associated with high-risk human papillomavirus infection and epidermal growth factor receptor overexpression in penile squamous cell carcinoma. Hum Pathol. 2016; S0046–8177(16)30290–8.

22. Zhang L, Wei P, Shen X, Zhang Y, Xu B, Zhou J, Fan S, Hao Z, Shi H, Zhang X, Kong R, Xu L, Gao J, et al. MicroRNA Expression Profile in Penile Cancer Revealed by Next-Generation Small RNA Sequencing. PLoS One. 2015; 10:e0131336.

23. Akavia UD, Litvin O, Kim J, Sanchez-Garcia F, Kotliar D, Causton HC, Pochanard P, Mozes E, Garraway LA, Pe’er D. An integrated approach to uncover drivers of cancer. Cell. 2010; 143:1005–17.

24. Cirilo PD, Marchi FA, Barros Filho MeC, Rocha RM, Domingues MA, Jurisica I, Pontes A, Rogatto SR. An integrative genomic and transcriptomic analysis reveals potential targets associated with cell proliferation in uterine leiomyomas. PLoS One. 2013; 8:e57901.

25. Chari R, Thu KL, Wilson IM, Lockwood WW, Lonergan KM, Coe BP, Malloff CA, Gazdar AF, Lam S, Garnis C, MacAulay CE, Alvarez CE, Lam WL. Integrating the multiple dimensions of genomic and epigenomic landscapes of cancer. Cancer Metastasis Rev. 2010; 29:73–93.

26. Laurila EM, Kallioniemi A. The diverse role of miR-31 in regulating cancer associated phenotypes. Genes Chromosomes Cancer. 2013; 52:1103–13.

27. Oberg AL, French AJ, Sarver AL, Subramanian S, Morlan BW, Riska SM, Borralho PM, Cunningham JM, Boardman LA, Wang L, Smyrk TC, Asmann Y, Steer CJ, et al. miRNA expression in colon polyps provides evidence for a multihit model of colon cancer. PLoS One. 2011; 6:e20465.

28. Vera-Badillo FE, Templeton AJ, de Gouveia P, Diaz-Padilla I, Bedard PL, Al-Mubarak M, Seruga B, Tannock IF, Ocana A, Amir E. Androgen receptor expression and outcomes in early breast cancer: a systematic review and meta-analysis. J Natl Cancer Inst. 2014; 106:djt319.

29. Shtivelman E, Beer TM, Evans CP. Molecular pathways and targets in prostate cancer. Oncotarget. 2014; 5:7217–59. doi: 10.18632/oncotarget.2406.

30. Lombard AP, Mudryj M. The emerging role of the androgen receptor in bladder cancer. Endocr Relat Cancer. 2015; 22:R265–77.

31. Feber A, Arya M, de Winter P, Saqib M, Nigam R, Malone PR, Tan WS, Rodney S, Lechner M, Freeman A, Jameson C, Muneer A, Beck S, et al. Epigenetics markers of metastasis and HPV-induced tumorigenesis in penile cancer. Clin Cancer Res. 2015; 21:1196–206.

32. Landskroner-Eiger S, Moneke I, Sessa WC. miRNAs as modulators of angiogenesis. Cold Spring Harb Perspect Med. 2013; 3:a006643.

33. Ghosh G, Subramanian IV, Adhikari N, Zhang X, Joshi HP, Basi D, Chandrashekhar YS, Hall JL, Roy S, Zeng Y, Ramakrishnan S. Hypoxia-induced microRNA-424 expression in human endothelial cells regulates HIF-α isoforms and promotes angiogenesis. J Clin Invest. 2010; 120:4141–54.

34. McDaniel AS, Hovelson DH, Cani AK, Liu CJ, Zhai Y, Zhang Y, Weizer AZ, Mehra R, Feng FY, Alva AS, Morgan TM, Montgomery JS, Siddiqui J, et al. Genomic Profiling of Penile Squamous Cell Carcinoma Reveals New Opportunities for Targeted Therapy. Cancer Res. 2015; 75:5219–27.

35. Ali SM, Pal SK, Wang K, Palma NA, Sanford E, Bailey M, He J, Elvin JA, Chmielecki J, Squillace R, Dow E, Morosini D, Buell J, et al. Comprehensive Genomic Profiling of Advanced Penile Carcinoma Suggests a High Frequency of Clinically Relevant Genomic Alterations. Oncologist. 2016; 21:33–9.

36. Brown A, Ma Y, Danenberg K, Schuckman AK, Pinski JK, Pagliaro LC, Quinn DI, Dorff TB. Epidermal growth factor receptor-targeted therapy in squamous cell carcinoma of the penis: a report of 3 cases. Urology. 2014; 83:159–65.

37. Men HT, Gou HF, Qiu M et a. A case of penile squamous cell carcinoma treated with a combination of antiepidermal growth factor receptor antibody and chemotherapy. Anticancer Drugs. 2014; 25:123–5.

38. Carthon BC, Ng CS, Pettaway CA, Pagliaro LC. Epidermal growth factor receptor-targeted therapy in locally advanced or metastatic squamous cell carcinoma of the penis. BJU Int. 2014; 113:871–7.

39. Di Lorenzo G, Buonerba C, Ferro M, Calderoni G, Bozza G, Federico P, Tedesco B, Ruggieri V, Aieta M. The epidermal growth factor receptors as biological targets in penile cancer. Expert Opin Biol Ther. 2015; 15:473–6.

40. Fröhlich E, Wahl R. Chemotherapy and chemoprevention by thiazolidinediones. Biomed Res Int. 2015; 2015:845340.

41. Magi-Galluzzi C, Maddala T, Falzarano SM, Cherbavaz DB, Zhang N, Knezevic D, Febbo PG, Lee M, Lawrence HJ, Klein EA. Gene expression in normal-appearingtissue adjacent to prostate cancers are predictive of clinical outcome: evidence for a biologically meaningful field effect. Oncotarget. 2016; 7:33855–65. doi: 10.18632/oncotarget.8944.

42. Baba Y, Ishimoto T, Kurashige J, Iwatsuki M, Sakamoto Y, Yoshida N, Watanabe M, Baba H. Epigenetic field cancerization in gastrointestinal cancers. Cancer Lett. 2016; 375:360–6.

43. Sachdeva M, Mo YY. miR-145-mediated suppression of cell growth, invasion and metastasis. Am J Transl Res. 2010; 2:170–80.

44. Zhang X, Wu J. Prognostic role of microRNA-145 in prostate cancer: A systems review and meta-analysis. Prostate Int. 2015; 3:71–4.

45. Zheng L, Pu J, Qi T, Qi M, Li D, Xiang X, Huang K, Tong Q. miRNA-145 targets v-ets erythroblastosis virus E26 oncogene homolog 1 to suppress the invasion, metastasis, and angiogenesis of gastric cancer cells. Mol Cancer Res. 2013; 11:182–93.

46. Gunia S, Kakies C, Erbersdobler A, Hakenberg OW, Koch S, May M. Expression of p53, p21 and cyclin D1 in penile cancer: p53 predicts poor prognosis. J Clin Pathol. 2012; 65:232–6.

47. Ferrandiz-Pulido C, Masferrer E, Toll A, Hernandez-Losa J, Mojal S, Pujol RM, Ramon y Cajal S, de Torres I, Garcia-Patos V. mTOR signaling pathway in penile squamous cell carcinoma: pmTOR and peIF4E over expression correlate with aggressive tumor behavior. J Urol. 2013;190:2288–95.

48. Gregoire L, Cubilla AL, Reuter VE, Haas GP, Lancaster WD. Preferential association of human papillomavirus with high-grade histologic variants of penile-invasive squamous cell carcinoma. J Natl Cancer Inst. 1995; 87:1705–9.

49. Djajadiningrat RS, Jordanova ES, Kroon BK, van Werkhoven E, de Jong J, Pronk DT, Snijders PJ, Horenblas S, Heideman DA. Human papillomavirus prevalence in invasive penile cancer and association with clinical outcome. J Urol. 2015 Feb;193:526–31.

50. Lont AP, Kroon BK, Horenblas S, Gallee MP, Berkhof J, Meijer CJ, Snijders PJ. Presence of high-risk human papillomavirus DNA in penile carcinoma predicts favorable outcome in survival. Int J Cancer. 2006; 119:1078–81.

51. Pfaffl MW. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001; 29:e45.

52. Wright GW, Simon RM. A random variance model for detection of differential gene expression in small microarray experiments. Bioinformatics. 2003; 19:2448–55.

53. Efron B. Logistic Regression, Survival Analysis, and the Kaplan-Meier Curve. In: Taylor & Francis L. Journal of the American Statistical Associationed. 1988; 83:12.

54. Dweep H, Gretz N, Sticht C. miRWalk database for miRNA-target interactions. Methods Mol Biol. 2014; 1182:289–305.

55. Chou CH, Chang NW, Shrestha S, Hsu SD, Lin YL, Lee WH, Yang CD, Hong HC, Wei TY, Tu SJ, Tsai TR, Ho SY, Jian TY, et al. miRTarBase 2016: updates to the experimentally validated miRNA-target interactions database. Nucleic Acids Res. 2015.

56. Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, Kong L, Gao G, Li CY, Wei L. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011; 39:W316–22.

57. Brown KR, Otasek D, Ali M, McGuffin MJ, Xie W, Devani B, Toch IL, Jurisica I. NAViGaTOR: Network Analysis, Visualization and Graphing Toronto. Bioinformatics. 2009; 25:3327–9.

58. Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M, Mueller R, Nolan T, Pfaffl MW, Shipley GL, Vandesompele J, Wittwer CT. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin Chem. 2009; 55:611–22.

59. Rosa FE, Caldeira JR, Felipes J, Bertonha FB, Quevedo FC, Domingues MA, Moraes Neto FA, Rogatto SR. Evaluation of estrogen receptor alpha and beta and progesterone receptor expression and correlation with clinicopathologic factors and proliferative marker Ki-67 in breast cancers. Hum Pathol. 2008; 39:720–30.