Identification of candidate miRNAs in early-onset and late-onset prostate cancer by network analysis

The incidence of patients under 55 years old diagnosed with Prostate Cancer (EO-PCa) has increased during recent years. The molecular biology of PCa cancer in this group of patients remains unclear. Here, we applied weighted gene coexpression network analysis of the expression of miRNAs from 24 EO-PCa patients (38–45 years) and 25 late-onset PCa patients (LO-PCa, 71–74 years) to identify key miRNAs in EO-PCa patients. In total, 69 differentially expressed miRNAs were identified. Specifically, 26 and 14 miRNAs were exclusively deregulated in young and elderly patients, respectively, and 29 miRNAs were shared. We identified 20 hub miRNAs for the network built for EO-PCa. Six of these hub miRNAs exhibited prognostic significance in relapse‐free or overall survival. Additionally, two of the hub miRNAs were coexpressed with mRNAs of genes previously identified as deregulated in EO-PCa and in the most aggressive forms of PCa in African-American patients compared with Caucasian patients. These genes are involved in activation of immune response pathways, increased rates of metastasis and poor prognosis in PCa patients. In conclusion, our analysis identified miRNAs that are potentially important in the molecular pathology of EO-PCa. These genes may serve as biomarkers in EO-PCa and as possible therapeutic targets.

