Role of IGF2 in the Study of Development and Evolution of Prostate Cancer.

Prostate Cancer (PC) is commonly known as one of the most frequent tumors among males. A significant problem of this tumor is that in early stages most of the cases course as indolent forms, so an active surveillance will anticipate the appearance of aggressive stages. One of the main strategies in medical and biomedical research is to find non-invasive biomarkers for improving monitoring and performing a more precise follow-up of diseases like PC. Here we report the relevant role of IGF2 and miR-93-5p as non-invasive biomarker for PC. This event could improve current medical strategies in PC.

Introduction

It is well known that prostate cancer (PC) is a heterogeneous disease, which makes it difficult the identification of any clinical and molecular biomarker in disease management. PC is one of the most frequently diagnosed tumors among men in Europe and reaches the second position when comparing data worldwide, with over 1.4 million diagnoses recorded in 2020 (WHO, 2020). The use of robust biomarkers; mainly focused on molecular non-invasive ones; is still a challenge in this tumor. Several germline variants have been suggested as relevant in PC such as those in ATM (ataxia-telangiectasia mutated), BRCA1 (breast cancer), BRCA2, MSH2 (MutS Homolog 2), MLH1 (mutL homolog 1), MSH6 (MutS Homolog 6), PMS2 (PMS1 homolog 2), EPCAM (epithelial cellular adhesion molecule) and HOXB13 (Homeobox B13) genes (Saunders et al., 2021). Additionally, recent data support the role of several SNPs in IL-6 (Interleukin 6) gene (rs1800795, rs1800796 and rs1800797) as biomarkers of an increased cancer risk in several tumors. Specially, variants rs1800795 and rs1800796 are associated with an overall increased risk of PC (Harun-Or-Roshid et al., 2021).

Here, we focus on the role of IGF2 (insulin-like growth factor 2) as a novel marker for PC management. IGF2 encodes a member of the insulin family of polypeptide growth factors, which are involved in development, cancer biology and growth. It binds to type-1 insulin-like growth factor receptor (IGF1R), which activates downstream members of the PI3K (phosphatidylinositol 3-kinase)/AKT (alpha serine/threonine-protein kinase) and MAPK (mitogen-activated protein kinases)/ERK (extracellular signal-regulated kinases) pathways. This gene is important for cell survival and tumorigenesis; moreover it has been suggested that IGF2, in combination with SSTR2 (Somatostatin Receptor 2), plays an important role in PC survival. Previous data have focused on the role of IGF2 expression patterns or epigenetic imprinting in PC (Cao et al., 2014; GeneCard, 2021); and other tumors such as colon cancer, detecting recurrent fusions in this gene (Yun et al., 2020). Furthermore, IGF2 messenger RNA binding protein 3 (IMP3) has been reported to be over-expressed in PC and strongly correlated to poor prognosis. The main role of IMP3 (U3 small nucleolar ribonucleoprotein) has been included in PI3K/AKT/mTOR signalling pathway (Zhang et al., 2020). The relevance of IGF2 in PC has been denoted not only in mRNA, but also in lncRNA and SNPs. That is the case of IGF2AS (IGF2 Antisense RNA), which has been proposed as an epigenetic tumor suppressor in human PC. Moreover, the relationship between IGF2AS/and IGF2 has been included as a possible marker for future therapeutic targets in PC treatment or gastric cancer (Chen et al., 2019; Xing et al., 2021). Both, IGF2 and its receptor IGF1R constitute desirable therapeutic targets; mainly due to these evidences showing that targeting either IGF2 or its receptor IGF1R, blocks cancer progression and displays significant antitumor activity (Xing et al., 2021).

There are several SNPs in IGF2 that have been previously reported to have a role in cancer or other diseases; such as rs1004446. This SNP has been previously proposed as a marker of decreased endometrial cancer risk (McGrath et al., 2011); or increased PC risk (Cao et al., 2014). Another IGF2 SNP, rs4320932, has been associated with a decreased risk of ovarian cancer in G allele carriers (Pearce et al., 2011).

