Type 2 Diabetes-Related Variants Influence the Risk of Developing Prostate Cancer: A Population-Based Case-Control Study and Meta-Analysis

Simple Summary We investigated the influence of GWAS-identified variants for T2D in modulating prostate cancer (PCa) risk through a meta-analysis of our data with those from the UKBiobank and FinnGEn cohorts and four large European cohorts. We found that genetic variants within the FTO, HNF1B, and JAZF1 loci were associated with PCa risk. Our results also suggested, for the first time, a potentially interesting association of SNPs within NOTCH2 and RBMS1 genes that need to be further explored and validated. This study also shed some light onto the functional mechanisms behind the observed associations, and demonstrated that the HNF1Brs7501939 polymorphism correlated with lower levels of SULT1A1, an enzyme responsible for the sulfate conjugation of multiple endogenous and exogenous compounds. Furthermore, we found that SNPs within the HFN1B, NOTCH2, and RBMS1 genes impacted PCa risk through the modulation of mRNA gene expression levels of their respective genes. However, given the healthy nature of the subjects included in the cohort used for functional experiments, the link between the HNF1B locus and SULT1A1 should be considered still speculative and, therefore, requires further validation. Abstract In this study, we have evaluated whether 57 genome-wide association studies (GWAS)-identified common variants for type 2 diabetes (T2D) influence the risk of developing prostate cancer (PCa) in a population of 304 Caucasian PCa patients and 686 controls. The association of selected single nucleotide polymorphisms (SNPs) with the risk of PCa was validated through meta-analysis of our data with those from the UKBiobank and FinnGen cohorts, but also previously published genetic studies. We also evaluated whether T2D SNPs associated with PCa risk could influence host immune responses by analysing their correlation with absolute numbers of 91 blood-derived cell populations and circulating levels of 103 immunological proteins and 7 steroid hormones. We also investigated the correlation of the most interesting SNPs with cytokine levels after in vitro stimulation of whole blood, peripheral mononuclear cells (PBMCs), and monocyte-derived macrophages with LPS, PHA, Pam3Cys, and Staphylococcus Aureus. The meta-analysis of our data with those from six large cohorts confirmed that each copy of the FTOrs9939609A, HNF1Brs7501939T, HNF1Brs757210T, HNF1Brs4430796G, and JAZF1rs10486567A alleles significantly decreased risk of developing PCa (p = 3.70 × 10−5, p = 9.39 × 10−54, p = 5.04 × 10−54, p = 1.19 × 10−71, and p = 1.66 × 10−18, respectively). Although it was not statistically significant after correction for multiple testing, we also found that the NOTCH2rs10923931T and RBMS1rs7593730 SNPs associated with the risk of developing PCa (p = 8.49 × 10−4 and 0.004). Interestingly, we found that the protective effect attributed to the HFN1B locus could be mediated by the SULT1A1 protein (p = 0.00030), an arylsulfotransferase that catalyzes the sulfate conjugation of many hormones, neurotransmitters, drugs, and xenobiotic compounds. In addition to these results, eQTL analysis revealed that the HNF1Brs7501939, HNF1Brs757210, HNF1Brs4430796, NOTCH2rs10923931, and RBMS1rs7593730 SNPs influence the risk of PCa through the modulation of mRNA levels of their respective genes in whole blood and/or liver. These results confirm that functional TD2-related variants influence the risk of developing PCa, but also highlight the need of additional experiments to validate our functional results in a tumoral tissue context.


Introduction
Prostate cancer (PCa) is the second most common cancer worldwide and one of the first leading causes of cancer-related deaths in men in developed countries [1]. It accounts for 7.3% of all cancers, with an incidence of 37.5 per 100,000 individuals [2]. Despite the refinements in prevention strategies, a total of 1.4 million new cases were diagnosed in 2020 [2], and the incidence of the disease is increasing over the world, likely due to the interaction of both inherited and modifiable factors [3,4].
Although PCa has a high prevalence among males, only age, family history, and ethnicity have been established as major risk factors for the disease, with an attributable effect ranging from 5 to 9% of the cases [5]. In addition, rare highly penetrant mutations in specific genes, high levels of endogenous androgens, smoking, alcohol consumption, exposure to chemical compounds, sexually transmitted infections, diet, obesity, insulin-like growth factors, and type 2-diabetes (T2D) have been suggested as important modulators of the prostatic tumorigenesis [6]. Among the modifiable risk factors, T2D has attracted significant attention, since it has been consistently identified as a protective factor for PCa development [7][8][9], but it seems to induce disease progression. Several studies have suggested that the use of anti-diabetic drugs such as metformin might account for this protective effect of T2D on PCa risk [10,11], but other studies were not able to confirm these results in larger cohorts [12], which suggested that the protective effect attributed to T2D on PCa depend on common molecular pathways between these traits rather than the use of anti-diabetic drugs.
Although research on the specific pathways interfering in the development of T2D and PCa traits is still under way, a mounting body of evidence suggests that these diseases might share a genetic component [13][14][15][16]. Recent studies have reported that common genetic polymorphisms within T2D-related genes have an important role in modulating the risk of many cancers [17][18][19], but, so far, only a few studies have investigated the impact of diabetogenic variants on PCa risk, showing controversial results [20][21][22][23]. Considering this background, we decided to conduct a population-based case-control study including 994 subjects (304 PCa patients and 686 controls) to evaluate whether 57 diabetogenic variants identified through genome-wide association studies (GWAS) are associated with the risk of developing Pca. In order to validate the association of T2D-related markers with Pca risk, when possible, we performed a meta-analysis with data from previous genetic studies. Finally, we analyzed whether the most interesting markers correlated with absolute numbers of 91 blood-derived cell populations, 106 immunological serum proteins, 7 steroid hormones, and 9 cytokines (IFNγ, IL1β, IL1Ra, IL6, IL8, IL10, IL17, IL22, and TNFα) after stimulation of whole blood, peripheral blood mononuclear cells (PBMCs), and monocyte-derived macrophages with LPS, PHA, Pam3Cys, and Staphylococcus aureus.

