Exploration of Extracellular Vesicle miRNAs, Targeted mRNAs and Pathways in Prostate Cancer: Relation to Disease Status and Progression

Simple Summary Prostate cancer lacks non-invasive specific biomarkers for aggressive disease. Urinary extracellular vesicles (uEV) could provide such markers; however, due to technical challenges, little is known regarding the pathogenesis pathways reflected in uEV. We performed a miRNA, target mRNA and pathway study focused on uEV, exploring the differences between cancer (1) status groups (Gleason score) and (2) progression groups. The uEV provided a surprisingly comprehensive presentation of differentially expressed miRNAs, target mRNAs and pathogenesis pathways. The miRNAs associated with prostate cancer status or progression were mostly unique, but still targeted overlapping sets of signalling, resistance, hormonal and immune pathways. Interestingly, mRNA targets of the key miRNAs (miR-892a, miR-223-3p, miR-146a-5p) were widely expressed in both uEV and plasma EV from PCa patients. The study thus suggests that uEV carry a vast presentation of PCa status and progression-linked RNAs that are worth further exploration in large personalized medicine trials. Abstract Background: Prostate cancer (PCa) lacks non-invasive specific biomarkers for aggressive disease. We studied the potential of urinary extracellular vesicles (uEV) as a liquid PCa biopsy by focusing on the micro RNA (miRNA) cargo, target messenger RNA (mRNA) and pathway analysis. Methods: We subjected uEV samples from 31 PCa patients (pre-prostatectomy) to miRNA sequencing and matched uEV and plasma EV (pEV) from three PCa patients to mRNA sequencing. EV quality control was performed by electron microscopy, Western blotting and particle and RNA analysis. We compared miRNA expression based on PCa status (Gleason Score) and progression (post-prostatectomy follow-up) and confirmed selected miRNAs by quantitative PCR. Expression of target mRNAs was mapped in matched EV. Results: Quality control showed typical small uEV, pEV, RNA and EV-protein marker enriched samples. Comparisons between PCa groups revealed mostly unique differentially expressed miRNAs. However, they targeted comprehensive and largely overlapping sets of cancer and progression-associated signalling, resistance, hormonal and immune pathways. Quantitative PCR confirmed changes in miR-892a (Gleason Score 7 vs. ≥8), miR-223-3p (progression vs. no progression) and miR-146a-5p (both comparisons). Their target mRNAs were expressed widely in PCa EV. Conclusions: PCa status and progression-linked RNAs in uEV are worth exploration in large personalized medicine trials.


