From the Broad Institute GDAC (Genome Data Analysis Center) (Broad Institue, 2021), we extracted all available Gene expression (mRNA-Seq) data of PRAD (a total of 550 cases), containing both tumoral (T) and non-tumoral (NT) tissue samples at preprocess level [prostate PRAD (cancer type), RNASeqV2 (data type), level 3 (archive type) and 2016-02-13 (data version)]. Data were generated based on Illumina HiSeq 2000 platform and annotated to reference transcript set of UCSC hg19 gene standard track. One single sample was removed from the set to prevent possible disturbances in the results as it corresponded to a metastatic sample. In addition, a total of 531 Isoform Expression Quantification (miRNA-Seq) files containing both tumoral (T) and non-tumoral (NT) tissue samples, as well as clinical data for each patient/sample was obtained from TCGA data portal (NIH, 2021a). All data are controlled; the access has been requested through the GDAC of the National Institutes of Health (NIH).

Differential expression analyses have been carried out using edgeR (version 3.28.0) Bioconductor package (Robinson et al., 2009; McCarthy et al., 2012). Quasi-likelihood F-test (QL) from the GLM (Generalized Linear Models) was used to determine differentially expressed genes (DEG) related to tumor aggressiveness and treatment effectiveness. For that purpose, two different analyses were performed, in which different groups were established based on the clinical information of each case:

1) DEG related to AD (androgen deprivation) therapy response: Three categories were established [DR1: treatment based on a single drug target (43 cases)], DR2: treatment in which a new drug target has been prescribed due to failure of the first one (25 cases), DR3: chemotherapy (5 cases). Four drug targets apart from the chemotherapy were considered: LHRH agonists, LHRH antagonists, antiandrogens, and CYP17 inhibitors. NT (non tumor) samples, as well as those lacking information on treatment, were excluded from this analysis.

To determine differentially expressed (DE) miRNA related to tumor aggressiveness, three categories were defined [G0 = NT samples (51 cases), G1 = Gleason score equal or lower than 7 (285 cases), G2 = Gleason score higher than 7 (195 cases)]. As for mRNA-Seq, Quasi-likelihood F-test (QL) from the GLM (Generalized Linear Models) was used to perform the analysis.

Low expression filter was applied for every single analysis by using as representative threshold the number of samples of the group that has less expression values. The normalization of the samples was calculated using the “calcNormFactors” function and the trimmed mean of M-values (TMM) method (Robinson and Oshlack, 2010), while the dispersions were estimated using the “estimateGLMCommonDisp,” “estimateGLMTagwiseDisp,” and “estimateGLMTrendedDisp” functions. Three contrasts were carried out in each analysis:

- DR2 vs. DR1, DR3 vs. DR1, and DR3 vs. DR2 for AD therapy response (mRNA-Seq).

In all our DE analyses, the p value was adjusted by Benjamini–Hochberg false discovery rate (FDR) procedure (Benjamini et al., 2001).

We used STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) database (Szklarczyk et al., 2019) to select DEGs along with IGF2. To do so, we performed the IGF2 interactor search according to the default parameters within “single protein by name” option and selecting no more than 20 interactors in first shell. From all of them, we choose those that exceeded a score of 0.980. Additionally, IGF2 interactors with particular significance in PC such as VEGFA, STAT3, FLT1, KDR, NRP1, NRP2, and HIF1A were manually added to complete 20 interactors.

On the other hand, the selection of miRNAs of interest was carried out using MirTarBase (miRTarBase, 2021) through the search of IGF2 by target gene according to default parameters.

SNVs (single nucleotide variants)/mutations may affect gene function occasionally leading to a total loss of function (LOF). This is related to the development and prognosis of a tumor. For this reason, we study the alterations at the mutation level of IGF2 gene in TCGA-PRAD cohort. To this end, we obtained 503 annotated somatic mutation files (MuTect2 annotation type) corresponding to TCGA-PRAD from TCGA data portal (NIH, 2021a). First, each of the Variant Call Format (VCF) file was parsed to extract information regarding the mutated genes presented in each sample. In a subsequent step, the result of each VCF file was checked against IGF2 gene.