Total RNA was extracted from the primary tumor tissue and matched control normal tissue samples, which were obtained from formalin-fixed paraffin-embedded tissue blocks from prostatectomies. The small RNA profile was generated using the Illumina Human Whole-Genome DASL (cDNA-mediated annealing, selection, extension, and ligation), while the miRNAs were sequenced on the Illumina HuSeq 2500 platform.
Identification of differentially expressed miRNAs. The comparison between transcriptomes of tumor and normal prostate samples employed stringent criteria of a fold change (FC) greater than 2 and less than -2 and a false discovery rate (FDR) less than 0.01. In the LO-PCa group, 43 miRNAs were identified as differentially expressed, including 26 upregulated and 17 downregulated miRNAs. In the EO-PCa group, 55 miRNAs were identified as differentially expressed, including 28 upregulated and 27 downregulated miRNAs. Subsequently, the two lists of differentially expressed miRNAs included 69 miRNAs DE in total with 29 miRNAs in common, and further analyses were conducted ( Fig. 1  Functional enrichment analysis of miRNAs with dysregulated expression. KEGG pathway enrichment analysis was successively predicted by miRNet and aimed to validate that these miRNAs DE are involved in the prostate cancer pathway. This analysis revealed that in EO-PCa samples, 23 and 44 pathways www.nature.com/scientificreports/ were present in upregulated and downregulated miRNAs, respectively. In LO-PCa, 10 and 33 pathways were upregulated and 33 downregulated, respectively (Fig. 2). Genes in the PCa pathway as annotated by KEGG were overrepresented as targets in upregulated and downregulated miRNAs in EO-PCa and LO-PCa. The detected pathways are dysregulated by upregulated and downregulated miRNAs in EO-PCa and by downregulated miR-NAs in LO-PCa. These pathways had roles in carcinogenesis, such as increasing the cellular proliferation rate, reducing cellular focal adhesion, and alteration of signaling pathways, such as MAPK, p53, Jak-STAT, neurotrophin, Wnt and ErbB. These results are similar to previous reports in renal cell carcinoma and thyroid cancer. The identified pathways dysregulated by upregulated miRNAs in LO-PCa included cellular proliferation, p53 signaling pathway, protein processing in the endoplasmic reticulum, adherens junction formation, and amino acid lysine degradation. Supplementary File 1 shows the target genes present in each of the dysregulated pathways.

Weighted coexpression networks.
To capture most biological differences in the compared groups, we used less stringent criteria to select differentially expressed miRNAs: FC greater than 1.5 or less than -1.5 and FDR less than 0.05. With those parameters, we identified in total 157 miRNAs differentially expressed (DE), 102 miRNAs DE in EO-PCa, 121 miRNAs DE, and 66 miRNAs DE in common. Based on the expression levels of the 157 miRNAs DE, two coexpression networks were calculated. To generate comparable Weighted Coexpression Networks (WGCA), we used the same list of miRNAs in each of the two networks generated. First, the 24 samples from young patients were used, and the similarity threshold was calculated using the maximum local method. A Pearson correlation coefficient greater than or equal to 0.57 differentiates the distributions of the correlations between the miRNAs and those of a random population. In total, 35 miRNAs had a correlation coefficient greater than the threshold. A similar process was performed with the 25 samples from the oldest patients using the same list of miRNAs. In this case, a Pearson correlation coefficient greater than or equal to 0.62 differentiates the distributions of the correlations between the miRNAs and those of a random population. In LO-PCa, 33 miRNAs were coexpressed over the threshold and were used for the following analysis.
miRNA-miRNA interactions to detect hub miRNAs. The networks calculated were analyzed using Cytoscape. The top 20 nodes ranked by the metrics MCC implemented CytoHubba were chosen 22 . In young and elderly samples, two well-defined networks composed of downregulated and upregulated miRNAs were observed ( Fig. 3 and Table 1). In EO-PCa, the network had 35 nodes, one connected component, a clustering coefficient of 0.695, network centralization of 0.250, and an average of 6.97 several partners of neighbors. On the other hand, in LO-PCa, the network had 33 nodes, two connected components, a clustering coefficient of 0.740, network centralization of 0.180, and an average of 5.394 several partners of neighbors. The score of the top 10 hub miRNAs was from 1560 to 7716 in the young network and from 60 to 182 in the old network. In the EO-PCa top 10 miRNAs, all miRNAs were downregulated. However, in the LO-PCa network, the top 6 miRNAs were upregulated, and the following 4 were downregulated. In the comparison between the top 10 and the top group, hsa-miR-3065-3p and hsa-miR-676-3p were exclusively identified in young patients, whereas these genes did not exhibit statistical significance in the elderly population. Additionally, in the top 20, hsa-miR-488-3p was www.nature.com/scientificreports/ exclusively identified in young patients, but this gene was not statistically significant in elderly patients. In the top 20 miRNAs, nine downregulated miRNAs were present in the two groups (hsa-miR-31-5p, hsa-miR-224-5p, hsa-miR-205-5p, hsa-miR-205-3p, hsa-miR-3545-3p, hsa-miR-224-3p, hsa-miR-1911-5p, hsa-miR-1912, and hsa-miR-452-5p). In addition, only three miRNAs were upregulated in young patients (hsa-miR-150-5p, hsa-miR-142-5p, and hsa-miR-146b-3p). This result is in contrast to elderly patients, where eight miRNAs were upregulated (hsa-miR-32-5p, hsa-miR-96-5p, hsa-miR-182-3p, hsa-miR-183-5p, hsa-miR-375, hsa-miR-183-3p, hsa-miR-32-3p, and hsa-miR-5096). The exclusively upregulated and downregulated hub miRNAs in the young patients were not statistically significant in the elderly patients.
To identify the relevance of the hub miRNAs in the EO-PCa network we used miRNet to perform a KEGG pathway enrichment analysis with the predicted transcriptional targets of the 20 hub miRNAs. The analysis Assessment of prognostic significance of EO-PCa cohort. The prognostic significance of the 20 hub miRNAs of the EO-PCa network were analyzed via PROGmir V2, which used the PRAD dataset 23 . Of the 20 hub miRNAs, four had prognostic value using survival data from general PCa. The hub miRNAs that are exclusively dysregulated in EO-PCa were hsa-miR-3065 (Hazard ratio (HR): 1.3) and hsa-miR-146b (HR: 1.34), which are associated with poor relapse-free survival. In addition, miR-676 (HR: 1.84) was related to poor overall survival. On the other hand, two miRNAs exclusively identified in old patients, namely hsa-miR-32 and hsa-miR-96 (HR: 1.46), were associated with poor relapse-free survival. In addition, hsa-miR-10a (HR: 2.31) was related to poor overall survival. Figure 4 presents the Kaplan-Meir survival plots and the number of events in each analysis.

Correlated mRNA genes with hub miRNAs from the EO-PCa coexpression network.
To molecularly explain the effects of dysregulation of those miRNAs in EO-PCa, we identified genes that were coex- www.nature.com/scientificreports/ pressed with the hub miRNAs. We used mRNA microarray data generated using the same samples from the same patients who were used to generate the miRNA expression data. After normalization and batch effect correction of the microarray dataset, six samples were identified as outliers, three in EO-PCa patients and three in LO-PCa patients. Those samples were not included in the additional analysis.
The Pearson coefficients of the correlation of the 20 hub miRNAs in EO-PCa with the genes in the microarray were calculated for the EO-PCa samples and similarly for the LO-PCa samples. A permutation test revealed that correlation coefficients greater than or less than ± 0.614 were statistically significant for EO-PCa samples and correlation coefficients greater than or less than ± 0.673 were statistically significant for LO-PCa samples.
Different numbers of coexpressed genes were identified for each hub miRNAs with numbers ranging from 0 to 167 genes. Supplementary file 2 presents the numbers, names and correlation coefficients of the coexpressed genes for EO-PCa and LO-PCa samples. The most remarkable result was that EO-PCa upregulated hub miRNAs had more coexpressed genes than downregulated miRNAs, and most of these miRNAs exhibited positive correlation coefficients (Supplementary file 2).
To provide a biological meaning of the lists of coexpressed genes, a pathway analysis with the statistically significant coexpressed genes for each hub miRNA was performed. We performed an overrepresentation analyses against the 21 cancer prostate pathways that we collected from Molecular Signatures Database (MSigDB) 24  . We found that the three upregulated miRNAs (hsa_miR_142_5p, hsa_miR_146b_3p and hsa_miR_146b_3p) were coexpressed with the DEGs in EO-PCa versus normal tissue (Ding Early onset prostate cancer 2016) 9 . Two of the upregulated miRNAs were coexpressed with DEGs upregulated in the more aggressive prostate cancers of African-Americans compared with the less aggressive prostate cancers of European-American patients (WALLACE PROSTATE CANCER RACE UP) 25 ( Table 2). Among the downregulated miRNAs, hsa_miR_3545_3p and hsa-miR-224-5p were significantly coexpressed with genes in LIU PROSTATE CANCER DN 26 (Table 2). This pathway was obtained from microarray analysis of 31 PCa samples. The Gleason score was variable, and information on patient age was not provided. Table 3 shows all the genes with significant correlation coefficients for Ding Early onset prostate cancer 2016 and WALLACE PROSTATE CANCER RACE UP.
A similar analysis was performed using the LO-PCa data, demonstrating that none of the miRNAs had more coexpressed genes in the pathways of early onset prostate cancer or more aggressive cancer in African-Americans than expected by chance. Four miRNAs (hsa_miR_31_5p, hsa_miR_205_5p, hsa_miR_224_3p, and hsa_miR_3545_3p) were significantly coexpressed with genes in the pathway LIU PROSTATE CANCER DN 26 (Supplementary Table 2). Assessment of prognostic significance of genes coexpressed with hub miRNAs. Using GEPIA 27 , which employed the PRAD dataset, the genes correlated with hub miRNAs that exhibited statistical signifi- www.nature.com/scientificreports/ Table 2. Over-representation analysis of co-expressed genes with hub miRNA using EO-PCa data. Only statistically significant associations are shown. Size means the number of genes in the pathway. Expect means the expected number of genes in the pathway co-expressed by chance with the miRNA. Ratio means the additional number of times that there are more genes co-expressed with the miRNA compared with the expected number. Overlap means the number of co-expressed genes that are part of the pathway. FDR is the false discovery rate of the overlap. www.nature.com/scientificreports/ cance were analyzed (Table 3)

Discussion
EO-PCa is a subtype of PCa, which is currently receiving high interest due to its impact on clinical behavior and pathobiological differences with the "classical" or elderly PCa (LO-PCa). In this study, novel data analysis was performed using transcriptomic data from patients with EO-PCa who were < 45 years old and LO-PCa who were 71 to 74 years old. The data analyzed were generated by Ding et al. using 49 PCa patients 9 . The tumor samples were GS 7 (3 + 4) and grade T2 (T2a or T2c). Samples were obtained from different ethnic groups, including 88% whites, 4% African-Americans, 4% Hispanics and 4% Asians. We identified 55 miRNAs DE in EO-PCa, including 28 upregulated and 27 downregulated. In addition, 26 of these genes were exclusively dysregulated in EO-PCa.
Using an overrepresentation analysis with the predicted targets genes of the miRNAs DE, we identified several pathways commonly dysregulated between EO-PCa and LO-PCa. These pathways are related to adherence junctions, cell cycle and p53 signaling. In addition, the neurotrophin signaling pathway was identified as dysregulated, and members of this pathway are expressed in PCa, i.e., trk receptors and neurotrophins (NGF, BDNF, and/or NT-3). 28. Strikingly, in upregulated miRNAs in elderly patients, the lysine degradation pathway exhibited statistical significance. Lysine modification is associated with carcinogenesis in different types of tumors 29 . Among the hub miRNAs exclusively deregulated in EO-PCa without significance in LO-PCa, three were downregulated (hsa-miR-3065-3p, hsa-miR-676-3p, and hsa-miR-488-3p), and three were upregulated (hsa-miR-150-5p, hsa-miR-142-5p, and hsa-miR-146b-3p). The role of hsa-miR-3065-3p and hsa-miR-676-3p in PCa is unknown, but reduced expression of these genes is observed in tumors, such as esophageal squamous cell carcinoma 30 , hsa-3065-3p in clear cell renal cell carcinoma 31 , and hsa-miR-676-3p in breast cancer cell lines 32 . In survival analysis, hsa-miR-3065 was associated with poor relapse-free survival, and hsa-miR-676 was related to poor overall survival. In the PCa cell lines (LNCaP and C4-2B), hsa-mir-488-3p inhibits the androgen receptor, blocks proliferation and induces apoptosis 33 . Among the exclusively upregulated hub miRNAs in EO-PCa, hsa-miR-142-5p and hsa-miR-150-5p had been described as deregulated [34][35][36][37] . The last may act as antitumor miRNA targeting SPOCK1 35 . hsa-miR-150 has been reported as upregulated with a role in proliferation and invasion by targeting p53 36 , and its expression is associated with poor overall survival (HR: 1.87, CI: 1.19-2.94) 37 . In addition, Ding et al. 9 found that hsa-miR-146b-3p exhibited the highest level of overexpression in young PCa patients. However, the specific role of hsa-miR-146b-3p in PCa its unknown, it was associated with poor relapse-free survival in the survival analysis (Fig. 4). hsa-miR-146b-3p is member of the miR-146 a/b family. hsa-miR-146a www.nature.com/scientificreports/ is dysregulated in PCa and other tumors. In androgen-independent PCa, its downregulation is involved in apoptosis through regulation of the ROCK/Caspase 3 pathway 38,39 .
To understand how dysregulation of hub miRNAs modulates the normal behavior of prostatic tissue, we searched the targets genes of each hub miRNAs. The number of potential target genes identified was 112 for hsa-miR-3065-3p, 51 for hsa-miR-676-3p, 131 for hsa-miR-488-3p, 534 for hsa-miR-150-5p, 226 for hsa-miR-142-5p, and 56 for hsa-miR-146b-3p. These hub miRNAs are annotated as members in different KEGG pathways that are relevant in PCa biology. For example, regulatory targets of hsa-miR-3065-3p, such as GSK3B (Glycogen Synthase Kinase-3), are involved in apoptosis, cell cycle, DNA repair, tumor growth, invasion, and metastasis pathways. In recent years, it has become a targeted gene for therapy 40 . In the progression to androgen-independent PCa, GSK3B may act with PTEN 41 as a positive regulator of androgen receptor transactivation and growth independent of the Wnt/β-catenin pathway 42 . hsa-miR-3065-3p and hsa-miR-488-3p are related to p53, which is a very relevant gene that is mutated in greater than half of all cancers and is associated with progression in PCa 43 . hsa-miR-3065-3p is repress by mutant p53 44 , and hsa-miR-488-3p activates the p53 pathway through suppressing ZBTB2 45 . On the other hand, C-terminal Binding Protein 1 (CTBP1) is a target gene of hsa-miR-676-3p. This gene is a transcriptional corepressor of tumor suppressors genes involved in cell death, and dysregulated expression of this gene is associated with PCa progression 46 . Platelet-derived growth factor receptors-β (PGDFR-β) is regulated by hsa-miR-488-3p and hsa-miR-146b-3p), which are key regulators of cell growth and division 47 . In PCa, PGDFR-β is expressed in the early stage of the disease 48 . Its activation is associated with the loss of PTEN 49 , and high PGDFR-β expression is associated with prostate cancer recurrence 50 .
Among the upregulated hub miRNAs, several target genes involved in the KEGG PCa pathway were identified. The targets of hsa-miR-150-5p in this pathway include CDK2, EP300, and TP53. CDK2 is a key regulatory protein involved in cell cycle arrest upon DNA damage 51 . It is upregulation is associated with PCa progression, and it is a probable novel target gene in treatment 52 . hsa-miR-142-5p is modulator of important genes involved in the pathogenesis of PCa, including Cyclin D1 (CCND1), MAPK1, and PTEN. CCND1 is associated with aggressiveness 53 . MAPKs are serine/threonine kinases that mediate intracellular signaling associated with a variety of cellular activities, such as cell proliferation. In PCa, MAPKs are involved in apoptosis, survival, metastatic potential, and androgen-independent growth 54 . Additionally, PTEN is the most commonly lost tumor suppressor gene in primary disease. In most cases with PTEN loss, the gene is lost by genomic deletion. The loss of PTEN is associated with prostate tumor aggressiveness, progression, and poor prognosis (reduced disease-specific survival) 55 .
Given the limited number of studies on the effects of the hub miRNAs in the pathogenesis of PCa, we performed a correlation analysis based on the expression of hub miRNAs and all the genes in the genome in patients with EO-PCa. The overrepresentation analysis of the genes with significant correlations among hub miRNAs in the EO-PCa network revealed that the three upregulated hub miRNAs were significantly coexpressed with the Ding Early-onset prostate cancer 2016 pathway: two of the upregulated hub miRNAs were coexpressed with Wallace_Prostate_Cancer_Race_Up and one downregulated hub miRNA was coexpressed with LIU PROSTATE CANCER DN (Table 2). In the primary analysis of the data used in the present study, Ding et al. 9 reported differential expression of genes annotated in immunological pathways in the age:tissue interaction analysis (B Cell Development, iCOS-iCOSL Signaling in T Helper cells, CD28 Signaling in T Helper Cells, Primary Immunodeficiency Signaling, Calcium-induced T Lymphocyte Apoptosis), including genes such as complement family genes, immune-cell surface antigens, chemokines, interleukin receptors, natural killer cells and extracellular matrix remodeling genes. Moreover, the Wallace_Prostate_Cancer_Race_Up 25 dataset was generated from the comparison of gene expression profiles of PCa from 33 African-American patients with 36 European-American patients. The genes in this pathway are involved in immune response, stress response, cytokine signaling, and chemotaxis pathways. Several known metastasis-promoting genes, including autocrine mobility factor receptor, CXCR4 (chemokine (C-X-C motif) receptor 4), and MMP9, were more highly expressed in tumors from African-Americans than European-Americans. The expression profiles of two upregulated hub miRNAs, namely hsa-miR-150-5p and hsa-miR-142-5p, were correlated with the Wallace_Prostate_Cancer_Race_Up dataset.
The genes shared in Wallace and Ding pathways and statistically correlated with upregulated hub miRNAs (hsa-miR-150-5p and hsa-miR-142-5p) are genes involved in the immunology response, such as CCR7, IL7R, ITK, PTPRC, MMP9, APOE, CCL19, and CD3D. These genes were upregulated in PCa of African-American patients and in EO-PCa 25 . CCR7 and MMP9 are genes associated with PCa progression and metastases [56][57][58] . MMP-9 is involved in several hallmarks of PCa progression, such as proliferation, angiogenesis, epithelial to mesenchymal transition, apoptosis, and metastasis 59 . MMP-9 expression is associated with the risk of PCa (OR = 7.91; 95% CI: 5.27-11.89; P < 0.00001) 60 . A primary PCa cell line derived from an African-American patient (E006AA) exhibited increased MMP9 expression compared to other studied cell lines (LNCaP, C4-2, and MDAPCa2b) 61. CCR7 is a chemokine receptor that is associated with lymph node metastasis in other tumors, such as breast cancer 62 , non-small cell lung cancer 63 , and gastric carcinoma 64 . The CCR7 ligand CCL21 is expressed selectively in high endothelial venules at the entry point into the lymph node and promotes cancer progression 56,62. Polymorphisms in CCR7 (rs3136685) are present in African-American PCa patients 65 . IL7R and IL7 are highly expressed in PCa and are associated cancer cell invasion and migration probably by activating the AKT/NF-κB pathway and upregulating MMP-3 and MMP-7 expression 66 .
All the coefficients of correlation among the upregulated hub miRNAs and the mRNAs genes were positive. It could be counterintuitive that the upregulation of a specific miRNA causes the up-regulation of a specific mRNA. We proposed that it may explained by the activation of the transcription of target genes by binding of miRNAs with reverse complementary sequences in promoter regions of genes 14 . For example, overexpression of hsa-miR-205 increases IL-35 and IL-24 expression 67 . In PCa cells (DU145 and PC3), overexpression of miR-3619-5p induces CDK1N1 gene expression via direct interaction with the promoter region 68 . On the other hand, one gene target can be targeted by several miRNAs 69 . Therefore, the targets of this gene are potentially regulated by miRNAs DE. Thus, it is necessary study the regulation of this gene's targets. www.nature.com/scientificreports/ In the analysis of overall survival of genes with high expression in Wallace_Prostate_Cancer_Race_Up 25 , high expression of DOCK10 and ITGB2 was associated with poor prognosis. Both genes are coexpressed with hsa-miR-150-5p and hsa-miR-142-5p. The roles of these genes in PCa are unknown. In cancer, DOCK10 (Dedicator Of Cytokinesis 10) is involved in the regulation of the epithelial to mesenchymal transition 70 . ITGB2 (Integrin beta-2 (CD18)) combines with integrin alpha to form the integrin lymphocyte function-associated antigen-1 (LFA-1). This gene is involved in tumor growth and metastasis 71 .
A limitation of this study was the high homogeneity of the patients, who all had Gleason scores of 7 (3 + 4). Patients with low-grade PCa have different molecular characteristics, clinical behavior and treatment 72 . However, the expression profile of miRNAs in low-grade PCa does not exhibit significant differences. Walter B et al. 73 compared miRNA expression in 26 patients with low-grade and 15 patients high-grade PCa. The results did not reveal any specific signatures to differentiate the two groups. Another limitation involves the survival analyses of DE miRNAs and DEGs because the data were obtained from TCGA database. In this analysis, we could not perform separate evaluations for EO-PCa and LO-PCa.
In conclusion, this is the first study that analyzed the expression of miRNAs in EO-PCa and LO-PCa patients using network analysis. Connections among miRNA expression, target genes, and molecular pathways for EO-PCa and LO-PCa were identified. Furthermore, specific miRNAs with clinical significance in young patients may explain molecular differences, and the different biological processes in young and elderly patients were identified. In addition, we found coexpression of genes and hub miRNAs that play important roles in PCa progression and metastasis and genes associated identified in Afro-American PCa patients. Most of these genes are involved in the immunology response. As a recommendation, the constructed network of biomarkers should be further assessed in EO-PCa, and candidate miRNAs and gene targets should be validated for patient diagnosis and prognosis.