Several miRNAs have been also evaluated with a relevant role in IGF2 regulation. That is the case of miR-100 and miR-125b, which play a proven tumor suppressor role in hepatocellular carcinoma, by inhibiting IGF2 expression and activating AKT/mTOR pathway (Seol et al., 2020). miR-141 down-regulation blocks VEGF (Vascular Endothelial Growth Factor) and IGF2 expression; and also interacts with osteoblasts proliferation, which is relevant in osteosarcoma (He et al., 2016). In pancreatic cancer, it has also been proven the role of miR-141 in IGF2BP2 (IGF2 mRNA-binding protein 2); which is known to play oncogenic roles. Genomic amplification and silencing of miR-141 also contribute to IGF2BP2 activation; opening promise molecular targets in pancreatic tumors (Xu et al., 2019).

Concerning to somatic mutations, there is scarce data in PC, mainly due to current methodological strategies limitations in their analysis of because of their location in non-coding regions. In PC, their effects on driving tumorigenesis and progression have not been systematically explored (Wang and Li, 2021). We have previously published the role of several somatic mutations, just discovering an incipient role of c.1621A > C (rs3822214) in KIT (tyrosine kinase), c.38G > C (rs112445441) in KRAS (kirsten rat sarcoma virus) and c.733G > A (rs28934575) in TP53 (tumor protein) genes among patients with PC; although with weak associations (Martinez-Gonzalez et al., 2018). Others authors suggested the association of increased expression patterns of homeobox B13 (HOXB13), a gene related to normal prostate development, with worse outcomes after PC surgery (Weiner et al., 2020).

Here, we perform an integrated analysis combining bioinformatic and experimental analyses proving the role of IGF2 as a marker in PC. We focus on two main SNPs and miRNAs interactions in PC. The use of both biomarkers (SNPs and miRNAs) could be easily developed in clinical routine practice, mainly by the low cost of these methodologies.

Materials and Methods

Study Population for Experimental Analysis

Present study includes data from 199 men with prostate specific antigen (PSA) values above 4 ng/ml and histological confirmed PC. A total of 30 EDTA (ethylenediamine tetraacetic acid) blood samples and 38 buccal swabs from these subjects were collected by the Urology Service from “Hospital Universitario Virgen de las Nieves de Granada, Spain.” Samples were stored at −20°C until they were processed. Both, blood samples and buccal swabs, were used for SNPs genotyping analysis.

For RNA expression analysis, 131 fresh tissue samples were collected from nearly 66% of the patients of present study; these samples were stored at −80°C until they were processed. Concerning mRNA analysis, just 78 samples (39.2% of the total samples) were available (mainly limited by the quality of mRNA). Moreover, for miRNAs analysis we included all 131 fresh tissue samples of PC patients and 28 controls. Several clinical data of samples were collected such as age, Gleason Score, minimum PSA value, and treatment follow-up (Table 1). All study participants provided a written informed consent before being enrolled, and the study was previously approved by the Research Ethics Committee of Granada Center (CEI-Granada internal code 1638-N-18) following Helsinki ethical declaration. A supplementary figure explaining this sample distribution is included in Supplementary Figure S1.

TABLE 1

Descriptive variables of PC samples.

Variables	PC n* (%)
Age (years)
<60	5 (7.69%)
60–69	19 (29.23%)
70–79	30 (46.16)
≥80	11 (16.92%)
PSA level (ng/ml)
<20	37 (47.44%)
≥20	41 (52.56%)
Gleason score
≤7	44 (53.01%)
>7	39 (46.99%)
D’Amico Risk Classification
Low	8 (10.53%)
Medium	20 (26.31%)
High	48 (63.16%)
Treatment Response
Sensitivity	35 (38.04%)
Resistance	57 (61.96%)
Metastasis
Yes	45 (48.91%)
No	47 (51.09%)

n*(some reports data are missing; for that reason, the total number of samples do not sum the same total).

Bioinformatic Analysis

This analysis was performed by the access to “The Cancer Genome Atlas (TCGA)” which was initiated in 2005 and, as of today, it has over 2.5 petabytes data of 20,000 primary cancers and matched normal samples from 33 cancers types. TCGA was created as an easy way of exploring the entire spectrum of genomic changes involved in human cancer (Cancer Genome Atlas, 2021). This makes TCGA repositories an extraordinary-value source of data in studies like the present one.

TCGA Data of Prostate Adenocarcinoma

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).

TCGA Differential Expression Analyses

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.

2) DEG related to Gleason score: Three categories were defined [G0 = NT samples (52 cases), G1 = Gleason score equal or lower than 7 (292 cases), G2 = Gleason score higher than 7 (205 cases)].

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).