Study Population
The study cohort consisted of 304 Caucasian PCa patients and 686 male healthy controls recruited in the Virgen de las Nieves University hospital (Granada, Spain) and the Complejo Hospitalario de Cáceres (Cáceres, Spain). Only patients without any prior history of malignancy, and who were not treated before blood withdrawal, were enrolled in this study. Patient characteristics are included in Table 1. The diagnosis of PCa was assigned by physicians, and fulfilled the international criteria [24]. Male controls with a mean age of 58.92 were blood donors from the Centro Regional de Transfusiones Sanguíneas de Granada (CRTS) and were selected from the same geographical region of the cases. In accordance with the Declaration of Helsinki, all participants gave their written informed consent to participate in the study and the ethical committees of the participant institutions approved the study.

SNP Selection and Genotyping
An extensive literature search concerning the mechanism of action of T2D-related genes was performed to select candidate genes that might affect the risk of developing PCa. SNPs were assessed on the basis of NCBI data, and were selected according to their known or putative functional consequences, i.e., their modifying influence on the structure of proteins, transcription level, or alternative splicing mechanisms. In total, 57 SNPs in 49 genes were selected for this study (Table 2).  Selected variants for T2D were genotyped using KASPar ® assays (LGC Genomics, London, UK) according to the manufacturer's instructions. For internal quality control, 5% of samples were randomly included as duplicates. Concordance between the original and the duplicate samples for the 57 SNPs was ≥99.0%. Call rates for all SNPs were ≥90.0%.

