A Multigene-Panel Study Identifies Single Nucleotide Polymorphisms Associated with Prostate Cancer Risk

The immune system plays a critical role in modulating cancer development and progression. Polymorphisms in key genes involved in immune responses are known to affect susceptibility to cancer. Here, we analyzed 35 genes to evaluate the association between variants of genes involved in immune responses and prostate cancer risk. Thirty-five genes were analyzed in 47 patients with prostate cancer and 43 healthy controls using next-generation sequencing. Allelic and genotype frequencies were calculated in both cohorts, and a generalized linear mixed model was applied to test the relationship between prostate cancer risk and nucleotide substitution. Odds ratios were calculated to describe the association between each single nucleotide polymorphism (SNP) and prostate cancer risk. Significant changes in allelic and genotypic distributions were observed for IL4R, IL12RB1, IL12RB2, IL6, TMPRSS2, and ACE2. Furthermore, a generalized linear mixed model identified statistically significant associations between prostate cancer risk and SNPs in IL12RB2, IL13, IL17A, IL4R, MAPT, and TFNRS1B. Finally, a statistically significant association was observed between IL2RA and TNFRSF1B and Gleason scores, and between SLC11A1, TNFRSF1B and PSA values. We identified SNPs in inflammation and two prostate cancer-associated genes. Our results provide new insights into the immunogenetic landscape of prostate cancer and the impact that SNPs on immune genes may have on affecting the susceptibility to prostate cancer.


Introduction
Prostate cancer (PCa) is one of the most widespread noncutaneous cancers in men. The incidence of PCa appears to be affected by several risk factors, among which inflammation plays a prominent role [1,2]. Inflammation appears to play a special role in promoting the onset and progression of tumors by enabling an immunosuppressive tumor environment [3]. Inflammation is the physiological response of the body to deleterious stimuli such as pathogens and tissue damage, which results in tissue repair. In some conditions, dysregulation of this finely regulated process can result in a chronic inflammatory state, which, in turn, can sustain tumorigenesis. Cancer cells originate from the accumulation of genetic alterations, such as single nucleotide polymorphisms (SNPs), copy number variation, chromosomal rearrangements, loss of heterozygosity, and allelic imbalance, which result in dysregulation of cellular homeostasis and consequent expansion. Furthermore, cancer cells develop several systems to evade immune surveillance, such as loss of antigenicity, programmed death-ligand 1, and limited access to the primary tumor through the extracellular matrix. Furthermore, the immunological state of the tumor microenvironment (TME) is crucial for the implementation of cancer immunotherapy [4].

Allele and Genotype Frequency Distribution in PCa Patients and HCs
The present study analyzed a total of 47 patients with PCa and 43 HCs using NGS, and comprehensive gene profiling was conducted on a panel of 35 genes. Table 1 summarizes the distribution of significant allele and genotype frequencies between PCa patients and HCs. Among the nucleotide substitutions analyzed, seven SNPs, located within IL4R, IL12RB1, IL12RB2, IL6, TMPRSS2, and ACE2 were found associated with PCa risk. Both the C allele and the CC genotype of IL4R rs2301807 were inversely associated with PCa risk (OR: 0.007, CI: 0.000-0.118, p < 0.0001; OR: 0.014, CI: 0.001-0.244, p < 0.0001, respectively). The AA genotype in the 3 UTR of IL12RB1 rs3746190, along with the C allele of IL12RB2 rs2307145, were found associated with a lower risk of PCa (OR: 0.062, CI: 0.003-1.243, p = 0.013; OR: 0.194, CI: 0.039-0.971, p = 0.038, respectively). The G allele of the intronic variant IL6 rs2069832 was positively associated with a higher risk of PCa (OR: 2.583, CI: 1.092-6.114, p = 0.028). In addition, both the A allele and CA heterozygous genotype of the intronic variant TMPRSS2 rs140141551 were both associated with a higher risk of PCa (OR: 7.535, CI: 0.922-61.586, p = 0.0382; OR: 8.205, CI: 0.922-61.586, p = 0.024, respectively). Conversely, the G allele and the AG genotype of TMPRSS2 rs11701576 were found to be associated with a decreased risk of PCa (OR: 0.309, CI: 0.097-0.977, p = 0.037; OR: 0.263, CI: 0.077-0.899, p = 0.027, respectively). Finally, a synonymous variant was identified within the ACE2 gene, and the T allele of rs35803318 was found to be associated with a higher risk of PCa (OR: 11.743, CI: 0.659-209.216, p = 0.022).