- G1 vs. G0, G2 vs. G0, and G2 vs. G1 for Gleason score (mRNA-Seq).

- G1 vs. G0, G2 vs. G0, and G2 vs. G1 for Gleason score (miRNA-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.

TCGA Somatic Mutations

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.

In Silico Analysis

Based on data available in Genome Browser of UCSC (University of California, Santa Cruz), we obtained a total of 60 different SNPs of IGF2 gene (Kent, 2002). An analysis of these IGF2 variants was carried out in “The variant effect predictor” (McLaren et al., 2016). This software was used to calculate changes in transcripts and malignancy of variants. We also used ClinVar tool (Landrum, 2018) for data validation (only 6 of the 60 SNPs were available in this software). These analyses were performed to confirm the interactions in this gene with other germline or somatic variants, as well as miRNAs. Main results are included in Supplementary Table S1.

Functional Analysis

Pathway analysis in IGF2 gene was evaluated using DAVID (Database for Annotation, Visualization and Integrated Discovery) Bioinformatics Resources v6.8 (https://david.ncifcrf.gov/, accessed on November 1, 2021) to obtain the role of gene pool, clinical implication, ontology and involved metabolic pathways. Moreover, STRING search tool was used to calculate Interacting Genes with our target and interaction among our target genes (https://string-db.org/, accessed on November 1, 2021) (Szklarczyk et al., 2019).

Molecular Analyses

This section pretends to validate main results obtained by bioinformatic analysis. We performed several analytic procedures contains SNPs genotyping and mRNA and miRNA expression analysis.

SNPs Selection and Genotyping

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 Analysis

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 (CT) ≥ 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 (ΔCT = CT IGF2 − CT HPRT1). Normalization was done using the mean of reference group; treatment sensitivity group for therapy response analysis; and Gleason ≤7 group for aggressiveness study (ΔΔCT = ΔCT − ΔCT reference). Relative quantification parameter (RQ or 2−ΔΔCt) was estimated for each case and used in statistical analysis.

miRNA 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.

Statistical Analysis

All analyses were performed using the SPSS v.22 statistical package (IBM Corporation, United States). Relationships between different genotypes and clinical variables were studied using the χ2 test. Odds Ratios (OR) and 95% confidence intervals (95% CI) were calculated by binary logistic regression. Genotypes analyses, Hardy–Weinberg equilibrium and Linkage disequilibrium (LD) analyses were performed using the online SNPStats software (Solé et al., 2006). SNPs are considered to be in LD when they have a value of r2 > 0.5. Present SNPs were in LD. For expression analysis, Shapiro-Wilks test was used to test the normalization of the samples. This test revealed that our results did not follow a gauss distribution, therefore a non-parametric test (U-Mann Whitney test) was performed for all variables. The level of statistical significance used was p < 0.05.

Results

The main features of the study population stratified by PC and treatment response are shown in Table 1. First, we conducted a bioinformatic analysis and most relevant data, such as rs1004446 (IGF2) and mir-93-5p, were selected for being 199 samples (more details in Supplementary Figure S1).

Bioinformatic Analysis (TCGA)

Gene Expression Analysis

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.

FIGURE 1

STRING network (default k-means clustering method) performed by the introduction of IGF2 with 20 interactors.

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.

TABLE 2

Differential Expression of IGF2 and IGF2 interacting protein genes using Gleason score.

Gene	LogFC	FDR	p value	Test
NRP2	−1.3309	1.08E-18	4.55E-20	G0 vs. G1
IGFBP2	0.9156	2.74E-09	5.53E-10	G0 vs. G1
KDR	−0.8654	2.98E-09	6.02E-10	G0 vs. G1
IGFBP3	−0.6879	1.08E-06	3.27E-07	G0 vs. G1
IGFBP4	−0.5516	4.89E-05	1.95E-05	G0 vs. G1
IGF2R	0.4125	0.0002	9.81E-05	G0 vs. G1
IGF2	−0.9725	0.0008	0.0004	G0 vs. G1
HIF1A	−0.3241	0.0015	0.0008	G0 vs. G1
IGFBP5	−0.4570	0.0022	0.0012	G0 vs. G1
IRS1	−0.4073	0.0042	0.0023	G0 vs. G1
VEGFA	−0.5247	0.0145	0.0089	G0 vs. G1
IGF1R	0.2770	0.0411	0.0278	G0 vs. G1
NRP2	−1.0033	2.46E-10	4.61E-11	G0 vs. G2
KDR	−0.7398	1.28E-06	4.17E-07	G0 vs. G2
IGFBP4	−0.6542	2.88E-06	9.91E-07	G0 vs. G2
IGFBP2	0.6512	2.29E-05	9.07E-06	G0 vs. G2
NRP1	0.5826	0.0003	0.0001	G0 vs. G2
IGFBP6	−0.7334	0.0007	0.0004	G0 vs. G2
IRS1	−0.4930	0.0007	0.0004	G0 vs. G2
IGF2	0.9849	0.0045	0.0026	G0 vs. G2
HIF1A	−0.2941	0.0057	0.0033	G0 vs. G2
VEGFA	−0.4553	0.0431	0.0296	G0 vs. G2
IGF2	1.9574	2.63E-23	1.13E-25	G1 vs. G2
IGFBP3	0.6200	2.61E-11	6.90E-13	G1 vs. G2
NRP1	0.6150	4.31E-11	1.18E-12	G1 vs. G2
IGFBP5	0.3610	0.0002	4.10E-05	G1 vs. G2
IGF2R	−0.2329	0.0008	0.0002	G1 vs. G2
IGFBP6	−0.4501	0.0017	0.0004	G1 vs. G2
NRP2	0.3275	0.0017	0.0005	G1 vs. G2
IGFBP2	−0.2644	0.0043	0.0013	G1 vs. G2

Gene, Gene symbol; logFC, logarithmic fold change; FDR, False Discovery Rate; Test, contrast. Here, comparisons were developed, including tissue samples (549 cases), comparing G0 = NT samples (52 cases), G1 = Gleason score equal or lower than 7 (292 cases), G2 = Gleason score higher than 7 (205 cases).

miRNA Expression Analysis

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.

TABLE 3

IGF2 target miRNAs analysis comparing Gleason scores.

miRNA	logFC	p value	FDR	Test
miR-93-5p	1.6018	1.8904e-37	2.5096e-35	G0 vs. G1
miR-200c-3p	1.56318	2.1194e-32	1.2504e-30	G0 vs. G1
miR-100-5p	−0.53748	3.6629e-09	1.5315e-08	G0 vs. G1
miR-320a	0.6363	2.3654e-08	8.5443e-08	G0 vs. G1
let-7a-5p	0.4601	3.5651e-07	1.0942e-06	G0 vs. G1
miR-339-3p	−0.2659	0.0044	0.0076	G0 vs. G1
miR-125b-5p	−0.2182	0.0055	0.0093	G0 vs. G1
miR-200b-3p	0.4704	0.0112	0.0180	G0 vs. G1
miR-320b	0.3298	0.0179	0.0277	G0 vs. G1
miR-93-5p	1.9519	5.3193e-48	9.4151e-46	G0 vs. G2
miR-200c-3p	1.6185	4.8291e-33	1.6027e-31	G0 vs. G2
miR-100-5p	−0.6493	7.6119e-12	2.9940e-11	G0 vs. G2
miR-125b-5p	−0.4584	1.7001e-08	4.8534e-08	G0 vs. G2
let-7a-5p	0.5252	2.0643e-08	5.8619e-08	G0 vs. G2
miR-320a	0.5621	1.2643e-06	2.9706e-06	G0 vs. G2
miR-3200-3p	0.6489	0.0004	0.0007	G0 vs. G2
miR-200b-3p	0.5366	0.0052	0.0079	G0 vs. G2
miR-429	0.5581	0.012	0.0178	G0 vs. G2
miR-320b	0.3375	0.0187	0.0261	G0 vs. G2
miR-150-5p	0.3765	0.0305	0.0413	G0 vs. G2
miR-93-5p	0.3502	2.2858e-08	3.0344e-07	G1 vs. G2
miR-125b-5p	−0.2401	1.2726e-06	1.1855e-05	G1 vs. G2
miR-3200-3p	0.3295	0.0012	0.0045	G1 vs. G2
miR-339-3p	0.1654	0.0044	0.0143	G1 vs. G2

miRNA, miRNA symbol; logFC, logarithmic fold change; FDR, False Discovery Rate; Test, contrast. Here, comparisons were developed, including tissue samples (531 cases), comparing G0 = NT samples (51 cases), G1 = Gleason score equal or lower than 7 (285 cases), G2 = Gleason score higher than 7 (195 cases).

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.

TABLE 4

Summary of the main somatic mutations in IGF2 (TCGA-PRAD cohort).

Mut position	Mut id	Total	G1 Gleason score ≤7	G2 Gleason score >7
2133567	rs758164144	7	6	1
2136949	rs3213216	4	2	2
2159830	rs3842753	2	0	2
2160994	rs689	2	2	0

A functional analysis in IGF2 and its 20-interactors genes was performed IGF2 using STRING (Szklarczyk et al., 2019) and DAVID (Database for Annotation, Visualization and Integrated Discovery) Bioinformatics Resources v6.8 (DAVID, 2020), to obtain the role of gene pool, clinical implication, ontology and involved metabolic pathways. As a result, we found that Proteoglycans in cancer pathway is the most enriched one according to p-value (4.6e-08) and FDR (1.6e-06) values. A simple diagram of the IGF2 pathway can be seen in Figure 2.

FIGURE 2

IGF2, IGF1 and insulin bind their specific receptors, which include IGF1R, IGF2R, IR, and hybrid receptors. Ligand binding results in autophosphorylation of the tyrosine residues of each receptor, leading to recruitment of the adaptor proteins IRS and Shc to the intracellular domains of the receptor’s β-subunits. This process activates different signalling cascades through the PI3K-AKT and RAS/RAF/MEK/ERK/ERK pathways, resulting in stimulation of translation and cell cycle progression, increased proliferation and growth, and inhibition of apoptosis.

Molecular Analysis

Once the bioinformatic analyses were performed, we tested the obtained results by molecular analysis in blood samples of our population described in Table 1.

Linkage Disequilibrium Analysis

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).

Association of rs1004446 (IGF2) Genotype With Aggressiveness and Treatment Response

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).