Statistical Analysis
The Hardy-Weinberg Equilibrium (HWE) tests were performed in the control group by a standard observed-expected chi-squared (χ 2 ) test. Logistic regression analyses were used to assess the effects of the genetic polymorphisms on PCa risk using dominant, recessive, and log-additive models of inheritance. Overall analyses were adjusted for age, and conducted using Stata (v12.1). Statistical power was calculated using the Quanto software (vs. 12.4; log-additive model).
In order to account for multiple testing, we calculated an adjusted significance level using data from the SNPclip Tool (https://ldlink.nci.nih.gov/?tab=snpclip, accessed on 8 May 2020), which consider the number of independent marker loci (n = 52). Given the high correlation between the log-additive and dominant models of inheritance, we corrected by log-additive and recessive models, resulting in a significant threshold for the main effect analysis of 0.00048 (0.05/52 SNPs/2 inheritance models). Since a study-wide significance threshold considering all these factors is generally perceived as a too conservative test, we also assessed the magnitude of observed associations between selected SNPs and risk of PCa through a quantile-quantile (QQ) plot generated from the results of the study population. The observed association p-values were ranked in order from smallest to largest on the y-axis and plotted against the expected results from a theoretical~χ 2 -distribution under the null hypothesis of no association on the x-axis.

Meta-Analysis
With the aim of assessing the consistency of the association between T2D-related SNPs and the risk of developing PCa, we performed a meta-analysis of our data with those from publicly available GWAS. We downloaded association estimates from the PheWeb site (https://pheweb.sph.umich.edu/, accessed on 11 May 2020) for 6311 PCa cases; 74,685 controls from the FinnGen research project; and 5993 PCa cases and 168,999 controls from the UK Biobank project (UKBiobank TOPMed-imputed). Details on genome-wide associations have been previously reported [70]. Briefly, analyses on binary outcomes were conducted using the SAIGE generalized mixed logistic regression model, adjusting for genetic relatedness, sex, birth year, and the first four principal components. For White British participants of the UK Biobank, endpoint definitions were generated from electronic health-records-derived ICD billing codes, and endpoint definitions for the FinnGen data can be found at risteys.finngen.fi (Risteys = intersection in Finnish). We also validated the association of genetic markers using data from previously published studies that were selected according to the following criteria: (1) GWAS or candidate-gene association studies found in PUBMED (https://www.ncbi.nlm.nih.gov/pubmed, accessed on 13 May 2020) using the following key words: prostate cancer, case-control association study, type 2 diabetes, genetic polymorphisms; (2) Studies using Caucasian populations; (3) Availability of association estimates according to a log-additive model of inheritance; (4) Hardy-Weinberg equilibrium in the control group; and (5) Written in English ( Figure 1). We pooled the Odds Ratios (ORs) using a fixed-effect model. Coefficients with a p-value ≤ 0.05 were considered significant. I 2 statistic was used to assess heterogeneity between studies. All statistics were calculated using STATA (v. 12).
the UK Biobank project (UKBiobank TOPMed-imputed). Details on genome-wide associations have been previously reported [70]. Briefly, analyses on binary outcomes were conducted using the SAIGE generalized mixed logistic regression model, adjusting for genetic relatedness, sex, birth year, and the first four principal components. For White British participants of the UK Biobank, endpoint definitions were generated from electronic health-records-derived ICD billing codes, and endpoint definitions for the FinnGen data can be found at risteys.finngen.fi (Risteys = intersection in Finnish). We also validated the association of genetic markers using data from previously published studies that were selected according to the following criteria: (1) GWAS or candidate-gene association studies found in PUBMED (https://www.ncbi.nlm.nih.gov/pubmed, accessed on 13 May 2020) using the following key words: prostate cancer, case-control association study, type 2 diabetes, genetic polymorphisms; (2) Studies using Caucasian populations; (3) Availability of association estimates according to a log-additive model of inheritance; (4) Hardy-Weinberg equilibrium in the control group; and (5) Written in English ( Figure  1). We pooled the Odds Ratios (ORs) using a fixed-effect model. Coefficients with a pvalue ≤ 0.05 were considered significant. I 2 statistic was used to assess heterogeneity between studies. All statistics were calculated using STATA (v. 12).

cQTL Analysis of the T2D-Related Variants
Cytokine stimulation experiments were conducted in the 500 Functional Genomics (500FG) cohort from the Human Functional Genomics Project (HFGP; http://www.humanfunctionalgenomics.org/, accessed on 7 July 2020). The HFGP study was approved by the Arnhem-Nijmegen Ethical Committee (no. 42561.091.12), and biological specimens were collected after informed consent was obtained. We investigated whether any of the 57 T2D-related SNPs correlated with cytokine levels (IFNγ, IL1β, IL1Ra, IL6, IL8, IL10, IL17, IL22, and TNFα) after the stimulation of peripheral blood mononuclear cells (PBMCs), macrophages, or whole blood from 172 healthy men with LPS (1 or 100 ng/mL), PHA (10 μg/mL), Pam3Cys (10 μg/mL), and Staphylococcus aureus. After log transformation, linear regression analyses adjusted for age were used to determine the correlation of selected SNPs with cytokine expression quantitative trait loci (cQTLs). All analyses were performed using R software (http://www.r-project.org/, accessed on 8 May 2020). In order to account for multiple comparisons, we used a significant threshold of 0.000106 (0.05/52 independent SNPs × 9 cytokines).

cQTL Analysis of the T2D-Related Variants
Cytokine stimulation experiments were conducted in the 500 Functional Genomics (500FG) cohort from the Human Functional Genomics Project (HFGP; http://www.humanfunctionalgenomics.org/, accessed on 7 July 2020). The HFGP study was approved by the Arnhem-Nijmegen Ethical Committee (no. 42561.091.12), and biological specimens were collected after informed consent was obtained. We investigated whether any of the 57 T2D-related SNPs correlated with cytokine levels (IFNγ, IL1β, IL1Ra, IL6, IL8, IL10, IL17, IL22, and TNFα) after the stimulation of peripheral blood mononuclear cells (PBMCs), macrophages, or whole blood from 172 healthy men with LPS (1 or 100 ng/mL), PHA (10 µg/mL), Pam3Cys (10 µg/mL), and Staphylococcus aureus. After log transformation, linear regression analyses adjusted for age were used to determine the correlation of selected SNPs with cytokine expression quantitative trait loci (cQTLs). All analyses were performed using R software (http://www.r-project.org/, accessed on 8 May 2020). In order to account for multiple comparisons, we used a significant threshold of 0.000106 (0.05/52 independent SNPs × 9 cytokines).
Details on PBMCs isolation, macrophage differentiation, and stimulation assays have been reported elsewhere [72][73][74]. Briefly, PBMCs were washed twice in saline, and suspended in medium (RPMI 1640) supplemented with gentamicin (10 mg/mL), L-glutamine (10 mM), and pyruvate (10 mM). PBMC stimulations were performed with 5 × 10 5 cells/well in round-bottom 96-well plates (Greiner) for 24 h in the presence of 10% human pool serum at 37 • C and 5% CO 2 . Supernatants were collected and stored in −20 • C until used for ELISA. LPS (100 ng/mL), PHA (10 µg/mL), and Pam3Cys (10 µg/mL) were used as stimulators for 24 or 48 h. Whole blood stimulation experiments were conducted using 100 µL of heparin blood that was added to a 48-well plate and subsequently stimulated with 400 µL of LPS and PHA (final volume 500 µL) for 48 h at 37 • C and 5% CO 2 . Supernatants were collected and stored in −20 • C until used for ELISA. Concentrations of human cytokines were determined using specific commercial ELISA kits (PeliKine Compact, Amsterdam, The Netherlands or R&D Systems, Minneapolis, MN, USA), following the manufacturer's instructions.

Correlation between T2D-Related Polymorphisms and Cell Counts of 91 Blood-Derived Immune Cell Populations and 103 Serum/Plasmatic Immunological Proteins
We also investigated whether selected polymorphisms had an impact on blood cell counts by analyzing a set of 91 manually annotated immune cell populations and genotype data from the 500 FG cohort that consisted of 172 healthy men (Supplementary Table S1). Cell populations were measured by 10-color flow cytometry (Navios flow cytometer, Beckman Coulter) after blood sampling (2-3 h), and cell count analysis was performed using the Kaluza software (Beckman Coulter, v. 1.3). In order to reduce inter-experimental noise and increase statistical power, cell count analysis was performed by calculating parental and grandparental percentages, which were defined as the percentage of a certain cell type within the cell populations one or two levels higher in the hierarchical definitions of cell sub-populations [75]. Detailed laboratory protocols for cell isolation, reagents, gating, and flow cytometry analysis have been reported elsewhere [76], and the accession number for the raw flow cytometry data and analyzed data files are available upon request to the authors (http://hfgp.bbmri.nl, accessed on 8 May 2020). A proteomic analysis was also performed in serum and plasma samples from the 500 FG cohort. Circulating proteins were measured using the commercially available Olink ® Inflammation panel (Olink, Sweden), which resulted in the measurement of 103 different biomarkers (Supplementary Materials  Table S2). Proteins levels were expressed on a log2-scale as normalized protein expression values, and normalized using bridging samples to correct for batch variation. Considering the number of proteins (n = 103) and cell populations (n = 91) tested, p-values of 9.33 × 10 −6 and 1.05 × 10 −5 were set as significant thresholds for the proteomic and cell-level variation analysis, respectively.

Correlation between Steroid Hormone Levels and T2D-Related SNPs
We also measured serum levels of seven steroid hormones (androstenedione, cortisol, 11-deoxy-cortisol, 17-hydroxy progesterone, progesterone, testosterone, and 25 hydroxy vitamin D3) in the 500 FG cohort. Complete protocol details of steroid hormone measurements have been reported elsewhere [74]. Hormone levels and genotyping data were available for a total of 167 subjects. After log-transform, correlation between steroid hormone levels and T2D-related SNPs was evaluated by linear regression analysis adjusted for age. In order to avoid a possible bias, we excluded from the analysis those subjects that were using oral contraceptives, or those subjects in which this information was not available. The significance threshold was set to 0.000137 considering the number of independent SNPs tested (n = 52) and the number of hormones determined (n = 7).

In Silico Functional Analysis
Once we assessed the correlation of T2D-related SNPs with cytokine and steroid hormone levels, we used the HaploReg SNP annotation tool to further investigate the functional consequences of each specific variant (http://www.broadinstitute.org/mammals/ haploreg/haploreg.php, accessed on 8 July 2020). We also assessed whether any of the potentially interesting markers correlated with mRNA expression levels of their respective genes using data from public eQTL browsers (GTex portal; www.gtexportal.org/home/, accessed on 8 July 2020; https://genenetwork.nl/bloodeqtlbrowser/, accessed on 8 July 2020) [77].

Overall Associations of Selected SNPs with PCa Risk
All SNPs were in HWE in the control group (p > 0.001). Logistic regression analysis adjusted for age showed that carriers of the IGF2BP2 rs4402960T/T , TCF7L2 rs12255372T/T , and TSPAN8|LGR5 rs7961581C/C genotypes had an increased risk of PCa (p = 0.037; 0.005 and 0.024), whereas those carrying the CDKAL1 rs7754840C , FLJ39370 rs17044137A , FTO rs9939609A , HNF1B rs7501939T , HNF1B rs757210T , JAZF1 rs10486567A , KCNQ1 rs2237897C , and KCNQ1 rs2237892C alleles showed a decreased risk of developing the disease (p = 0.022, 0.021, 0.046, 0.030, 0.024,  Table 3). Although none of the reported associations remained statistically significant after a stringent correction for multiple testing (p = 0.00048), the QQ plot showed a pronounced and early deviation of identity line, which confirmed that the effect attributed to SNPs in T2D-related loci was more than expected under the null hypothesis and, therefore, might represent true associations ( Figure 2). Table 3. Association of T2D-related variants and risk of developing PCa in the discovery population.

Overall Associations of Selected SNPs with PCa Risk
All SNPs were in HWE in the control group (p > 0.001). Logistic regression analysis adjusted for age showed that carriers of the IGF2BP2rs4402960T/T, TCF7L2rs12255372T/T, and TSPAN8|LGR5rs7961581C/C genotypes had an increased risk of PCa (p = 0.037; 0.005 and 0.024), whereas those carrying the CDKAL1rs7754840C, FLJ39370rs17044137A, FTOrs9939609A, HNF1Brs7501939T, HNF1Brs757210T, JAZF1rs10486567A, KCNQ1rs2237897C, and KCNQ1rs2237892C alleles showed a decreased risk of developing the disease (p = 0.022, 0.021, 0.046, 0.030, 0.024, 0.011, 0.041, and 0.0002; Table 3). Although none of the reported associations remained statistically significant after a stringent correction for multiple testing (p = 0.00048), the QQ plot showed a pronounced and early deviation of identity line, which confirmed that the effect attributed to SNPs in T2D-related loci was more than expected under the null hypothesis and, therefore, might represent true associations ( Figure 2).   The identity line represents the null hypothesis (no significant association between T2D-related SNPs and PCa risk). Early deviation of the identity line might represent true associations.

Meta-Analysis
In order to confirm these potentially interesting associations, we conducted a metaanalysis with GWAS data from two large cohorts (UKBiobank and FinnGen) and previously published genetic association studies. After filtering all studies found in the literature according to selected key words, we found that four case-control studies met the eligibility criteria [20][21][22]71]. The meta-analysis of our data with those from all these studies confirmed that carriers of the FTO rs9939609A , HNF1B rs7501939T , HNF1B rs757210T , HNF1B rs4430796G , and JAZF1 rs10486567A alleles had a decreased risk of developing PCa (p = 3.70 × 10 −5 , 9.39 × 10 −54 , 5.04 × 10 −54 , 1.19 × 10 −71 , and 1.66 × 10 −18 ; Table 4). Although the effect of the FTO and HNF1B loci on PCA risk has been consistently validated in previous studies, this is the first validation study confirming the association of the JAZF1 variant with the risk of developing the disease. In addition, although the association did not remain signifi-cant after correction for multiple testing, the meta-analysis suggested modest associations with the risk of developing the disease for SNPs within the NOTCH2 and RBMS1 loci (p = 8.49 × 10 −4 and 0.004; Table 4). These associations are potentially interesting and need to be further investigated.

Functional Characterization of T2D-Related Variants in the HFGP Cohort
In order to test the possible functional relevance of the most interesting SNPs, we analyzed data from the HFGP cohort. The proteomic analysis of the immunological serum proteins showed that carriers of the HNF1B rs7501939T allele had decreased circulating levels of ST1A1 protein (p = 0.00030; Figure 3). Although this correlation did not survive multiple testing correction, this result supported the implication of the HNF1B locus in modulating PCa risk, likely by the regulation of the sulfatation of multiple compounds in the liver. In addition, given the healthy nature of the subjects included in the HFGP cohort, this result is still speculative and needs to be further confirmed in tumor samples of PCa patients. No significant correlation was found between HFN1B, FTO, JAZF1, and NOTCH2 SNPs and blood-derived cells populations, steroid hormones, or cQTL data, which suggested that these SNPs might affect PCa risk likely through the regulation of mRNA levels of their respective genes. and blood-derived cells populations, steroid hormones, or cQTL data, which suggested that these SNPs might affect PCa risk likely through the regulation of mRNA levels of their respective genes. Next, we assessed functional information from the HaploReg SNP annotation tool, and we also assessed whether any of the potentially interesting markers correlated with mRNA expression levels of their respective genes using data from public eQTL browsers. We found that, according to Haploreg data, the HNF1Brs7501939, HNF1Brs757210, and HNF1Brs4430796 SNPs were modestly associated with mRNA HNF1B expression levels in peripheral whole blood (p = 9.23 × 10 −4 , 1.00 × 10 −3 , and 2.00 × 10 −3 ) [77], and that they mapped among histone marks (H3K4me1, H3K4me3) in several tissues and changed motifs for CEBPB, DMRT5, p300, HES1, and Maf (Supplementary Table S3). Similarly, we found that the NOTCH2rs10923931 and RBMS1rs7593730 SNPs also correlated with NOTCH2 and RBMS1 mRNA expression levels in peripheral whole blood (p = 7.3 × 10 −6 and 3.31 × 10 −7 , respectively) [77]. On the other hand, we found that the FTOrs9939609 and JAZF1rs10486567 SNPs mapped among histone marks in multiple tissues and several immune cells, and changed regulatory motifs for multiple regulatory transcription factors (Supplementary Table S3). In addition, GTEx portal data suggested that the NOTCH2rs10923931 SNP is an eQTL in the pancreas and liver (Supplementary Figure S1).

Discussion
T2D has been consistently identified as protective factor for PCa development and disease progression [7][8][9]. Several studies have also suggested that both diseases might share a genetic component [13][14][15][16], and some others have attempted to demonstrate the impact of diabetogenic variants on PCa risk, showing controversial results [20][21][22][23]. With Next, we assessed functional information from the HaploReg SNP annotation tool, and we also assessed whether any of the potentially interesting markers correlated with mRNA expression levels of their respective genes using data from public eQTL browsers. We found that, according to Haploreg data, the HNF1B rs7501939 , HNF1B rs757210 , and HNF1B rs4430796 SNPs were modestly associated with mRNA HNF1B expression levels in peripheral whole blood (p = 9.23 × 10 −4 , 1.00 × 10 −3 , and 2.00 × 10 −3 ) [77], and that they mapped among histone marks (H3K4me1, H3K4me3) in several tissues and changed motifs for CEBPB, DMRT5, p300, HES1, and Maf (Supplementary Table S3). Similarly, we found that the NOTCH2 rs10923931 and RBMS1 rs7593730 SNPs also correlated with NOTCH2 and RBMS1 mRNA expression levels in peripheral whole blood (p = 7.3 × 10 −6 and 3.31 × 10 −7 , respectively) [77]. On the other hand, we found that the FTO rs9939609 and JAZF1 rs10486567 SNPs mapped among histone marks in multiple tissues and several immune cells, and changed regulatory motifs for multiple regulatory transcription factors (Supplementary  Table S3). In addition, GTEx portal data suggested that the NOTCH2 rs10923931 SNP is an eQTL in the pancreas and liver.   Abbreviations: SNP, single nucleotide polymorphism; OR, odds ratio; CI, confidence interval; C, cytosine; T, thymine; A, adenine; G, guanosine; n/s, not specified. a Estimates calculated according to a log-additive model of inheritance and adjusted for age. Meta-analysis was performed assuming a fixed-effect model. p < 0.05 in bold.  Authors report the effect found for the rs4411878 (a SNP in complete linkage disequilibrium with the rs4607103, r 2 = 0.97).  Authors report the effect found for the rs11257655 (a SNP in strong linkage disequilibrium with the rs12779790, r 2 = 0.82).  Authors report the effect found for the rs7756992 (a SNP in strong linkage disequilibrium with the rs7754840, r 2 = 0.75).  Authors report the effect found for the rs8050136 (a SNP in complete linkage disequilibrium with the rs9969309, r 2 = 0.98). Authors report the effect found for the rs1635852 (a SNP in complete linkage disequilibrium with the rs864745, r 2 = 0.98). * Authors report the effect found for the rs2641348 (a SNP in complete linkage disequilibrium with the rs10923931, r 2 = 1.00).  Authors report the effect found for the rs4607517 (a SNP in complete linkage disequilibrium with the rs1799884, r 2 = 1.00).  Authors report the effect found for the rs780094 (a SNP in complete linkage disequilibrium with the rs1260326, r 2 = 0.92). Ϯ Authors report the effect found for the rs13414140 (a SNP in complete linkage disequilibrium with the rs7578597, r 2 = 1.00).  Authors report the effect found for the rs1353362 (a SNP in strong linkage disequilibrium with the rs7961581, r 2 = 0.92).  Authors report the effect found for the rs10830963 (a SNP in moderate to high linkage disequilibrium with the rs1387153, r 2 = 0.67

Functional Characterization of T2D-Related Variants in the HFGP Cohort
In order to test the possible functional relevance of the most interesting SNPs, we analyzed data from the HFGP cohort. The proteomic analysis of the immunological serum proteins showed that carriers of the HNF1Brs7501939T allele had decreased circulating levels of ST1A1 protein (p = 0.00030; Figure 3). Although this correlation did not survive multiple testing correction, this result supported the implication of the HNF1B locus in modulating PCa risk, likely by the regulation of the sulfatation of multiple compounds in the liver. In addition, given the healthy nature of the subjects included in the HFGP cohort, this result is still speculative and needs to be further confirmed in tumor samples of PCa patients. No significant correlation was found between HFN1B, FTO, JAZF1, and NOTCH2 SNPs Abbreviations: SNP, single nucleotide polymorphism; OR, odds ratio; CI, confidence interval; C, cytosine; T, thymine; A, adenine; G, guanosine; n/s, not specified. a Estimates calculated according to a log-additive model of inheritance and adjusted for age. Meta-analysis was performed assuming a fixed-effect model. p < 0.05 in bold. η Authors report the effect found for the rs4411878 (a SNP in complete linkage disequilibrium with the rs4607103, r 2 = 0.97). ξ Authors report the effect found for the rs11257655 (a SNP in strong linkage disequilibrium with the rs12779790, r 2 = 0.82). # Authors report the effect found for the rs7756992 (a SNP in strong linkage disequilibrium with the rs7754840, r 2 = 0.75). δ Authors report the effect found for the rs8050136 (a SNP in complete linkage disequilibrium with the rs9969309, r 2 = 0.98). ℵ Authors report the effect found for the rs1635852 (a SNP in complete linkage disequilibrium with the rs864745, r 2 = 0.98). * Authors report the effect found for the rs2641348 (a SNP in complete linkage disequilibrium with the rs10923931, r 2 = 1.00). ∂ Authors report the effect found for the rs4607517 (a SNP in complete linkage disequilibrium with the rs1799884, r 2 = 1.00). ∏ Authors report the effect found for the rs780094 (a SNP in complete linkage disequilibrium with the rs1260326, r 2 = 0.92). Abbreviations: SNP, single nucleotide polymorphism; OR, odds ratio; CI, confidence interval; C, cytosine; T, thymine; A, adenine; G, guanosine; n/s, not specified. a Estimates calculated according to a log-additive model of inheritance and adjusted for age. Meta-analysis was performed assuming a fixed-effect model. p < 0.05 in bold.  Authors report the effect found for the rs4411878 (a SNP in complete linkage disequilibrium with the rs4607103, r 2 = 0.97).  Authors report the effect found for the rs11257655 (a SNP in strong linkage disequilibrium with the rs12779790, r 2 = 0.82).  Authors report the effect found for the rs7756992 (a SNP in strong linkage disequilibrium with the rs7754840, r 2 = 0.75).  Authors report the effect found for the rs8050136 (a SNP in complete linkage disequilibrium with the rs9969309, r 2 = 0.98). Authors report the effect found for the rs1635852 (a SNP in complete linkage disequilibrium with the rs864745, r 2 = 0.98). * Authors report the effect found for the rs2641348 (a SNP in complete linkage disequilibrium with the rs10923931, r 2 = 1.00).  Authors report the effect found for the rs4607517 (a SNP in complete linkage disequilibrium with the rs1799884, r 2 = 1.00).  Authors report the effect found for the rs780094 (a SNP in complete linkage disequilibrium with the rs1260326, r 2 = 0.92). Ϯ Authors report the effect found for the rs13414140 (a SNP in complete linkage disequilibrium with the rs7578597, r 2 = 1.00).  Authors report the effect found for the rs1353362 (a SNP in strong linkage disequilibrium with the rs7961581, r 2 = 0.92).  Authors report the effect found for the rs10830963 (a SNP in moderate to high linkage disequilibrium with the rs1387153, r 2 = 0.67

Functional Characterization of T2D-Related Variants in the HFGP Cohort
In order to test the possible functional relevance of the most interesting SNPs, we analyzed data from the HFGP cohort. The proteomic analysis of the immunological serum proteins showed that carriers of the HNF1Brs7501939T allele had decreased circulating levels of ST1A1 protein (p = 0.00030; Figure 3). Although this correlation did not survive multiple testing correction, this result supported the implication of the HNF1B locus in modulating PCa risk, likely by the regulation of the sulfatation of multiple compounds in the liver. In addition, given the healthy nature of the subjects included in the HFGP cohort, this result is still speculative and needs to be further confirmed in tumor samples of PCa patients.
Authors report the effect found for the rs13414140 (a SNP in complete linkage disequilibrium with the rs7578597, r 2 = 1.00). τ Authors report the effect found for the rs1353362 (a SNP in strong linkage disequilibrium with the rs7961581, r 2 = 0.92). ς Authors report the effect found for the rs10830963 (a SNP in moderate to high linkage disequilibrium with the rs1387153, r 2 = 0.67). ς Results from Lewis et al. Plos One 2010; 5: e13485 [78].

Discussion
T2D has been consistently identified as protective factor for PCa development and disease progression [7][8][9]. Several studies have also suggested that both diseases might share a genetic component [13][14][15][16], and some others have attempted to demonstrate the impact of diabetogenic variants on PCa risk, showing controversial results [20][21][22][23]. With this background, we decided to further investigate the association of diabetogenic variants identified through GWAS with the risk to PCa, and attempted to identify the biological mechanisms underlying the most interesting associations through the analysis of functional data from the HFGP cohort and eQTL browsers.
The meta-analysis of the Spanish cohort with those from the UKBiobank, FinnGen, and previously published studies [20][21][22] confirmed that carriers of the FTO rs9939609A , HNF1B rs7501939T , HNF1B rs757210T , HNF1B rs4430796G , and JAZF1 rs10486567A alleles showed a decreased risk of developing the disease. Although the association did not reach the stringent significance threshold, we also found that the NOTCH2 rs10923931 and RBMS1 rs7593730 SNPs associated with the risk of developing the disease. The strongest effect on PCa risk was observed for SNPs within the HNF1B locus (rs757210, rs7501939, and rs4430796), which showed a similar direction across all study populations. The HNF1B (TCF2) gene is located at chromosome 17q12, and it encodes for a transcription factor implicated in the control of regulatory networks related to pancreas and kidney development. It has been reported that the HNF1B locus participates not only in the generation of endocrine precursors, but also in the modulation of acinar cell identity and duct morphogenesis. In addition to these functions, it has been consistently reported that the HNF1B locus plays a key role in modulating tumorigenesis in solid [79] and hematological cancers [80], and that its methylation or mRNA expression levels can be used for patient stratification [81] and prediction of disease outcome [82]. The association of the HNF1B rs7501939 , HNF1B rs757210 , and HNF1B rs4430796 polymorphisms with PCa risk was in agreement with results recently reported in GWAS for PCa [14,83,84], whereas large-scale fine mapping studies have even found additional polymorphisms that might contribute to the development of PCa [71]. These findings, together with our functional results reporting that carriers of the HNF1B rs7501939T and HNF1B rs757210T alleles showed decreased levels of SULT1A1 protein (also known as ST1A1), suggest that the effect of the HFN1B locus on PCa risk might be mediated through the regulation of SULT1A1 expression levels. This protein is an enzyme that catalyzes the sulfate conjugation of many hormones, neurotransmitters, drugs, and xenobiotic compounds, among other compounds. It also has been demonstrated that SULT1A1 regulates the metabolic activation of carcinogenic N-hydroxyarylamines, leading to highly reactive intermediates capable of forming DNA adducts, which could result in mutagenesis [85]. In support of the hypothesis of a tumorigenic effect of the SULT1A1, several studies have shown that an increased expression of the SULT1A1 mRNA expression levels contributes to PCa development [86,87]. Although our functional experiments were conducted in a cohort of healthy donors and, therefore, cannot be directly translated to a disease context, our experimental data are in agreement with previous studies, and suggest that the protective effect attributed to the HNF1B rs7501939 and HNF1B rs757210 SNPs could be mediated by a reduction in the expression of the SULT1A1 protein. Furthermore, it has been reported that the HNF1B rs7501939 , HNF1B rs757210 , and HNF1B rs4430796 SNPs are modestly associated with mRNA HNF1B expression levels in peripheral blood, which might help to explain how these genetic variants may influence the risk of developing PCa [77]. However, despite these interesting data, we think that the biological link between the HNF1B locus and SULT1A1 is still speculative and needs to be further explored and validated, since, if confirmed, it might represent a potentially interesting therapeutic target. An option to confirm this hypothesis would be to measure SULT1A1 levels in tumoral tissues.
Besides these results, this study also confirmed the association of the FTO locus with PCa risk. The FTO gene is located on chromosome 16q12.2, and it has been implicated in determining not only obesity, but also other symptoms of the metabolic syndrome. In addition, it has been reported that the FTO gene acts as a tumor suppressor gene by regulating the proliferation, migration, and invasion of PCa cells, and the FTO expression level had a relevance with the development of PCa and the prognosis of PCa patients [88]. Although the association of the FTO rs9939609 polymorphism with PCa risk was weak in all previous studies, and might depend on different confounding factors, the meta-analysis performed in this study confirmed a strong and consistent association of this intronic variant with a decreased risk of PCa. A recent study also demonstrated that the association of the FTO rs9939609 SNP with a decreased risk of PCa was found in non-European populations, and that the presence of this genetic marker tended to be associated with disease severity in patients that were overweighted [89]. These findings, together with the lack of significant results in our functional studies, suggest that the role of the FTO locus in determining PCa might be mediated by complex obesogenic and/or diabetogenic mechanisms. In support of this hypothesis, we found that the association of the FTO rs9939609 SNP with a decreased risk of PCa showed a similar direction to the one observed in the GWAS for T2D.
Similarly, we also found that the presence of the JAZF1 rs10486567 SNP was inversely associated with the risk of developing PCa. These results were again in concordance with previous GWAS that have consistently reported that JAZF1 is a susceptibility locus for PCa [84,[90][91][92][93]. The JAZF1 gene is located at 7p15, and it encodes for a zinc finger protein that is overexpressed in the human prostate tissue where it induces cell proliferation, migration, and invasion, and tumor development [94]. Even though there is not much evidence about the functional role of the JAZF1 rs10486567 SNP in PCa, it has been demonstrated that the deletion of the JAZF1 locus is associated with reduced levels of IGF-1 and insulin resistance in mice [95], which suggests that the presence of functional polymorphisms within this locus might act to promote PCa development through diabetogenic mechanisms.
Finally, although the association was not statistically significant after correction for multiple testing, it seems to be reasonable to suggest that genetic variants within the NOTCH2 and RBMS1 genes could weakly influence the risk of developing PCa. In this regard, it has been reported that genes of the NOTCH family play a relevant role in multiple cancers, including PCa, and that their deregulation may be a key event in tumor onset and disease progression [96][97][98][99][100][101][102]. The human NOTCH2 locus is located in the chromosomal region 1p11.2, and it plays an important role in modulating prostate development and homeostasis, and its deregulation induces proliferation and expansion of both basal and luminal cells in the prostate [97]. In addition, it has been reported that NOTCH activity promotes prostate cancer cell migration [97], invasion [96,97], aggressiveness [98], and metastasis [99], and that its silencing induces apoptosis and increases the chemosensitivity of PCa cells [100]. However, a tumor suppressive role of the NOTCH pathway has also been suggested in the literature [103,104], which points towards the need of additional studies to elucidate the role of these genes in the etiopathogenesis of PCa. On the other hand, the meta-analysis of all study cohorts suggested that carriers of the RBMS1 rs7593730T allele had an increased risk of developing PCa. The RBMS1 gene is located at chromosome 2q24.2, and it encodes for a protein that binds single stranded DNA/RNA and plays an important role in DNA replication, cell cycle progression, gene transcription, and apoptosis. A recent study demonstrated that the RBMS1 locus acts by regulating the expression of the miR-106b [105], which has been found overexpressed in hepatocellular carcinoma [106], cervical cancer [107], renal carcinoma [108], and gastric cancer [109]. At the same time, Dankert and collaborators have also demostrated that miR-106b can regulate endogenous RBMS1 expression in PCa cell lines and, thereby, act as a tumor suppressor gene with inhibitory effects on colony formation and cell growth [105]. Despite the lack of evidence linking RBMS1 SNPs and risk of PCa, it seems to be plausible to suggest that the presence of genetic variants in the RBMS1 locus might control miR-106b levels and, therefore, favors tumorigenesis. In support of this hypothesis, haploreg data showed that the RBMS1 rs7593730 SNP is associated with different mRNA RBMS1 expression levels in several tissues and cells [77], and that it maps among histone marks in multiple tissues and several immune cells, and changed regulatory motifs for multiple regulatory transcription factors. Nonetheless, although interesting, the effect of NOTCH2 and RBMS1 SNPs on PCa risk must be considered as preliminary and, therefore, needs to be further confirmed in independent cohorts. This study has both strengths and weaknesses. The major strength of this study is the large number of genetic markers analyzed that allowed us to perform a well-powered meta-analysis of our data with those from previous studies. The meta-analysis of all study cohorts allowed us to not only confirm previous associations between T2D-related polymorphisms and PCa risk, but also to identify potentially interesting genetic markers for the disease. Although the discovery population was relatively small and the influence of diabetogenic variants on the risk of the disease was expected to be modest, our study was sufficiently powered to detect such small effects. Based on the genotype frequencies observed in our study cohort, we had 80% of power (log-additive model) to detect an odds ratio of 1.59 at alpha = 0.00048 (multiple testing threshold) for a polymorphism with a minor allele frequency of 0.25. Likewise, we comprehensively analyzed the impact of T2D-related SNPs in modulating immune responses, blood cell counts, steroid hormones, and serum and plasma metabolites in a relatively large cohort of healthy subjects. However, we also have important limitations. One of them was the fact that functional characterization of the most interesting SNPs was conducted in a healthy control cohort rather than in PCa patients. In addition, we did not have access to medication history, T2D status, and BMI for a substantial number of PCa cases included in the meta-analyses, which did not allow us to adjust our analyses for these confounding variables and, consequently, to rule out the possibility that some of the reported associations could arise as a result of a different distributions of diabetic and/or obese subjects between PCa cases and controls, or because of the effect of diabetes medication rather than the diabetes condition itself. Nonetheless, previous studies have reported that the effect of T2D-related variants on the risk of PCa was independent of T2D status and BMI.

Conclusions
In conclusion, our study indicates that T2D-related variants within the HNF1B, FTO, and JAZF1 genes influence the risk of PCa likely through the modulation of diabetogenic pathways, and suggests, for the first time, an association of SNPs within the NOTCH2 and RBMS1 loci that need to be validated in independent cohorts. This study also suggests that the effect of the HFN1B SNPs on PCa risk might be mediated by not only the ST1A1 protein, but also HFN1B mRNA expression levels, whereas the effect of the FTO, JAZF1, NOTCH2, and RBMS1 SNPs on PCa risk seem to be involved in the regulation of mRNA expression levels of their respective genes.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cancers14102376/s1. Table S1: Cell types analyzed either in whole blood or peripheral mononuclear blood cells; Table S2: Serum and plasma metabolites measured in the HFGP cohort; Table S3: Summary of functional and regulatory annotation of the most significant SNPs. Funding: This study was supported by grants from the FIBAO foundation (Granada, Spain) and from the Instituto de Salud Carlos III (PI12/02688, PI17/02256 and PI20/01845; Madrid, Spain).

Institutional Review Board Statement:
The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Institutional Review Boards of all institutions participating in patient recruitment (Virgen de las Nieves University Hospital, 18012, Granada, Spain; Hospital de San Pedro Alcántara, 10003, Cáceres, Spain).
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The genotype data used in the present study are available from the corresponding authors upon reasonable request. Functional data used in this project have been meticulously catalogued and archived in the BBMRI-NL data infrastructure (https://hfgp.bbmri.nl/, accessed on 7 July 2020) using the MOLGENIS open-source platform for scientific data [52]. This allows flexible data querying and download, including sufficiently rich metadata and interfaces for machine processing (R statistics, REST API), and using FAIR principles to optimize findability, accessibility, interoperability, and reusability [53,54].