Association between Nucleotide Substitution for Each Locus Examined and PCa Estimation
Multivariate regression analysis was performed using backward selection to assess the relationship between nucleotide substitutions at a given locus and PCa. To test this relationship, a generalized linear mixed model (GLMM) was implemented. The first GLMM runs were applied individually to each gene to identify loci statistically associated with the occurrence of PCa. Statistical significance was set at p < 0.05, with a few exceptions for borderline values, which were retained to avoid premature selection in intermediate models ( Table 2). Odds ratio (ORs) calculations were performed as well on statistically significant variables representing useful predictors (Table 3). Statistically significant variables identified were then used to build a new dataset, and GLMMs runs were applied.    Using the predictors identified in Table 3, the analysis was performed considering statistically significant variables together (Table 4). ORs were calculated to assess their association with PCa risk (Table 5) ACE, ATG7, CD28, IFNG, IFNGR1, IFNGR2, IL12A, IL12B, IL18, IL18R1, IL2,  IL2RA, IL4R, IL5, IL7R, IRF5, LAG3, R1I2, PACGR, PTPN22, TNFRS, TRAF6, and VDR.  Finally, GLMM analysis was performed to explore the relationship between statistically relevant variants identified as useful predictors in the initial model (Table 2) and the clinical characteristics of patients with PCa (Table 5). Particularly, the following clinical characteristics were used for the analysis: Gleason score (GS), serum prostate-specific antigen (PSA) levels, age, clinical stage (cTNM), and body mass index (BMI).
Two polymorphisms (IL2RA rs7093069 and TNFRSF1B rs2275416) were significantly associated with GS. Accordingly, the variables rs7093069 and rs2275416 (within IL2RA and TNFRSF1B genes, respectively) represent good predictors for the PCa severity. Indeed, with GS with values of around 6, rs7093069 and rs2275416 do not appear as informative variables, while with values of around 9, they appear strongly associated with the severity of the pathology. In addition, SLC11A1 rs17221959 and both TNFRSF1B rs636964 and rs2275416 were found to be significantly associated with serum PSA levels. Overall, only GS and PSA have been positively selected by the model, while other variables (age, cTNM, and BMI) were removed during the backward selection. Indeed, in the dataset here, analyzed age, cTNM, and BMI demonstrated to have no association with PCa, nor alone or taken together.