IGF2 Gene Expression Analysis by qPCR

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.

FIGURE 3

IGF2 expression analysis comparing aggressiveness (A) and treatment response (B). Mean is represented by a plus symbol. 2−ΔΔCt (mean ± SD): G1 = 1.185 ± 1.247; G2 = 1.642 ± 1.542; Sensitivity = 1.8710 ± 2.150; Resistance = 0.8775 ± 0.534.

miRNAs 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.

FIGURE 4

miR-93-5p expression analysis comparing Gleason score. dCt (mean ± SD): Gleason score 0 (NT controls) = 10.258 ± 1.129; Gleason score 6 = 7.692 ± 1.643; Gleason score 7 = 8.280 ± 3.209; Gleason score 8 = 4.622 ± 2.096; Gleason score 9 = 5.477 ± 2.156.

Discussion

Here, we focus on the aim of reinforcing the interesting role of bioinformatic analysis using TCGA database for searching genetic markers in PC. First of all, according to our bioinformatic analysis in PRAD (TCGA), IGF2 is denoted as one of the most expressed gene in prostate tissue samples. The physiological roles of IGF2, as well as its dependence on GH (growth hormone) production, are still controversial. Both IGF1 and IGF2 activate a common receptor, the IGF1 receptor (IGF1R), which stimulates mitogenic signals, antiapoptotic and pro-survival activities (Werner et al., 2021). Furthermore, an over-expression of IGF2BP2 (mRNA binding proteins 2) has been associated with a poor prognosis of the disease in multiple human cancers, as well as, with a shorter survival and poor prognosis in acute myelocytic leukemia, low-grade gliomas, breast , esophageal, hepatocellular, head and neck squamous cell, pancreatic ductal adenocarcinoma and gallbladder carcinomas (Wang et al., 2021). Also, IGF1R inhibitors are suggested as anti-cancer drugs because of their effects on proliferation inhibition (Tsui et al., 2021).