Introduction
Prostate cancer (PCa), the most common malignancy of men in Finland and other Western countries, is one of the leading causes of cancer-related deaths. For PCa diagnosis, tumour tissue biopsies are the gold standard. However, biopsies cause problems due to their invasiveness, low yields of tissue, e.g., from bone metastasis, or insufficient representation across the heterogeneous tumour foci [1]. In contrast, the widely applied PCa biomarker, prostate-specific antigen (PSA), can easily be measured from serum for clinical diagnostics, prognostics and disease or therapy monitoring. The specificity of PSA is unfortunately poor, leading to over-diagnosis or -treatment [2,3]. Thus, in need of better PCa biomarkers, the PCa research field has turned its focus to liquid biopsies.
Urinary extracellular vesicles (uEV) are membrane-enclosed vesicles derived mainly from the cells of the genitourinary system. They carry internal and external molecular cargo such as nucleic acids, proteins, lipids and metabolites [4,5]. UEV are of interest, e.g., due to their potential functions in intercellular communication and cargo, which is considered as a non-invasive source of biomarkers for different pathological conditions, including PCa [6,7]. Particularly, the microRNA (miRNA) content of uEV from patients with PCa has raised considerable interest [8]. Most studies have evaluated the potential of uEV miRNAs in PCa detection (diagnosis), and included uEV technology development [9][10][11][12]. However, less is known about how PCa progression or therapy resistance could change the contents of the uEV in clinical samples. Few studies indicate changes in urine or uEV miRNAs associated with biochemical recurrence [13,14] or between low vs. high grade and localized vs. metastatic cancers [15][16][17]. A recent study elucidated the effect of androgen manipulation (dihydrotestosterone and enzalutamide) on the EV contents from LnCaP cells, including miRNAs [18]. However, the discovered miRNA biomarker candidates have differed between studies, rendering it difficult to evaluate uEV's potential to provide information regarding the PCa pathogenesis pathways. Equally little is known regarding whether analysis of miRNA with their target mRNA from EV adds valuable insights into PCa pathogenesis, or how EV from urine and plasma differ in their RNA cargo. These are interesting questions because, along with their own challenges, both sample types also have unique potential in PCa research and diagnostics, e.g., urine as a source of EV from primary tumour and blood as a source of EV from distant metastasis [19].
Technically, RNA research into uEV is challenging. The best practices in EV pipelines are currently still developing, and they include pre-analytics (sample collection, processing and storage), EV isolation, RNA detection and data normalization/analysis [5,7,20,21]. The International Society of Extracellular Vesicles (ISEV) has launched several standardization efforts to expedite the development, including, e.g., Minimal Information for Studies of Extracellular Vesicles (MISEV) guidelines, and targeted ISEV position papers [5,22]. However, so far, the heterogeneity of the uEV pipelines, on top of the heterogeneity of the PCa tumours and study designs, likely explains why the discovered miRNA biomarker candidates differ [8]. For example, heterogeneity of the isolated uEV and copurified non-EV materials from different isolation workflows impact the miRNA results [7,23,24]. As the average sizes, biogenesis routes and cargo of small and large EV differ [22,[25][26][27], targeting small EV for biomarker analysis is one strategy to limit heterogeneity of the EV and dissect exosomal messages. Prior studies have further shown that PCa cells secrete large quantities of exosome-sized EV, e.g., due to hypoxia and low pH, typical for the solid tumour microenvironment [28,29]. Finally, targeting small EV also allows better sample purification i.e., preclearing steps to remove cells, bacteria or other bigger contaminants before EV purification.
In this study, we focused on small uEV enriched samples to explore whether miRNAs and their target mRNAs could provide biomarker candidates and information regarding the different pathogenesis pathways linked to PCa status or progression.

Study Participants and Groups
Urine or EDTA plasma samples and/or clinical data were obtained from consented donors of the Helsinki Biobank, Helsinki Urological Biobank, and the SalWe "Get It Done" research project and from the FinnProstata IX study.
For the main study, PCa patients (n = 30, age range  were classified into three cancer status groups according to the following criteria: Group A had a Gleason score of ≥8; Group B had a Gleason score of 7 and were lymph node metastasis positive and/or had margin positive prostatectomy tissue sample and/or PSA > 0.2 ng/mL after radical prostatectomy (post-RP); Group C patients had a Gleason score of 7, were lymph node and margin negative and their PSA remained < 0.2 ng/mL post-RP. PSA measurements were taken ≥33 days post-RP. Healthy technical/biological controls; Group D (n = 10) were young (<45 years) asymptomatic men.
The outcomes of the patients were followed for 4-7 years post-RP. Then, PCa patients were divided into new groups according to disease progression. Group I contained individuals who died due to PCa and/or had metastasis and/or were hormonally treated, Group II obtained secondary treatment and/or had a biochemical recurrence and Group III had no progression during follow-up.
For the correlation study, the PCa patients (n = 3) were classified according to the above status criteria with the addition of group E having metastatic castration resistant PCa at the time of urine and plasma sample collection. The Gleason scores of Group E were determined through needle biopsies as the patients did not undergo RP.

Urine and Plasma Sample Collections and EV Isolation
For the main study, all urine samples were from spot mid-stream urine collections obtained before RP (pre-RP, 0-43 days before, urine collection dates available for 25 out of 30 PCa patients). Samples were kept cold (ice, cold package or +4 • C) during storage and delivery to the laboratory (average 4-6 h), where they were centrifuged at 1800× g, 10 min, +4 • C. Supernatants were frozen at −80 • C or liquid nitrogen vapor phase. Samples were isolated using a modification of protocol described producing small uEV enriched samples [4,7,30]. Briefly, urine samples were vortexed for 90 s and centrifuged at 8000× g at +4 • C for 15 min using a fixed angle AG-6512C rotor (Kubota Corp., Tokyo, Japan). Supernatants were filtered with Whatman 1.2 µm cellulose acetate syringe filters (GE Healthcare, Buckinghamshire, UK) and 18 mL of processed urine diluted to 30 mL with PBS, which was then centrifuged at 4 • C for 90 min at 100,000× g (maximum breaking) in a SW28 or SW32 rotor, k-factor 251 and 266, respectively (Beckman Coulter, Inc., Brea, CA, USA). Supernatants were discarded and the pellets suspended in filtered PBS (0.22 µm PES filter; Jet Bio-Filtration, Guangzhou, China) and stored in protein or DNA LowBind tubes (Eppendorf) at −80 • C. For uEV quality control with nanoparticle tracking analysis (NTA), Western blotting, electron microscopy (EM) or sequencing, 9-18 mL of urine available for six patient samples, or 18 mL of pooled healthy control samples (HcI, n = 11 men, and HcII, n = 9, 2 men and 7 women) were subjected to uEV isolation following the same protocol.
For the correlation study, matched urine and EDTA plasma samples were collected from all non-fasting donors on the same day. P33 gave samples on the day of RP (pre-RP) and 6 months post-RP, and P34 and P35 gave samples 2-3 weeks after and before needle biopsy, respectively. Urine was collected, processed and stored similarly as above and uEV isolated from 17-30 mL urine as described for small uEV enriched samples [4,7,30]. Blood was drawn in K2F EDTA BD Vacutainer ® tubes (BD Biosciences, Franklin Lakes, NJ, USA) and centrifuged at 2000× g for 10 min, and plasma aliquots were frozen at −80 • C or the liquid nitrogen vapor phase. ExoEasy Maxi kit (Qiagen, Hilden, Germany) was used for EV isolation from 0.9-1.5 mL of plasma, according to the manufacturer's instructions, including the pre-filtration (0.8 µm pore size, Millex AA, Merck Millipore, Burlington, MA, USA) of plasma. Eluted EV were frozen at −80 • C in DNA LowBind tubes (Eppendorf, Hamburg, Germany). All samples of HC11 were processed as three technical replicates. For uEV quality control with EM, Western blotting, or sequencing, plasma from a healthy Cancers 2022, 14, 532 4 of 29 man, HC12, pooled healthy donor urine samples (HcIII, n = 3 men) and matched plasma and urine samples from a healthy woman, HC13, were collected and processed similarly as correlation study samples.
The particle concentration and size distribution were analysed from uEV preparations by an NTA instrument LM14C equipped with a violet (405 nm, 70 mW) laser (Malvern Instruments Ltd., Malvern, UK) and an sCMOS camera (Hamamatsu Photonics K.K., Hamamatsu, Japan). Capture settings were temperature 18.2-19.9 • C, viscosity 1.035-1.043 cP and camera level 14. EV samples were diluted with filtered PBS, and five 30s videos were recorded with 36-102 particles/frame. Analysis was carried out with auto settings and detection threshold 4 in Nanosight software 3.1 (Malvern Instruments Ltd., Malvern, UK).

RNA Extraction and miRNA Sequencing in the Main Study
Isolated uEV samples were submitted to Qiagen Genomic Services (Hilden, Germany) for miRNA sequencing (miRNAseq). Total uEV RNA was extracted using Qiagen's miRNeasy serum/plasma kit according to the manufacturer's instructions. The library preparation was carried out using the QIAseq miRNA Library Kit with unique molecular identifiers (UMI) tagging (Qiagen). The cDNA was amplified using PCR (22 cycles). Library quality was monitored using either Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA) or TapeStation 4200 (Agilent). Sequencing was performed with NextSeq500 and 75 basepairs (bp) single-end reads (Illumina Inc., San Diego, CA, USA).
Raw data was de-multiplexed and FASTQ files for each sample were generated using the bcl2fastq software (Illumina Inc., San Diego, CA, USA). FASTQ data were checked using the FastQC tool. Sequence annotation was performed using Homo sapiens reference genome GRCh37 and miRBase 20 as an annotation reference. Bowtie2 (2.2.2) was used for mapping the reads. Quality was monitored through UMI correction, quality score (Q-score) of incorrect base call probability using 30 as a cut-off (an error probability < 0.001), read length (>15 nucleotides (nt) as cut-off) and mapping to reference genome or miRbase, as well as through read numbers (Table S1).
Differential expression analysis was performed using the CLC Genomic Workbench version 20.0.4 (Qiagen) and EdgeR statistical software package (Bioconductor, http:// bioconductor.org/, (accessed on 12 October 2020)). For normalisation, counts per million (CPM) or the trimmed mean of the M-values method (TMM) based on log-fold and absolute gene-wise changes in expression levels between samples were used [31]. Principal component analysis (PCA) was performed with CLC Genomic Workbench or R using normalized quantifications. The Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathways were analysed with miRWalk [32] using all nominal significantly differentially expressed (DE) miRNAs (FDR (false discovery rate) p or p < 0.05) as input and filtering for 3 UTR (untranslated region) interactions and for targets found using TargetScan, miRDB and miRTarBase. Tam 2.0 tool overrepresentation analysis [33,34] was also performed with all nominally significant DE miRNAs as inputs and the settings of all curated miRNAs as background, limited to a miRNA set size of two, including up-and down-regulated miRNA sets.

Quantitative PCR and Analysis
For technical confirmation of the miRNAseq results, the quantitative PCR (qPCR) of selected targets and reference candidate miRNAs were carried out at Qiagen (all tested assays listed in Table S2). The DE miRNAs were selected based on two or more of the following criteria: transcripts per million (TPM) ≥10 in all samples of a group, FDR p or p < 0.05, log2FC (fold change) >0.6 or <−0.6, as well as visual evaluation of variance and separation between groups. RNA (1.7 µL) was reverse transcribed using the miRCURY LNA RT Kit (Qiagen), and 50-fold dilutions of cDNA were assayed once for each miRNA on the miRCURY LNA miRNA PCR Custom panel using SYBR Green master mix. The amplification was performed in a LightCycler ® 480 Real-Time PCR System (Roche Diagnostics, Mannheim, Germany). The amplification curves were analysed using the Roche LC software. Quality control was performed at reverse transcription (RT) step by inclusion of RNA and DNA spike-ins, by melting curve and amplification efficiency analysis as well as by comparisons of Cq (quantification cycle) values to background level in the negative control samples (no template in the RT step). Inclusion criteria for analysis was that assays yielded signals ≥3 Cq lower than the negative control in at least 20% of samples, or <37 Cq in the case of there being no signal from negative control.

RNA Extraction and miRNA and mRNA Sequencing for the Corelation Study of Three Patients
Total RNA from uEV and pEV samples and 300 µL of plain plasma were isolated using a miRNEasy micro kit (Qiagen) according to the protocol for animal cells without DNAse treatment and with the exception that lysis was performed using Trizol LS (Life Technologies Corp., Carlsbad, CA, USA) due to its suitability for larger sample volumes. The quality and quantity of the RNA was monitored with a Bioanalyzer 2100 (Agilent) Pico kit.
The mRNA sequencing (mRNAseq) was performed as described [20]. Briefly, for sequencing libraries, 400 pg of total RNA were prepared with SMART-Seq v4 Ultra Low Input RNA Kit for Sequencing (Takara BIo Inc., Mountain View, CA, USA) and Nextera XT kit (Ilumina Inc., San Diego, CA, USA) according to the manufacturer's protocols. Libraries were sequenced on the Illumina Hiseq 2000 platform (Illumina Inc.) as 100 bp cycle paired-end reads. The sequences of adaptors were removed from the reads and the reads were aligned on human genome (version GRCh38) using gene annotation from ensemble database version 81 using STAR [41]. The miRNA sequencing libraries were prepared using 1 ng of total RNA for SMARTer ® smRNA-Seq Kit for Illumina ® (Takara Bio USA, Inc., San Jose, CA, USA), and 15 pM loading concentration was then used in sequencing with MiSeq (Illumina Inc.) as 75 bp paired-end reads. The miRNAseq data analysis was carried out essentially as above but without UMI correction. The mRNA targets of miR-146a-5p, -223-3p and -892a were downloaded from miRWalk [32], including all targets found in TargetScan, miRDB and Cancers 2022, 14, 532 6 of 29 miRTarBase (accession date 22 December 2021). Expression level differences were calculated using CPM data from both miRNA-and mRNAseq (Tables S3 and S4).

Statistical Testing and Venn Analysis
Statistical testing for miRNAseq read numbers was carried out with Student's t-test. The statistical significance of patient demographic and clinical data and DE miRNAs in miRNAseq, qPCR and pathway analyses were assessed using t-test or ANOVA with and without Benjamini-Hochberg FDR correction [42]. In overrepresentation analysis with Tam 2.0, significance was additionally tested with multiple correction via the Bonferroni method. In addition, in qPCR, the normality of the data was tested with the Shapiro-Wilk method [43]. In case of non-normally distributed data, a non-parametric unpaired Wilcoxon test [44] was used for assessing statistical significance; p-values < 0.05 were considered statistically significant. Venn analyses were carried out with a Venn diagram tool (https: //www.vandepeerlab.org/?q=tools/venn-diagrams, (accessed on 30 December 2021)).

Design of the Main Study
Our study aimed to clarify whether miRNA contents in uEV from PCa patients differed (1) between PCa status groups and (2) between PCa progression groups ( Figure 1). For these aims, we prepared uEV and carried out uEV quality control and miRNAseq from the urine samples of 30 PCa patients prior to their radical prostatectomy (RP) and from 10 healthy non age-matched men as additional technical and biological controls. The Gleason score and other histological characteristics of cancer in the prostatectomy tissue, or PSA levels post-RP, were used for dividing the patients into cancer status groups (Table 1). Group A, with the most aggressive disease, had a Gleason score of ≥8 (n = 10). Group B (n = 9) had a Gleason score of 7 and possible or verified metastasis or elevated postoperative PSA. Group C (n = 11) had a Gleason score of 7 and no indication of metastasis or elevated postoperative PSA. The healthy group D consisted of younger men who had the minimal possibility of asymptomatic PCa. After prostatectomy, we followed the disease progression in the patients for 4-7 years and re-grouped the patients to form progression groups (Table 1). Patients in Group I (n = 5) had the most severe disease progression with death, metastasis and/or hormonal treatment, Group II (n = 11) obtained secondary treatment and/or had a biochemical recurrence and Group III (n = 14) had non-progressive disease. For both status and progression groups, we analysed DE miRNAs and their regulated pathways and confirmed the selected miRNAs by qPCR using the same set of samples.

uEV and miRNA Sequencing Quality
We isolated uEV using a modified ultracentrifugation (UC) protocol due to the low sampling volumes of the collected urine samples (18 mL as compared to 30 mL protocol used in, e.g., [7]). Thus, we performed a thorough quality control of the uEV from representative individual patient and pooled control samples by EM, NTA, Western blotting and RNA profiling analysis (Figure 2, n = 3-8 individual samples per method). The EM images showed similar small uEV enriched samples as before [4,7,20,30] with variable levels of THP, but no obvious other contaminants such as cellular remnants or intracellular organelles, including mitochondria. By NTA, particle mean (166-189 nm) and mode (118-139) sizes were similar between samples, while concentrations varied almost 10-fold (1.6 × 10 9 -1.1 × 10 10 particles per ml urine), with both the minimum and maximum concentration detected in patient samples. Western blotting further confirmed the presence of uEV as the typical commonly enriched uEV markers; CD9, CD63, CD59, TSG-101 and PDX were detected in the samples, albeit in variable quantities. In addition, we detected the presence of THP and also some calnexin, an endoplasmic reticulum marker ( Figure S1). Extracted total RNAs were analysed by BioAnalyzer 6000 Pico RNA assay and showed a small RNA peak between 25-300 nt and no or a small amount of 18S and 28S rRNA. Study outline shows urine sample collection from prostate cancer (PCa) patients (time 1, t1) before radical prostatectomy (RP) and subsequent urinary extracellular vesicle (uEV) isolation with ultracentrifugation (UC) and miRNA sequencing (miRNAseq). The uEV quality control was carried out via RNA profiling, electron microscopy (EM), nanoparticle tracking analysis (NTA) and Western blotting (WB). After miRNAseq, differentially expressed (DE) miRNAs were compared between PCa status groups (A-C) and healthy technical controls (D). Groups A-C were based on Gleason scores (GS) 7-9, determined in prostatectomy tissue and other findings. Group B differed from C by more severe histological findings, metastasis or higher prostate specific antigen levels post-RP. Cancer progression was followed for 4-7 years post-RP and then (t2) patients were reclassified to groups (I-III) according to disease aggressiveness for reanalysis of the miRNAseq data. Groups I-III were based on death due to prostate cancer, metastasis or hormonal treatment (I), secondary treatment or biochemical recurrence (II) or no events (III) during followup time. Selected DE miRNAs from both analyses were confirmed by quantitative PCR (qPCR). Created with BioRender.com. Figure 1. Outline of the main study. Study outline shows urine sample collection from prostate cancer (PCa) patients (time 1, t1) before radical prostatectomy (RP) and subsequent urinary extracellular vesicle (uEV) isolation with ultracentrifugation (UC) and miRNA sequencing (miRNAseq). The uEV quality control was carried out via RNA profiling, electron microscopy (EM), nanoparticle tracking analysis (NTA) and Western blotting (WB). After miRNAseq, differentially expressed (DE) miRNAs were compared between PCa status groups (A-C) and healthy technical controls (D). Groups A-C were based on Gleason scores (GS) 7-9, determined in prostatectomy tissue and other findings. Group B differed from C by more severe histological findings, metastasis or higher prostate specific antigen levels post-RP. Cancer progression was followed for 4-7 years post-RP and then (t2) patients were reclassified to groups (I-III) according to disease aggressiveness for reanalysis of the miRNAseq data. Groups I-III were based on death due to prostate cancer, metastasis or hormonal treatment (I), secondary treatment or biochemical recurrence (II) or no events (III) during follow-up time. Selected DE miRNAs from both analyses were confirmed by quantitative PCR (qPCR). Created with BioRender.com. Table 1. Clinical characteristics of groups in the main study. Data of prostate cancer patients in the status Groups A-C or healthy technical controls, D, or the same patients classified into prostate cancer progression Groups I-III after 4-7 years follow-up. Numbers in parenthesis denote the subject number for whom the information was available-in the absence of parenthesis, the information was available for all patients. After prostatectomy (post-RP), prostate-specific antigen (PSA). centration detected in patient samples. Western blotting further confirmed the presence of uEV as the typical commonly enriched uEV markers; CD9, CD63, CD59, TSG-101 and PDX were detected in the samples, albeit in variable quantities. In addition, we detected the presence of THP and also some calnexin, an endoplasmic reticulum marker ( Figure  S1). Extracted total RNAs were analysed by BioAnalyzer 6000 Pico RNA assay and showed a small RNA peak between 25-300 nt and no or a small amount of 18S and 28S rRNA. The RNAs were subjected to miRNAseq, which gave good quality sequences (Q-score > 30) and resulted in 17.4M average raw reads (range, 13.7-24.1 M) and 3.2 M average UMI corrected reads (range, 1.5-6.4 M, Figure S2). There were no statistically significant differences between read numbers from healthy controls and patients. The sequencing read distribution showed that an average of 30% (range, 7-57%) aligned to miRNAs, 8% (1-25%) to other small RNAs and 40% (17-66%) were either out-or unmapped. Under 1% aligned to predicted or putative miRNAs. The uEV from PCa patients yielded 2.2-, 2.7-and 8.3-fold more reads of miRNA, small RNA and predicted miRBase miRNA, respectively, than uEV from controls (p < 0.05). On average, we identified 461 (range, 341-545) miRNAs that were expressed at the ≥1 TPM level in all 40 samples. Out of these, 247 (range, 181-247) were expressed robustly at the ≥10 TPM level. Despite the larger miRNA read numbers from PCa uEV samples, the number of identified miRNAs did not differ significantly between PCa (average 467/251 miRNAs with 1/10 TPM, respectively) and healthy controls (average 445/238 miRNAs with 1/10 TPM, respectively). The two pools (HCpI and HCpII) used for quality control of the uEV isolation workflow (Figure 2) also produced good sequencing quality, with over 29 M raw reads, 7 M UMI corrected reads, 450/250 identified miRNAs (expressed at 1/10 TPM level) and a similar read distribution to the study samples ( Figure S2D).

uEV from Prostate Cancer Patient Status Groups Differed in the Quantities of miRNAs Targeting Cancer and Progression-Linked Signaling, Resistance and Hormonal Pathways
We started miRNAseq data exploration of the PCa patient status groups (A-C) and healthy controls (D) by PCA focused on the 50 most variable miRNAs ( Figure S3). PCA did not show clustering according to PCa, age or health status groups, but revealed some non-clustering PCa samples including an outlier, patient 30 (P30), in status group A. However, no clear status-based groupings were evident, even when P30 was removed from the analysis ( Figure S3A,B). Hence, we decided to include this patient in the differential expression comparisons (Table S5), but also provide the data without P30 as additional supplementary information (Table S6).
DE analysis found 25 differential miRNAs between PCa patients with a Gleason score of ≥8 vs. those with a score of 7 (A vs. BC combined) (FDR p or p < 0.05, Table S5). The two Gleason score 7 status groups, B and C, differed by only 11 miRNAs (FDR p or p < 0.05, Table S5). We additionally checked the most DE miRNAs, 21, between all PCa patients (ABC combined, n = 30) and the healthy group (D, n = 10) (FDR p < 0.05, Table S5). Venn analysis indicated that over 80% of the changed miRNAs were unique to one of the comparisons, and none were common to all comparisons ( Figure 3A, Table S5).  According to gene set enrichment analysis for KEGG pathways, the DE miRNAs from these status group comparisons (A vs. BC, B vs. C) targeted mRNAs from a comprehensive set of cancer and progression-linked signalling, resistance, hormonal and pancancer pathways ( Figure 4A,B). A similar result was obtained from the ABC vs. D comparison. According to the numbers, we found 13 pathways between Gleason score ≥ 8 vs. 7 groups (A vs. BC), whereas the Gleason score of 7 status groups (B vs. C) differed in 64 pathways and PCa patients and the healthy group (ABC vs. D) by 69 pathways (p < 0.05, Table S7). However, overall >70% of the pathways from these comparisons overlapped with each other ( Figure 4C). They included the EGFR tyrosine kinase inhibitor and the endocrine resistance pathways and the PI3K-Akt, p53, FoxO, HIF-1, JAK-STAT, TNF, According to gene set enrichment analysis for KEGG pathways, the DE miRNAs from these status group comparisons (A vs. BC, B vs. C) targeted mRNAs from a comprehensive set of cancer and progression-linked signalling, resistance, hormonal and pan-cancer pathways ( Figure 4A,B). A similar result was obtained from the ABC vs. D comparison. According to the numbers, we found 13 pathways between Gleason score ≥ 8 vs. 7 groups (A vs. BC), whereas the Gleason score of 7 status groups (B vs. C) differed in 64 pathways and PCa patients and the healthy group (ABC vs. D) by 69 pathways (p < 0.05, Table S7). However, overall >70% of the pathways from these comparisons overlapped with each other ( Figure 4C). They included the EGFR tyrosine kinase inhibitor and the endocrine resistance pathways and the PI3K-Akt, p53, FoxO, HIF-1, JAK-STAT, TNF, TGF-beta, Wnt, MAPK, AMPK, Ras and Rap1 pathways, as well as the prolactin, neurotrophin, apelin and estrogen pathways (FDR p or p < 0.05, Figure 4A, Table S7). As expected, "miRNAs in cancer" and "prostate cancer" were among the top significant pathways ( Figure 4B). In addition, the targeted pathways included apoptosis, cellular senescence and choline metabolism as well as various viral pathways (Table S7).
We next performed qPCR to confirm the expression of selected DE miRNAs (n = 32 totally for all comparisons between groups A-C or A-C vs. D) and reference candidate miRNAs (n = 10, see methods) using all the samples. Here, in addition to miRNAs from the main comparisons (Table S5), we included some DE miRNAs from additional comparisons (Table S6), such as members of the miRNA 888-cluster-miR-891b, -888-5p, -892a, -892b, -891a-5p-differing particularly in comparisons that included group B (p < 0.05). In the PCA of the qPCR, despite some tendency of separation between PCa and healthy, a clear clustering of the status groups was not found ( Figure S3C).
For the comparison of status groups A, B and C, combined or individually, two miRNAs, miR-146a-5p and -892a, were significantly changed similarly to miRNAseq (FDR or p < 0.05 for pair-wise comparisons and ANOVA, Figure 5A, Table 2). The qPCR analysis additionally suggested dysregulation (1.3-4.2 fold) of miR-34a-5p, -424-5p, -892b, -1299 and -3065-5p between some of the groups, A-C (p < 0.05, Table S8), and 15 miRNAs between ABC vs. D, with the same changes as in the miRNAseq (FDR or p < 0.05, Tables S8 and S9). For nine of the 15 miRNAs, qPCR suggested changes also between some PCa groups, including miR-892a. We then conducted a ROC curve analysis using the qPCR data of the two most DE miRNA. Hsa-miR-892a, as a standalone marker for differentiating A from BC, produced an area under the curve of 0.894 (95% confidence interval 0.758-1.031, p = 1.46 × 10 −8 , Figure 6A), as well as 0.846 (0.546-0.981) specificity and 0.875 (0.473-0.997) sensitivity for the optimal cut-point. Similarly, hsa-miR-146a-5p produced an area under the curve of 0.770 (95% confidence interval, 0.587-0.953, p = 0.004, Figure 6A), as well as 0.900 (0.683-0.988) specificity and 0.600 (0.262-0.878) sensitivity for the optimal cut-point. Table 2. Fold changes of differentially expressed miRNAs between prostate cancer groups in qPCR. Table includes only the most significant miRNAs with statistically significant changes in both miR-NAseq analysis and in qPCR validation. Results marked in bold had a higher statistical significance (FDR p < 0.05 to <0.0001) than the rest (p < 0.05 to <0.0001). Fold change (FC), quantitative polymerase chain reaction (qPCR).      Table 2). Box-and-whiskers plots present the first quartile, median and third quartile, minimum and maximum values and outliers that are marked outside the whiskers range. Stars after p-values indicate the false discovery rate (FDR p value ** < 0.01, *** < 0.001). Delta quantification cycle (dCq).  Table 2). Box-and-whiskers plots present the first quartile, median and third quartile, minimum and maximum values and outliers that are marked outside the whiskers range. Stars after p-values indicate the false discovery rate (FDR p-value ** < 0.01, *** < 0.001). Delta quantification cycle (dCq).