Discussion
The importance of the immune system in cancer onset and progression has been well established in cancer research. In the last decade, the notion of engaging the immune system as a patient-personalized weapon against cancer cells has undoubtedly revolutionized cancer therapy. A growing body of evidence has investigated the presence of variants in genes implicated in anti-tumor immune responses and cancer risk. Therefore, we investigated the presence of polymorphisms in a panel of 35 genes involved in pro-and anti-inflammatory responses, autophagy-lysosome pathway and two PCa-associated genes (ACE2 and TMPRSS2) in a cohort of patients with PCa and a cohort of HCs. We identified several novel polymorphisms that contributed to PCa susceptibility. In particular, we found that six SNPs in several key genes, namely IL4R, IL12RB1, IL12RB2, TMPRSS2, and TNFRSF1B, were associated with a lower risk of PCa. In contrast, the eight SNPs identified in IL6, TMPRSS2, ACE2, IL12RB2, IL13, IL17A, MAPT, and TNFRSF1B were associated with a higher risk of PCa. PCa patients with specific polymorphisms in T-helper cytokine genes (IL13 rs20541 and IL4R rs2301807) may have an altered T-helper response and consequently, a potential imbalance in their immune system. Altered expression and increased circulating levels of both pro-and anti-inflammatory cytokines have been reported in several types of cancer [11]. Particularly, the role of interleukin 4 (IL-4) and its receptor (IL-4R) in facilitating a pro-metastatic phenotype has been described in epithelial cancer cells [12][13][14]. IL-4R appears as a plausible candidate to attenuate metastatic tumor growth, especially since previous evidence demonstrated how antibody-mediated IL-4R neutralization lessened metastatic lung tumor burden and IL-4Rα loss reduced metastatic tumor growth [13]. Furthermore, IL-4 is able to activate the androgen receptor in conditions of low androgen levels [15]. Our finding that IL6 rs2069832 was positively associated with a higher risk of PCa is consistent with that of Slattery et al. [16], in which rs2069832 was associated with an increased risk of breast cancer as well. Interleukin 6 (IL-6) is a pro-inflammatory cytokine expressed in PCa tissues, and its increased serum levels were observed in patients with metastatic and castration-resistant PCa (CRPC) [17,18]. Moreover, an association between IL-6 serum levels and unfavorable prognosis has also been observed in many cancers other than PCa [18,19]. IL12RB1 rs3746190, IL12RB2 rs2307145, and rs2229546 were associated with decreased PCa risk. Interleukin 12 (IL-12) plays a pivotal role in stimulating natural killer (NK) cell cytotoxicity, inflammation, and boosting T-helper type I cells [20]. No other studies have associated both polymorphisms with PCa, although Kundu et al. demonstrated higher serum levels of IL-12 p40 monomer in patients with PCa than in healthy individuals [21]. In particular, the IL-12 p40 monomer helped cancer cells escape cell death by suppressing IL12RB1 internalization; therefore, PCa cell apoptosis was induced after IL-12 p40 neutralization and consequent IL12RB1 internalization. With regard to IL12RB2 rs2229546, our results are in line with another study that found an association between rs2229546 and a reduced risk of cervical squamous cell carcinoma and HPV-18-positive cervical cancer [9]. IL17A rs7747909 belongs to the 3 -UTR region of interleukin 17A (IL-17A) and was found associated with a higher risk of PCa. Interestingly, Bedoui et al. reported a strong association between IL17A polymorphisms, including rs7747909, and the development of tolerance to colon rectal cancer therapy [22]. IL-17A is a pro-inflammatory cytokine primarily produced by Th17 cells whose elevated expression has been associated with an excessive inflammatory state and tumor growth in prostate and colon cancers [23][24][25][26]. Interestingly, gene amplification of IL17 has been reported more frequently in CRPC and neuroendocrine prostate cancer than in hormone-naïve PCa [27]. IL13 rs20541 polymorphism, located in exon 4, appears to be strongly associated with increased serum IgE levels in children [28]. Interleukin 13 (IL-13) is primarily produced by activated Th2 cells and exhibits strong antitumor activity in PCa cell lines [29]. Moreover, rs20541 has previously been associated with cancer risk, but this association remains uncertain. In fact, in our study, IL13 rs20541 was positively associated with the risk of PCa, whereas it was previously reported to be associated with a decreased risk of glioma [30,31]. In addition, TMPRSS2 rs140141551 and ACE2 rs35803318 were associated with a higher risk of PCa than TMPRSS2 rs11701576, which was negatively associated with PCa risk. TMPRSS2 encodes a transmembrane serine protease 2 (TMPRSS2), which is expressed by prostate epithelial cells. To date, TMPRSS2 physiological function remains unclear, albeit its expression appears to be regulated by androgens in PCa cells [32,33]. In contrast, ACE2 encodes the angiotensin-converting enzyme (ACE2), which is aberrantly expressed in several cancers [34,35]. Compelling evidence demonstrated the strong association between ACE2 expression and tumor progression, prognosis, antitumor immunity, and immunotherapy response [36][37][38][39][40]. In fact, increased expression of ACE2 has been linked to better survival, probably due to the antitumor effects on angiogenesis and heightened immune infiltration in different cancer cohorts. Moreover, a positive correlation was also found between ACE2 expression and PD-L1, a predictive biomarker of response to immune checkpoint inhibitors. This correlation, in turn, was linked to a favorable response to immunotherapy.
ACE2 rs35803318 is a silent polymorphism which may affect the expression of ACE2 increasing the susceptibility to PCa, especially taking into account recent studies demonstrating the effect of silent mutations on gene function and expression, splicing mechanisms, and phenotypes by affecting miRNA binding, protein folding, and mRNA stability [41][42][43][44]. Finally, the analysis displayed a significant association between MAPT rs2258689 and increased risk of PCa. MAPT encodes microtubule-associated protein tau (MAPT) whose expression is not restricted to neurons. In truth, MAPT overexpression was previously associated with lower GS, suggesting its potential as marker of PCa prognosis [45] and response to docetaxel [46]. In our analysis, we also considered whether the SNPs identified by our model were associated with the clinical information of patients with PCa, particularly the GS and PSA serum levels. GS is a system implemented to help urologists evaluate staging and predict prognosis in patients with PCa. We found that this parameter was affected by IL2RA rs7093069 and TNFRSF1B rs2275416, implying the strong impact of genetic variability within immune genes on cancer development, progression, and staging. Alterations in IL2RA relative expression have been reported in groups of patients with different GS [47]. On the other hand, TNFRSF1B rs2275416 was associated with serum PSA levels, GS, and PCa risk. TNFRSF1B encodes the receptor for tumor necrosis factor-α (TNFR2), which plays a crucial role in promoting tumor growth, invasion, and metastasis. TNFR2 is highly expressed in T regulatory cells (Tregs) and its presence is crucial in the TME [48]. Hence, specific TNFR2-targeted cancer therapies and TNFR2-antagonist antibodies have been tested to regulate tumor growth, showing promising results in vitro and in vivo [49][50][51]. As the presence of TNFR2 in the blood and TME has been associated with unfavorable prognosis in several types of cancer [52], it is possible that the presence of TNFRSF1B rs2275416 may also affect PCa prognosis, given its association with both GS and PSA levels. SLC11A1 rs17221959 was also associated with serum PSA levels. This finding is also consistent with the results recently reported by Zhu et al. in which SLC11A1 rs7573065 C > T was associated with increased risk of PCa and low overall survival; meanwhile, higher mRNA expression was found in PCa tissues compared to normal tissues [53]. Despite the small sample size, our study is the first to investigate these genes in PCa. For this reason, future analysis, performed on a larger scale, is necessary to further confirm our findings and explore the prognostic impact of these variants in patients diagnosed with PCa.