Out of the IGF2 gene SNVs annotated with clinical association in public databases such as The National Center for Biotechnology Information (NCBI, 2019), we selected two (rs1004446 and rs3741211) for the present study. Only those SNVs with an allele frequency higher than 20% on the minor allele (MAF) in the Caucasian population according to the Ensembl database (Ensembl, 2021) were taken into consideration, more details of the probes can be found in Supplementary Table S2.

Each buccal swab or blood sample DNA was extracted using organic extraction reagents (1 ml de Stain Extraction Buffer + Proteinase K). DNA extraction protocol was performed as described by Freeman et al. (2003) and optimizations developed by Gomez-Martín A. et al. (Gómez-Martín et al., 2015). All samples were standardized to 20 ng/μl using Nanodrop 2000/2000c (ThermoFisher, United States) quantification. DNA genotyping was performed using TaqMan® Genotyping Master Mix (Applied Biosystems, United States) which included all essential components (except probes, templates and water) for polymerase chain reaction (PCR). Allelic discrimination assays were carried out in a 7900HT Fast Real-Time PCR System (Applied Biosystems, United States). Results were analyzed using SDS software v.2.4 (Applied Biosystems, United States).

mRNA from a total of 78 fresh tissue samples was extracted using Trizol®/chloroform method and quality validated by A260/A280 in NanoDrop™ 2000c. Only those samples with the best quality and including all clinical records were selected. Reverse transcription was performed with TaqMan™ Advanced mRNA cDNA Synthesis kit (Applied Biosystem, Foster City, CA). Quantitative polymerase chain reaction (qPCR) was performed with SYBR Green designed probes (Life Technologies, Carlsbad, CA), on a 96-wells plate with QuantStudio 6 Flex Real-Time PCR System (Applied Biosystems). qPCR reactions were performed as follows: 95°C during 10 min for enzyme activation; followed by 40 cycles of 15 s at 95°C and 1 min at 60°C for denaturing and annealing/extension. Primers were designed using Primer-Blast (Ye et al., 2012) (NIH) software under the following conditions: they must span an exon-exon junction, have a PCR product size between 60–150 nucleotides and have a primer melting temperature within the range of 59–61°C. Sigma Aldrich company designed the primers with the following sequences UG_GX_IGF2_f (Forward): CGC​TGT​TCG​GTT​TGC​GAC, and UG_GX_IGF2_r (Reverse): GGA​TTC​CCA​TTG​GTG​TCT​GGA.

All samples were run in triplicates, with a NTC (non template control) in each plate. Threshold cycles (C_T) ≥ 35 were considered as undetermined values. mRNAs expression levels were quantified using the comparative threshold cycle (Ct) method (2^−ΔΔCt) relative to HPRT1 (hypoxanthine phosphoribosyltransferase 1) expression as an endogenous control. Firstly, difference between IGF2 and HPRT1 expression was calculated for each sample (ΔC_T = C_{T IGF2} ⁻ C_{T HPRT1}). Normalization was done using the mean of reference group; treatment sensitivity group for therapy response analysis; and Gleason ≤7 group for aggressiveness study (ΔΔC_T = ΔC_T ⁻ ΔC_{T reference}). Relative quantification parameter (RQ or 2^−ΔΔCt) was estimated for each case and used in statistical analysis.