Although there are not many data in PC, there are reports about the role of IGF2-mRNA and its peptide in PC, with a decrease of 80% in PC compared to non-neoplastic adjacent prostate (Kingshott et al., 2021). There are also data, including the over-expression patterns of IMP3 (messenger RNA binding protein 3 related to IGF2) in PC, related to patients’ poor prognosis (Zhang et al., 2020). The results obtained in the present work reveal a significant increase of IGF2 expression in patients with Gleason scores above 7 in comparation with controls or less aggressive phenotypes of the tumor. Moreover, IGF2 expression is decreased in treatment resistant PC patients compared with sensitive ones.

Based on our previous results, we searched for the most interesting germline and somatic mutations in IGF2 related to PC. The combined effect of germline variants, which alter the structure, expression or function of protein-coding regions of cancer-biology related genes; determines which and how many somatic mutations must occur for malignant transformations; that is the reason why we also analyzed them (Qing et al., 2020). Concerning to somatic mutations, we found that rs758164144 is predominantly presented in patients with Gleason scores ≤7 contrasting with rs3842753 clustered in Gleason scores >7. This is the first time that rs758164144 and rs3842753 are described in PC.

Among germline variants, rs1004446 is the top one, according to bioinformatic analysis. This SNP has been previously associated with cancer, such as endometrial cancer risk (McGrath et al., 2011) and PC survival (Cao et al., 2014); or type 1 diabetes (McGrath et al., 2011). However, there are not many details in PC effect. For that reason, we have developed an analysis in blood and buccal swabs samples for testing the effect of this SNP according to aggressiveness or treatment response, but not statistical results were found.