Study Population, Clinical Characteristics, and Blood Samples Collection
The ethical committee of the AOU, Sassari, approved this study. Two cohorts, comprising 47 patients with PCa and 43 HCs, were recruited. The demographic and clinical characteristics of the patients with PCa and HCs are summarized in Table 6. Patients newly diagnosed with PCa were recruited by the Urology Department of the University Hospital of Sassari, whereas HCs were recruited by the Transfusion Center (AOU, Sassari, Italy). In PCa patients, blood sample collection was performed at the time of the diagnosis. Peripheral blood samples were collected from both cohorts using K + -EDTA test tubes. PCa assessment was based on the clinical examination and histopathological analysis of isolated biopsies performed in all patients. Table 6. Clinical and demographic information of the PCa and HCs cohorts.

Peripheral Blood Mononuclear Cell Isolation and DNA Extraction
PBMCs were isolated by Ficoll-Histopaque gradient centrifugation (Sigma-Aldrich, St. Louis, MO, USA). PBMCs were washed twice using phosphate-buffered saline (PBS 1X) and stored at −80 • C in fetal bovine serum (FBS) and dimethyl sulfoxide (DMSO) (Sigma-Aldrich, St. Louis, MO, USA). Genomic DNA extraction from PCa and HCs samples was performed using the DNeasy Blood and Tissue Kit (Qiagen, CA, USA). After thawing, cells were washed twice in PBS 1X to remove DMSO and FBS traces and then resuspended in 200 µL PBS 1X before proceeding according to the manufacturer's instructions. The final DNA concentration of each sample was measured using Nanodrop One (Thermo Scientific, Waltham, MA, USA), and DNA quality was assessed by calculating the ratios of absorbance at A260/A280 and A260/A230.

Library Preparation and Genomic DNA Sequencing
Library preparation and DNA sequencing were performed at the CRS4 Research Center (Pula, Sardegna, Italy). Thirty-five genes involved in both innate and adaptive immune responses as well as PCa-associated genes were sequenced. Two AmplySeq primer pools were designed and purchased from Illumina for the gene panel analysis. The PCR protocol and reagents were adapted from QIAGEN Application Note (Qiagen 1104745 10/2016) to amplify the entire length of each gene. The PCR products were quantified using Qubit, and equimolar amounts were pooled. Gene-panel PCR pools were grouped to achieve the appropriate sequencing coverage. Libraries were obtained using Nextera DNA Flex with 100 ng of DNA and indexed with IDT for the Illumina Nextera DNA UD Indexes Primer Set (Illumina, San Diego, CA, USA). PCR products were purified with 1X AMPure XP beads (Beckman Coulter, Brea, CA, USA), and the libraries were quantified using Qubit. A loading pool consisting of 90 samples (47 PCa and 43 HCs) was diluted to 9 pM before sequencing using the MiSeq Reagent Kit v3 600-cycle (Illumina). The BaseSpace Sequence Hub (Illumina) was used to conduct demultiplexing and fastq file generation, and two distinct analyses were performed for each sample. Gene panel sequences were analyzed using an in-house pipeline.

Statistical Analysis
The allele and genotype frequency estimates were calculated using (http://ihg2.hel mholtzmuenchen.de/cgi-bin/hw/hwa1.pl, accessed on 1 January 2022). Fisher exact test, chi-square test, and both allelic and genotypic odds ratios (ORs) were used to compare genetic variants between PCa and HC cohorts. To test the relationship between cancer onset and nucleotide substitution for a given locus, multivariate regression analysis was carried out by applying backward selection, in which non-significant covariates were removed step by step until the final model was obtained. A GLMM implemented in the R package lme4 [54] in the R environment (version 4.1.1; R Core Team available at https://www.r-project.org/, accessed on 1 January 2022), was used to perform the analysis. Odds ratios with 97.5% confidence interval were calculated using the R-package "oddsratio" (available at https://cran.r-project.org/web/packages/oddsratio/, accessed on 1 January 2022).