Total RNA of 159 fresh tissue biopsies (including patients and controls) were extracted using Trizol®/chloroform method and quality validated by A260/A280 in NanoDrop™ 2000c. Reverse transcription was performed with TaqMan™ Advanced miRNA cDNA Synthesis kit (Applied Biosystem, Foster City, CA). Quantitative polymerase chain reaction (qPCR) was performed with TaqMan™ probes (Life Technologies, Carlsbad, CA), according to manufacturer’s protocol, on a 96-wells plate with QuantStudio 6 Flex Real-Time PCR System (Applied Biosystems). qPCR reactions were performed as follows: 5°C during 20 s for enzyme activation; followed by 40 cycles of 1 s at 95°C and 20 s at 60°C for denaturing and annealing/extension. For liquid biopsy analysis, plasma of 60 samples was isolated from blood; this process was carried out, at most, 4 h after collection. Total RNA of the samples was extracted using the miRNeasy Serum/Plasma Kit (Qiagen GE). All samples were run in triplicate, with a NTC in each plate.

We included the analysis of miR-93-5p, as one of the most interesting miRNAs according to bioinformatic analysis. miRNAs expression levels were quantified using the comparative Ct method (2^−ΔΔCt) relative to RNU6B (U6B small nuclear RNA) expression as an endogenous control.

Out of all the DEGs obtained in differential expression analyses, we focused on IGF2 and some of its closest interactors obtained by using the STRING database. The protein-protein interaction network obtained can be seen in Figure 1.

According to the obtained results by G0 vs. G1, G0 vs. G2, and G1 vs. G2, the most interesting genes in terms of p-values and FDR are NRP2, KDR, and IGF2. These results can be seen in Table 2. Data obtained from differential expression analysis based on AD therapy response are shown in Supplementary Table S3. Except for IGF2, which is under-expressed in patients with treatment resistance, any other gene showed any statistically significant result.

MiRTarBase (mirTarBase, 2021) is one of the most comprehensively annotated and experimentally validated miRNA–target interaction databases. According to this database, there are 46 miRNAs which may target IGF2.

By the search of these miRNAs, according to Gleason-based differential expression analysis performed on the TCGA-PRAD miRNA database, it was determined that miR-93-5p and miR-200c-3p are those which could have a potential influence on IGF2 modulation in PC. However, none of them showed any statistical significance in G1 vs. G2. More details can be observed in Table 3.

According to the in silico analysis and TCGA comparisons, we have obtained several remarkable data in the interactions of IGF2 with other pathogenic effect variants. That is the case of rs1114167321, rs553443857, rs1057518115, rs1064794050, and rs869320620. When conducting somatic analysis, rs758164144 is the most frequent variant in G1 cluster (Gleason scores ≤7) and, even with low presence, rs3842753 is only present in Gleason scores >7. See more details in Table 4. rs1004446 has also located as a somatic mutation in G1 clustering, data not shown.

For IGF2 gene, both SNPs (rs1004446 and rs3741211) were linked with a statistic r = 0.9798, therefore we have only analyzed one of them with TaqMan probes (i.e., rs1004446).

In relation to aggressiveness, we compared Gleason scores ≥ or <7; as well as the value of D’Amico risk. None of them showed statistical significance. In the case of treatment response classification, we grouped patients according to sensibility or resistance to treatment, but unfortunately, we could not prove any data with statistic power (Supplementary Table S4).

Results achieved by qPCR of fresh tissue samples were analyzed following ΔΔCt method and using a non-parametric test (Mann Whitney). The level of statistical significance used was p < 0.05. The value of genetic expression of each patient was calculated as the average ±SD of three different replicates. A Tukey’s range test was performed to detect anomalous values. To increase the statistical significance and verify our results tendency, analysis was repeated including all replicates from patients as individual values. As can be seen in Figure 3, when comparing aggressiveness, we found similar statistically significant patterns as in the TCGA analysis. Although when analyzing treatment response, we could not observe any significant differences, we can see the same patterns that those observed in bioinformatic analysis.

We found by experimental analysis in plasma and tissue samples, that when comparing G1 vs. G2, miR-93-5p is over-expressed according to aggressiveness (Gleason score), the same patterns are repeated with TCGA data. In Figure 4, we can see how miR-93-5p expression changes according to Gleason score. Furthermore, we found that miR-93-5p follows the same expression patterns in both plasma and tissue samples.