Moreover, when conducting bioinformatic analysis in expression patterns, we discovered that NRP2 (Neuropilin-2) and KDR (Kinase Insert Domain Receptor) genes have the top positions for screening searching (G0 vs. G1 and G0 vs. G2). NRP2 is a member of the neuropilin receptor family and it is reported to regulate autophagy and mTORC2 signalling in PC. It has been identified as an important prognostic marker for worse clinical outcome especially in patients with high PC risk (Borkowetz et al., 2020). There is scarce data according to PC but in other tumors such as bladder cancer, high messenger RNA expression of NRP2, NRP1, PDGFC, and PDGFD are associated with a more aggressive disease (i.e., a high T stage, positive lymph node status and reduced survival) (Förster et al., 2021). Present data reports similarities as previously described in bladder cancer. Concerning KDR, there are not many published data KDR in PC. Just A. Fraga et al. demonstrated that KDR−604 T > C was correlated with protein level, accounting for a potential gene-environment effect in the activation of hypoxia-driven pathways in PC (Förster et al., 2021). In colorectal cancer, for example, a significant association was found between KDR expression, disease stage and lymph status (Fraga et al., 2017). Here we report higher expression patterns in PC in contrast to controls, as well as differential expression in treatment management.

Finally, we conducted a miRNA analysis, highlighting the role of miR-93-5p and 200c-3p. miR-200c-3p has previously been associated with PC aggressiveness, by its epithelial traits that leads to the anticipation of molecular reprogramming of Zeb1-Slug/vimentin axis (Basu et al., 2020). Recent data also indicated the role in PC progression of miR-200b-3p/200c-3p and XBP1 (X-box binding protein 1) as critical upstream regulators of PRKAR2B (type II-beta regulatory subunit of PKA) (Xia et al., 2020). Here we found that miR-200c-3p is situated in the top position according to bioinformatic analysis when comparing Gleason scores classification and PC absence. Thus, this suggest this miRNA as a good screening biomarker option. In relation to miR-93-5p, we have combined bioinformatic analyses with experimental ones, with promising results according to PC aggressiveness and non-invasive biomarkers. miR-93-5p has been previously reported in PC associated with lymphatic dissemination in locally advanced PC (Pudova et al., 2020), or combined with E2F2 (E2F transcription factor 2), RRM2 (ribonucleotide reductase regulatory subunit M2), and PKMYT1 (protein kinase, membrane associated tyrosine/threonine 1) genes and other three miRNAs (hsa-mir-17-5p, hsa-mir-20a-5p, hsa-mir-92a-3p), which marked this miRNA with promising therapeutic options in PC (Wei et al., 2020).

To sum up, here we report the role of IGF2 as an important marker for aggressiveness in PC. rs1004446 is, for the first time, included as a main somatic and germline mutation in this tumor. Although here, we just found statistically significance when comparing bioinformatic analysis, a deeper analysis with more samples will improve present data. Moreover, NRP2 and KDR have also been included as top screening biomarkers, according to bioinformatic analysis, which opens new strategies in the inclusion of these biomarkers in PC screening. Finally, we found that miR-93-5p could be an efficient strategy as an aggressiveness biomarker with non-invasive techniques.