Analysis of Prostate Cancer Progression Groups Uncovered Unique miRNA Signatures and Overlapping Cancer Progression-Linked Pathways
To understand which of the miRNAs in the uEV samples were linked to disease progression, we reanalysed the miRNAseq data using disease progression groupings formed based on clinical information 4-7 years post-RP (Table 1). Five patients had the most progressive disease during the follow-up time (PCa I group). Most of them (4/5) came from status group A (Table 1), including P30. In miRNAseq PCA, the P30 sample was among the non-clustering samples ( Figure S3D), whereas in PCA of the qPCR data, it clustered together with other samples ( Figure S3E). PCa group II included 3-4 patients from all status groups A-C (Table 1). Interestingly, the non-progressing group III also included patients from all the status groups, although eight out of 14 came from Group C (Table 1). Thus, the prior PCa status did not directly dictate the PCa progression.
In total, 151 unique miRNAs were DE between the PCa progression groups in the miRNAseq data (FDR p or p < 0.05, Table S10). Out of these, 59 (39%) and 66 (44%) miRNAs were DE between the most (I and II) compared to the least (III) progressed groups, respectively. Groups I and II differed in the expression of a lower number of miRNAs-23 (15%). Again, the proportion of miRNAs uniquely DE in different comparisons between progression groups was overall high-74% of the total ( Figure 3B). Interestingly, the miRNAs were also mostly different from those found in the status group analysis ( Figure 3C).
Despite the unique miRNAs, analysis for KEGG pathways again revealed most of the same cancer and progression-linked pathways found from previous comparisons between the PCa status groups or the PCa vs. healthy group (FDR p or p < 0.05, Table S11, Figure 4). Particularly, both progressed groups I and II differed from Group III of the non-progressors by 14 commonly known or potential cancer associated signalling, resistance and hormonal pathways ( Figure 4A), and altogether they shared 58 differential pathways ( Figure 4C). The pathways that differed uniquely only between Groups I and III, and not between other progression groups, included some hormonally regulated pathways such as those of relaxin, growth hormone, apellin and estrogen, as well as Rap1, mTOR and leukocyte pathways ( Figure 4A). Notably, out of all the comparisons, we found growth hormone, mTOR and leukocyte pathways here only. Group II differed uniquely from group III by Wnt, TNF, ErbB and neurotrophin terms. Analysis of KEGG pathways from other comparisons gave no pathways (I vs. II or I + II vs. III) or only one (I vs. II + III) significant cancer associated pathway (Table S11).
Cancers 2022, 14, Figure 6. ROC curves. Receiver operator characteristic (ROC) curve analysis using qPCR dat key differentially expressed miRNAs: (A) miR-892a and miR-146a-5p for separating Gleaso ≥ 8 (Group A) from the Gleason score 7 (B and C combined), as well as (B) miR-223-3p an 146a-5p for separating progressors (Groups I and II combined) from non-progressors (Gro Area under the curve (AUC). Table includes only the most significant miRNAs with statistically significant changes in bo NAseq analysis and in qPCR validation. Results marked in bold had a higher statistical sign (FDR p < 0.05 to <0.0001) than the rest (p < 0.05 to <0.0001). Fold change (FC), quantitative poly chain reaction (qPCR).  To further confirm that the DE uEV miRNAs between progression groups were cancerpathway regulating miRNA, we carried out overrepresentation analysis with the TAM 2.0 tool [33,34]. It confirmed that 32-61% of the DE miRNAs from all comparisons were identified as PCa related and overrepresented by 2.2-3.8-fold (Bonferroni adj p or FDR p < 0.01, Table 3). The top significant functions of the miRNAs differentiating groups I and II from III were linked to, e.g., epithelial-to-mesenchymal transition (EMT), inflammation, hematopoiesis and apoptosis (Bonferroni adj p or FDR p < 0.01). In addition, I differed from III (FDR p < 0.01) in the angiogenesis pathway. However, again no highly significant functions were found for the comparisons of I vs. II, or I and II combined vs. III.

MiRNA Name
Two of the 18 selected DE miRNAs from miRNAseq were confirmed by qPCR: miR-223-3p and -146a-5p (p < 0.05 for pairwise comparisons and/or ANOVA, Figure 5B). Both showed upregulation in the most progressed groups I and II relative to Group III (Table 2). In addition, qPCR suggested changes in 12 more miRNAs between some of the progression groups, including seven (miR-15a-5p, -16-5p, -31-3p, -31-5p, 92a-3p, -210-3p, 3065-5p) that were also DE between PCa status groups and/or the healthy group (Tables S8 and S9). ROC curve analysis using the qPCR data for differentiating progressors I and II from non-progressors III showed AUC 0.745 for miR-223-3p (p = 0.016) with a sensitivity of 0.714 (95% ci, 0.571-0.857) and a specificity of 0.786 (0.643-0.929) at the optimal cut-point ( Figure 6B). The AUC for miR-146a-5p was 0.768 (p = 0.003) with a sensitivity of 0.812 (0.544-0.960) and a specificity of 0.643 (0.351-0.872) in the optimal cut-point ( Figure 6B). Together, the findings suggested that uEV from different PCa status and progression groups display miRNAs regulating comprehensive and largely overlapping sets of cancer pathways. Despite the overlap, the individual miRNA signatures differentiating each of the PCa status or progression groups were unique. Table 3. Overrepresentation analysis using miRNAs differing between PCa progression groups. Table shows overrepresentation (fold) and miRNAs associating with the disease term "prostate cancer" and the most significant function term(s) (FDR p < 0.05, or p < 0.001, if no terms reached the FDR limit). Bonferroni adjusted p-value is also shown. Differentially expressed (DE).       Table 4). The EV quality was monitored with Western blotting, EM, total RNA profiles in Bioanalyzer Pico assays and sequencing using either the correlation study samples or healthy control samples collected and processed the same way ( Figure S4). UEV quality appeared similar to before [4,7,30] and to that of pEV in previous publications using ExoEasy [45,46]. We started the study by sequencing miRNAs in the samples from P33 with GS 7 disease, from whom we had collected urine and plasma pre-and post-RP. A healthy male (HC11, group D) and female donor (HC13) were included as technical or biological controls. In addition to uEV and/or pEV, we included plain plasma in this analysis. All samples passed the basic miRNAseq quality evaluation (Phread scores > 30) and resulted in an average of 0.9 M reads (range, 0.7-1.2 M, Table S12) and detection of 501, 407 and 633 miRNAs from uEV, pEV and plasma, respectively. MiRNA-146a-5p showed consistently higher expression in all of the pre-RP samples compared to post-RP samples or HC samples, with the exception of pEV of HC11 ( Figure S4). MiR-223-3p gave similar results from pEV and plasma, but it was not detected in the uEV samples. MiR-892a could not be detected in any sample using this platform.
The mRNAs were then compared with mRNA target lists for miR-146a-5p (n = 434), -892a (n = 427) and -223-3p (n = 370), obtained from miRWalk (Table S14). As 80%, on average, of the targets were expressed in any EV sample (Table S13), we restricted our As the miR-892a in particular had showed some differences in expression levels between status groups ABC and group D in the main study (Table S9), we additionally compared the expression of mRNA targets in the patient EV relative to HC11 EV (Table S16). We divided the mRNAs according to them showing a binary expression pattern (expressed in primary PCa samples, but not in post-RP samples of P33 or in HC11 samples), or mRNAs upregulated ≥5-fold relative to post-RP and ≥5-fold up-or downregulation relative to HC11. However, miR-892a targets in uEV were found in this analysis as often as the targets of miR-146a-5p and -223-3p (~11% of the listed mRNA targets for all three miRNAs). Thus, the correlation study suggested that uEV and pEV offer both unique and overlapping contents for studying both miRNA and their target mRNA in PCa pathogenesis.