Materials and methods
Data selection. In June 2019, an advance search was performed to identify studies that analyze miRNA expression in EO-PCa patients. The sources used included PubMed and the National Center for Biotechnology Information (NCBI) GEO database (https ://www.ncbi.nlm.nih.gov/geo/). The keywords 'young OR early-onset AND prostate cancer' were used in the search. The results were limited to Homo sapiens as the organism and expression profiles were determined using array dataset types.
The inclusion criteria for the systematic review were (1) miRNA expression was assessed in prostate tissue of young and elderly PCa patients in the same dataset, (2) raw data were available, and (3) data passed quality control. Two reviewers performed an eligibility assessment by screening titles and abstracts from the publications. Subsequently, the articles that did not meet the eligibility criteria were rejected. Additionally, we searched The Cancer Genome Atlas (TCGA); however, we could not use these data because the database only contained matched tumor and normal data from three patients diagnosed with PCa under the age of 55 years.
Raw count miRNA data were normalized using the trimmed mean of the M-value (TMM) method through the EdgeR package 74 . Because the original data were sequenced in two batches 9 , this nonbiological variability source was corrected using nonparametric Bayesian statistics methodology using the sva package. 75 . The extent of the batch effect correction was assessed by principal component analysis.
The samples were assigned to two different experimental groups to identify the differentially expressed miR-NAs in LO-PCa (71-74 years old) and EO-PCa (38-45 years old). miRNA expression levels were compared between tumor and normal samples in each group. The statistical significance level for this study was calculated using the Limma statistics package (Linear Models for Microarray and RNA-Seq Data) 76 . Limma uses a linear modelling to detect differentially expressed genes. The fold change was calculated, and the statistical significance was adjusted for multiple comparisons (False Discovery Rate (FDR)).

Weighted gene coexpression networks analysis (WGCNA).
The total levels of differentially expressed miRNAs between normal and tumor samples in older and young patients were collected in one list. Two coexpression networks were developed using this gene list. The first list was generated for the young samples, and the second list was generated for the older group. First, the similarity matrix was calculated by identifying the Pearson correlation coefficients of the expression levels for the samples based on all possible gene pairs. Then, the similarity threshold was calculated with the adjacency function, which was established according to the unique characteristics of each similarity matrix 77 . The method developed by Elo was used to select the threshold 78 . This method compared the tau values for the network grouping coefficient (Co) with that expected for a random network (Cr). It uses the clustering coefficient of the real graph in comparison to a rando graph. The threshold for significant similarities is chosen so that the obtained real graph is scale free. Finally, the adjacency matrix (2 × 2) of the network was established and allowed the representation of binary relationships. In this case, a pair of genes that exhibit coordinated gene expression activity (coexpression) is indicated by (1); otherwise, a (0) is reported. All WGCNA analyses were performed in an R unique environment using statistical functions (https ://www.r-proje ct.org/). www.nature.com/scientificreports/ Detection of hub miRNAs. The hub miRNAs were identified through network analysis using Cytoscape and its plugin (CytoHubba). This plugin accurately identifies hub genes by 12 topological analysis methods. For this study, the Maximal Clique Centrality (MCC) method proposed by CytoHubba was used; recently, this method exhibits improved performance to capture essential targets in the top rank list in both high-and lowgrade PCa 22 . In addition, MCC helped to identify the top 20 hub miRNAs. On the other hand, the network analyzer plugin was used to recognize the network parameters.
Functional annotations of hub miRNAs. The miRNet database to identify target genes (https ://www. mirne t.ca/) was used to facilitate the interpretation of biological mechanisms related to hub miRNAs. This tool integrates data from eleven different miRNA databases: TarBase, miRTar-Base, miRecords, miRanda, miR2Disease, HMDD, PhenomiR, SM2miR, PharmacomiR, EpimiR, and starBase 79 . The following information was provided for miRNet analysis: organism name (H. sapiens), ID type (miRBASE ID), and tissue origin (Not specified). No degree or betweenness filter was used for network visualization. Additionally, miRNet was applied to identify biological pathways and processes, molecular functions, and cellular components that are statistically enriched for the corresponding miRNA target genes. For the functional evaluation of the miRNAs, the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted in mirNet. Only statistically significant annotation categories (P value < 0.05) were retained.
Survival analyses of miRNAs with dysregulated expression. PROGmiR V2 is an online free tool and is available at https ://www.compb io.iupui .edu/progm ir. This program combines the prognostic data of miR-NAs for different types of cancers from TCGA dataset. This tool was used to compare the overall, relapse-free, and metastasis-free survival of prostate adenocarcinoma patients with DE of miRNAs in young and old cohorts. It also divides samples based on high and low expression and calculates the hazard ratio (HR) with relative confidence intervals (CI) and P values for the proportional hazards model 23 .

Determination of mRNA expression levels of all genes in normal and tumor tissues from EO-PCa and LO-PCa for the diagnosis of prostate cancer.
The level expression of mRNAs from cancer and normal tissues from EO-PCa and LO-PCa patients were downloaded from the same dataset GSE88808 available in the GEO OMNIBUS database (https ://www.ncbi.nlm.nih.gov/geo/query /acc.cgi?acc=GSE88 808). These data were generated in parallel with the miRNA dataset. The same patients and same tissues were used to identify the level expression of mRNAs and miRNAs. A detailed description of the RNA obtention and determination of mRNA expression level are provided in the primary paper 9 . In summary, total RNA was extracted from approximately 5 mg of unsectioned formalin-fixed paraffin-embedded core samples using the RecoverAll Total Nucleic Acid Isolation kit (Life Technology, Inc.). The Illumina HumanHT-12 WG-DASL V4.0 expression beadchip was used for mRNA expression profiling of 29,000 genes in the human genome. The levels of intensity were normalized using the quantile normalization method 80 . The batch effects secondary to different times of hybridization were corrected using the empirical Bayes methods as is implemented in the ComBat in sva package 81 .
Identification of outliers was performed using the Pearson correlation measurements of the level expressions of all genes in the microarray between all the samples. Samples with correlation coefficients less than 0.9 compared with the other samples were excluded from additional analyses.
Analysis of the correlation of the expression of hub miRNAs with the expression of genes involved in the pathogenesis of prostate cancer. The hub miRNAs in the EO-PCa coexpression network were included in additional correlation analyses to determine whether they were coexpressed with genes involved in prostate cancer.
First, we collected the 183 differentially expressed genes (DEGs) identify by Ding et al.in the primary analysis of the mRNA expression data 9 ; they identified differences in tumor vs. normal tissues between samples from young and old patients. We refer to this list of genes as Ding Early-onset prostate cancer 2016.
Second, the Molecular Signatures Database (MSigDB) 24 was interrogated to collect the pathways related to prostate cancer. In total, 22 different pathways were identified (Supplementary file 3). This collection of pathways represents the state of knowledge about transcriptomic modifications between tumor tissues compared with normal tissues from prostate cancer patients (19 pathways) and tumor samples from African-American compared with European-American patients with primary prostate cancer 25 .
Third, Pearson correlation coefficients were calculated among the expression profiles of selected hub miRNAs and mRNA expression levels from all the genes in the Illumina microarray. We initially performed correlation analysis only for young or old samples. Using a permutation test, the confidence intervals were calculated. Correlations coefficients with P values less than 0.002 were selected as statistically significant.
Finally, overrepresentation analyses were performed using the hypergeometric test as implemented in WebGestalt 82 . For each selected hub miRNA, statistically significant coexpressed genes were interrogated against the genes in the 22 prostate cancer pathways to identify whether more (overrepresentation) genes coexpressed with miRNAs are present than expected by chance.

Survival analyses of genes coexpressed with hub miRNAs. Gene Expression Profiling Interactive
Analysis (GEPIA; https ://www.gepia .cance r-pku.cn) 27 was used to calculated disease-free survival and overall survival between DEGs coexpressed with hub miRNAs. The lower and upper 50% of gene expression levels were