Discussion
The quest for detecting or predicting aggressive prostate cancer and its progression has been ongoing for decades. The most sought-after goal involves non-invasive biomarkers, in which regard urine and, lately, uEV in particular have been investigated [8,11,49]. However, despite the interest, a complete picture is missing regarding cancer pathways that could be detected via uEV. Our analysis of uEV miRNA in prostate cancer status and progression groups, as well as target pathways and mRNAs, adds new insight regarding EV as a liquid biopsy.
Our results from uEV pointed to changes in unique sets of miRNAs in different PCa status and progression groups. These included several miRNAs that have been previously reported to regulate processes contributing to PCa development and metastasis, including both tumour suppressive and promoting roles in primary tumour-cancer associated fibroblasts, extracellular matrix, angiogenesis, EMT-and in migration, extravasation and colonization of metastatic sites [50]. Our pathway analysis consistently indicated that the DE miRNAs regulated the majority of the most important cancer progression promoting pathways and processes ( Figure 4, Table 3). This was evident when comparing both PCa status and progression groups and also between PCa and healthy groups, which supports the idea that uEV offer a good sample for studying miRNAs from these pathways. The high number of shared pathways was still a surprising finding given the unique miRNA signatures (Figures 3 and 4). The findings therefore suggests that specific miRNA signatures regulating the common pathways may uniquely associate with specific PCa disease and progression stage.
Pathways showing changes in all or most comparisons, and are therefore robustly associated with a higher Gleason score, advanced status or progression, included receptor tyrosine kinase -linked pathways (EGFR, PI3K-Akt, Ras, MAPK, JAK-STAT), p53, and the TGF-beta, AMPK and HIF-1 pathways. All these pathways have been linked to PCa progression and bone metastasis [51]. More unique changes were found in the mTOR and immune/leukocyte pathways (in I or II vs. III), several hormonal pathways (commonly changed in I vs. III) as well as the Wnt, TNF, ErbB and neurotrophin pathways (II vs. III or B vs. C, Figure 4, Table 3). Here, the association between the leukocyte pathways and aggressive progression is of interest, because infiltration of some T-cell subtypes has been shown to associate with a higher risk of PCa-induced death [52]. Out of the hormonal pathways, relaxin and apelin upregulation have been linked to metastasis or androgen independence [53,54]. For instance, Thompson et al. found that relaxin expression increased in RP specimens obtained after 6 months of androgen ablation, in androgen independent tumours and in bone metastases [53]. Local estrogen signalling may connect to PCa progression and development of hormone refractory disease [55]. Finally, previous work suggests that neurotrophins are particularly expressed in metastatic PCa, and-with EGF signalling-can mediate autocrine signalling, which is important for the progression of PCa [56]. In our results, brain-derived neurotrophic factor (BDNF) appeared as the target neurotrophic factor of miR-204-5p, downregulated upon progression (Tables S10 and S11). BDNF was also higher expressed in the pEV sample of P35 compared to the controls (Table S4). In agreement, loss of miR-204 results in BDNF/TRKB overexpression and activation of the Akt/mTOR/Rac1 signalling pathway, cancer cell migration and invasion in many cancers [57]. In PCa, overexpressed BDNF promotes progression via induction of EMT and anoikis resistance [58]. Early data suggests that neurotrophins induce neuregulin 1 (NRG1) release [59]. Interestingly, the neuregulin 1 secreted from cancer associated fibroblasts was recently shown to promote antiandrogen resistance [60]. This links our neurotrophin pathway finding to a possible PCa resistance mechanism. As neuregulin activates the HER3 and then the PI3K/Akt pathway [60,61], some additional pathways found in our study (e.g., PI3K/Akt signalling) could be associated to this chain of events. Thus, the EV results warrant further studies in men with advanced PCa treated with novel antiandrogens.
On the level of single miRNAs, the best candidates-miRNA-892a, miRNA-223-3p and miRNA-146a-5p-were all upregulated in patients with GS ≥ 8 or in patients progressing via the most aggressive disease course post-RP (Figures 5 and 6). In agreement with this, upregulation of the miRNA-888 cluster was previously shown in expressed prostatic secretion (EPS) urine exosomes and PC3 cell lines from high-grade prostate cancer and downregulation in lower grade cancer [15]. Upregulation of miRNA-223-3p was demonstrated in prostate cancer tissues and cell lines and was found to target SEPTIN-6 [62]. However, SEPTIN-6 mRNA expression appeared steady in all patients and EV types in our correlation study (Table S4). Both miR-223-3p and miR-888 clusters were dysregulated in semen EV from PCa patients relative to healthy controls [63]. Even if also upregulated in benign prostate hyperplasia, a combination of miR-223-3p with two other miRNAs and PSA was able to discriminate PCa patients from hyperplasia controls [63,64].
Interestingly, miR-888 cluster and miRNA-223-3p reside within X-chromosome, where particularly Xq27-28 within the HPCX1 locus has been associated with hereditary PCa and susceptibility to PCa, especially in Finland [15,65,66]. The study by Mattila et al. [66] screened Finnish families with a strong linkage to hereditary PCa and to this locus: miR-223 and also miR-146a (chromosome 5) were dysregulated in the lymphoblastoid cells from patients with hereditary PCa compared to healthy controls. Hence, we may have identified these particular miRNAs as we studied a cohort from a Finnish biobank.
Even if dysregulation of miR-146a has been widely linked to cancer, its upregulation in uEV may appear contradictory to its tumour suppressive functions and downregulation in PCa tissues [50,[67][68][69][70]. For example, miR-146a-5p downregulation was observed in castration resistant PCa tissues compared to androgen dependent PCa tissues, and its upregulation in PCa cells inhibited anchorage-independent growth, migration, invasion and angiogenesis via targeting the EGFR pathway or Rac1 [68,70]. The results from uEV and tissues/cells may differ, because uEV could serve as a disposal route for some tissue components, thus decreasing their quantity in tissues. Disposing of tumour suppressors could bring a growth advantage to PCa cells [71]. In support of this, Zhang et al. observed this kind of a difference in the tumour suppressor miR-15a-5p levels of hepatocellular carcinoma EV and cells [72]. Additionally, in the case of mRNAs, the transcripts in parental cells and their EV may differ [73]. However, prior data from urine has shown downregulation of miR-146a-5p [13]. The contradictory results could be explained by the different sample type (urine vs. uEV) and compared groups: we compared PCa status and progression groups, while the other study compared PCa to benign prostate hyperplasia patients [13]. The results could equally well differ due to choices concerning any (pre)analytical stepsstandardization of methods is missing, particularly for uEV [5,7,20,24,64].
Our study targeted miRNAs and mRNAs that were detected in small uEV-enriched samples isolated by ultracentrifugation, or pEV samples isolated with ExoEasy without further steps of purification or enzymatic treatments (RNAse/proteinase). These steps attempt to remove miRNA in other carriers than EV or miRNA associated with the molecular corona on the EV surface [23,27,74,75]. Our detected miRNAs can therefore be located in or on EV or in other small-sized miRNA carriers remaining in the samples, including recently identified supermeres [27,76]. Other limitations of the study include relatively small sample numbers in the individual groups, which affects statistical power, and modest fold changes for some miRNAs. We had limited success in confirmation of the DE miRNAs by qPCR in the main study or by using another miRNAseq platform in the correlation study. While the minute quantities of RNA in EV pose a challenge for miRNAseq, the qPCR validation problem may be due to lack of standard reference miRNAs for uEV and PCa studies. Thus, even if we strived to identify stable miRNAs in our and previous datasets, currently accumulating knowledge of uEV miRNAs hopefully leads to identification of better, widely accepted reference miRNAs in the future [20,21].
The complexity of EV miRNA research for the notoriously heterogeneous PCa could be helped by conducting larger or personalized studies targeting not only miRNA but also the mRNAs within EV. In this regard, our simple correlation study of three patients uncovered a notable number of mRNA targets for miR-146a-5p, -223-3p and -892a in uEV and pEV. This approach appears promising based on the high number of genes with previous implications in PCa pathogenesis. For example, PCa patient EV expressed genes dysregulated or altered in PCa and/or metastasis (e.g., FNLA, MMP16, CDON, NRP2, PAX5, SULT1B1, L1CAM, SCHBP1) [77][78][79][80][81][82][83][84], in castration resistant prostate cancer (DPY19L2, NOVA1, NOVA2) [85,86] and in enzalutamide resistant (SLC6A15) [87] or androgen independent PCa cells (CALN1) [88]. They also expressed tumour suppressors (e.g., SULF1, EBF3) [89,90], many protocadherin family members, which have both suppressor and oncogenic roles in PCa [91], and solute carrier family members (e.g., SLC35F1, SLC26A2) implicated in drug uptake and efficacy modulation [92]. It is clear that the significance of their EV expression in PCa remains to be elucidated-some of the targets were expressed at low levels. However, as we discovered interesting, but only partly overlapping, targets from the matched uEV and pEV (Tables 5 and S15), our study may help to focus future efforts with multiple patient and control types to the best EV source for selected miRNAs, mRNAs or their integrative analysis.

Conclusions
In conclusion, this study revealed that uEV samples contain miRNAs regulating well-known and emerging cancer pathways across the axis from cancer development to metastasis and to therapy resistance. We further identified cancer status and progressionassociated miRNAs that were located in the X-chromosome and/or had been previously linked to hereditary prostate cancer, especially among Finns. As patient EV expressed a notable number of PCa-associated mRNA targets for these key miRNAs, EV appear to provide candidate prognostic biomarkers to be further explored in large, personalized medicine trials.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/cancers14030532/s1, Figure S1: Western blotting of uEV samples, Figure S2: Quality of miRNAseq, Figure S3: Principal component analysis of miRNAseq and qPCR data, Figure S4. Quality control and miRNA expression in the correlation study, Figure S5: Western blotting originals, Table S1: Raw count data from miRNAseq, Table S2: Quantitative PCR assay information for all assays and Normfinder stabilities in miRNAseq for the selected reference miRNAs, Table S3: Raw count data from miRNAseq in the correlation study, Table S4: Raw count data from mRNAseq in the correlation study, Table S5: Differentially expressed miRNAs between the main prostate cancer status groups and healthy controls from miRNAseq, Table S6: Additional comparisons of differentially expressed miRNAs between prostate cancer status groups and healthy controls from miRNAseq, Table S7: KEGG pathway analysis with differentially expressed miRNAs between prostate cancer status groups and healthy controls from miRNAseq, Table S8 Fold changes of differentially expressed miRNAs between prostate cancer groups in qPCR, Table S9: Fold changes of differentially expressed miRNAs between prostate cancer status groups and healthy controls in qPCR, Table S10: Differentially expressed miRNAs between prostate cancer progression groups from miRNAseq, Table S11: KEGG pathway analysis with differentially expressed miRNAs between prostate cancer progression groups from miRNAseq, Table S12: Read and miRNA counts in the miRNAseq of the correlation study, Table S13: Read, mRNA and miRNA targeted mRNA counts in the mRNAseq of the correlation study, Table  S14: Target mRNAs of miR-146a-5p, -223-3p and -892a, Table S15: Expression of target mRNAs of miR-146a-5p, -223-3p and -892a in uEV and pEV from correlation study patients, Table S16: Binary, higher or lower expression of miR-146a-5p, -223-3p and -892a target mRNAs in uEV and pEV from correlation study patients compared to controls.  The raw sequencing data analysed in this study is not publicly available due to local biorepository regulations. We provide raw miRNAseq and mRNAseq count data in the Tables S1, S3 and S4-the data from the healthy control group of the main study has been analysed in different context earlier [20].