Vitamin D receptor-binding site variants affect prostate cancer progression

Vitamin D is an important modulator of cellular proliferation through the vitamin D receptor (VDR) that binds to DNA in the regulatory sequences of target genes. We hypothesized that single nucleotide polymorphisms (SNPs) in VDR-binding sites might affect target gene expression and influence the progression of prostate cancer. Using a genome-wide prediction database, 62 SNPs in VDR-binding sites were selected for genotyping in 515 prostate cancer patients and the findings were replicated in an independent cohort of 411 patients. Prognostic significance on prostate cancer progression was assessed by Kaplan-Meier analysis and the Cox regression model. According to multivariate analyses adjusted for known predictors, HFE rs9393682 was found to be associated with disease progression for localized prostate cancer, and TUSC3 rs1378033 was associated with progression for advanced prostate cancer in both cohorts. Vitamin D treatment inhibited HFE mRNA expression, and down-regulation of HFE by transfecting small interfering RNA suppressed PC-3 human prostate cancer cell proliferation and wound healing ability. In contrast, vitamin D treatment induced TUSC3 expression, and silencing TUSC3 promoted prostate cancer cell growth and migration. Further analysis of an independent microarray dataset confirmed that low TUSC3 expression correlated with poor patient prognosis. Our results warrant further studies using larger cohorts. This study identifies common variants in VDR-binding sites as prognostic markers of prostate cancer progression and HFE and TUSC3 as plausible susceptibility genes.


INTRODUCTION
Epidemiological studies have demonstrated that low sunlight exposure and poor vitamin D status at higher latitudes account for an elevated risk of a number of cancers, including prostate cancer [1,2]. Vitamin D can be obtained from the diet; however, it is majorly synthesized in the skin using solar irradiation. The biologically active form of vitamin D 3 , 1α,25-dihydroxyvitamin D 3 (1,25-VD), is produced by 25-hydroxylase in the liver, followed by 1α-hydroxylase in the kidney. The actions of 1,25-VD are mediated by the vitamin D receptor (VDR), a ligandactivated transcription factor. Upon activation by 1,25-VD, VDR forms a heterodimer with the retinoid X receptor, and binds to vitamin D response elements (VDREs) in the promoters of vitamin D-responsive genes [3]. Prostate cells express VDR and vitamin D metabolizing enzymes, and can respond to 1,25-VD. A volume of data supports multipronged effects of 1,25-VD in the prevention of prostate cancer progression by induction of detoxifying enzymes [4], cell cycle arrest [5], and apoptosis [6], as well as inhibition of prostate cancer cell invasion [7] and angiogenesis [8].
Genetic variants in VDREs may affect VDR-VDRE interactions, thereby resulting in altered expression of target genes and consequent cancer progression. However, no study to date has investigated the single nucleotide polymorphisms (SNPs) in VDR binding sites and their relationship to the clinical outcomes of prostate cancer. Accordingly, we conducted a two-stage study to evaluate the associations of VDRE SNPs with prostate cancer progression, and further assessed the functional relevance of candidate genes of interest, as illustrated in Supplementary Figure 1.

RESULTS
The clinical characteristics of patients in the discovery and replication cohorts and the association with disease progression are shown in Table 1. For localized prostate cancer, 45 (30.0%) and 75 (43.9%) patients experienced disease progression after radical prostatectomy (RP) during the median follow-up of 23 and 30 months in the discovery and replication cohorts, respectively. Prostate-specific antigen (PSA) at diagnosis, pathologic Gleason score, and pathologic stage were significantly associated with cancer progression in both cohorts. In the advanced prostate cancer group, 271 (74.5%) and 180 (75.3%) patients had disease progression after androgen deprivation therapy (ADT) during the median follow-up of 61 and 57 months in the discovery and replication cohorts, respectively. PSA at ADT initiation, and PSA nadir were significantly associated with cancer progression in both cohorts. Gleason score, clinical stage at diagnosis, and treatment modality were also associated with progression in the discovery cohort, but only weakly associated in the replication cohort.
Of the 62 SNPs in VDREs analyzed in the discovery cohort, six SNPs were associated with time to progression (TTP) (Supplementary Table 1) in localized prostate cancer patients. rs9393682 was found to be significantly associated with TTP in the same direction as the discovery cohort in an independent replication cohort. In combined analysis, rs9393682 was associated with a per-allele hazard ratio (HR) of 1.79 [95% confidence interval (CI), 1.38-2.33; P < 0.001; Table 2 and Figure 1A]. This association remained significant (P = 0.006) after adjusting for age, PSA at diagnosis, pathologic Gleason score, and stage. Furthermore, the outcome prediction model based on clinical factors (age, PSA at diagnosis, pathologic Gleason score, and stage) plus rs9393682 was significantly improved over the model with clinical factors only, as indicated by the likelihood ratio test (χ 2 69.63, df 1, P < 0.001).
For the advanced prostate cancer group, four VDRE SNPs were associated with TTP during ADT (Supplementary Table 2). Only rs1378033 showed significant correlation with a decreased risk of disease progression in both the discovery and replication cohorts, and upon combined analysis (HR 0.75, 95% CI 0.63-0.89, P = 0.001; Table 3 and Figure 1B). After adjusting for known predictors, the association remained significant (P = 0.004). The model based on clinical factors plus rs1378033 was significantly improved over the model with clinical factors only, as indicated by the likelihood ratio test (χ 2 380.51, df 1, P < 0.001).
rs9393682 is located in the intergenic region between HIST1H1C (histone cluster 1 H1 family member c) and HFE (hemochromatosis), and rs1378033 is located in the intron of SGCZ (sarcoglycan zeta) and 5′ of TUSC3 (tumor suppressor candidate 3). To gain initial insight for further analysis, we investigated if rs9393682 and rs1378033 were associated with differential expression of nearby genes expression in prostate tissues. The Genotype-Tissue Expression (GTEx) database showed a significant trend for increased HFE expression in rs9393682 C www.impactjournals.com/oncotarget Abbreviations: IQR, interquartile range; PSA, prostate-specific antigen; RP, radical prostatectomy; CI, confidence interval; ADT, androgen-deprivation therapy; RT, radiotherapy.
a P value was calculated by the log-rank test or Cox regression for disease progression. b Median follow-up time and 95% CIs were estimated with the reverse Kaplan-Meier method. Oncotarget 74123 www.impactjournals.com/oncotarget allele carriers (P = 0.006; Figure 2A). rs9393682 and rs1378033 fall within putative VDREs, and may alter 1,25-VD-mediated gene regulation. The effects of 1,25-VD on the mRNA expression levels of HIST1H1C, HFE, SGCZ, TUSC3, and CYP24A1 (cytochrome P450 family 24 subfamily A member 1), a well-known 1,25-VD target gene, were examined using quantitative real-time polymerase chain reaction (qRT-PCR) in PC-3 human prostate cancer cells. As a positive control, 1,25-VD can strongly induce CYP24A1 expression ( Figure 2B). The expression of TUSC3 was also induced by 1,25-VD, but HFE was downregulated. We then assessed HFE and TUSC3 expression using a prostate cancer complementary DNA array containing 39 tissue samples. A significant downregulation of TUSC3 expression was observed in more advanced stage cancer samples (P = 0.038, Figure  2C). The prognostic values of HFE and TUSC3 on prostate cancer progression were further evaluated using a publicly available prostate cancer microarray dataset. Patients were dichotomized by HFE or TUSC3 gene expression using an optimization algorithm for the minimum P value. Low expression of TUSC3 was associated with shorter progression-free survival (P = 0.001; Figure 2D). Taken together, these observations indicate that 1,25-VD can stimulate TUSC3 expression, and down-regulation of TUSC3 in late stages of cancer may increase the risk of prostate cancer progression.
To better understand the roles of HFE and TUSC3 in prostate cancer, we used small interfering RNAs (siRNAs) to silence their expression in PC-3 cells, and evaluate the consequences of loss of HFE or TUSC3 function on cell proliferation and migration. As shown in Figure 3A, HFE or TUSC3 gene expression was effectively down-regulated by corresponding siRNAs, compared to that of the negative control (NC)-siRNA. Silencing of HFE inhibited cell proliferation, whereas knockdown of TUSC3 in PC-3 cells demonstrated a significant growth advantage over NC cells ( Figure 3B). In the wound healing assay, our data showed that HFE silencing retarded wound closure by 50%, but TUSC3 silenced PC-3 cells increased migration ability by 157% ( Figure 3C). These data demonstrated that both HFE and TUSC3 could regulate prostate cancer cell proliferation and migration.

DISCUSSION
We conducted a two-stage study to identify genetic variants in VDR binding sites corresponding with prognosis of patients with surgically resected localized prostate cancer and patients with ADT treated advanced disease. The associations of HFE rs9393682 and TUSC3 rs1378033 with disease progression were replicated across both stages of the study for localized and advanced prostate cancer, respectively. In addition, knockdown HFE and TUSC3 expression significantly influenced prostate cancer cell proliferation and migration. We also found that TUSC3 expression could be induced by 1,25-VD, and its expression was down-regulated in advanced stage cancer. These findings suggest that TUSC3 could be used as a prognostic marker for prostate cancer.
Based on the HaploReg data, rs9393682 and its correlated variants within a linkage disequilibrium block are positioned in the promoter elements, and rs1378033 is situated in a regulatory region containing DNase hypersensitivity peaks in different cell lines. The expression of two genes near these VDRE SNPs, HFE and TUSC3, can be regulated by 1,25-VD, suggesting that HFE and TUSC3 might be potential vitamin D Oncotarget 74124 www.impactjournals.com/oncotarget target genes and thus contribute to prostate cancer progression. HFE has been recognized as a receptor for β2-microglobulin (β2-M), and β2-M is a known growthpromoting gene for several human cancers, including prostate cancer [9,10]. β2-M interacts with HFE to modulate intracellular iron, activate iron responsive hypoxia inducible factor-1α signaling, and promote cancer bone metastasis through its induction of epithelialmesenchymal transition in cancer cells [11]. Consistent with these studies, our data showed that silencing HFE expression could markedly inhibit prostate cancer cell growth and migration (Figure 3). Although no association was observed between HFE expression and prostate cancer stage and progression ( Figure 2C and 2D), high HFE expression was significantly associated with shorter overall survival in multiple The Cancer Genome Atlas (TCGA) cancer sets, such as breast invasive carcinoma, brain lower grade glioma, and pancreatic adenocarcinoma (Supplementary Figure 2; prostate adenocarcinoma was not calculated because of low number of deaths). This discrepancy may be due to the relative small sample size of tumor tissues used to assess the association Oncotarget 74125 www.impactjournals.com/oncotarget between gene expression and patient prognosis. TUSC3 was identified as a potential tumor suppressor gene on chromosome 8p22, a common homozygously deleted region of the metastatic prostate cancer [12]. It has been described as a homologue of a subunit of the yeast oligosaccharyltransferase complex, which modulates glycosylation of proteins in the endoplasmic reticulum (ER) [13]. Loss of TUSC3 has been proposed to facilitate prostate cancer progression by increasing protein glycosylation, alleviating unfolded protein response and ER stress, and promoting Akt survival signaling [14]. Our data also showed that silencing TUSC3 expression could increase prostate cancer cell proliferation and migration (Figure 3). In addition, TUSC3 expression was downregulated in late stage cancers, and low TUSC3 expression was significantly associated with shorter TTP in prostate cancer, as well as decreased survival in several TCGA cancer sets (Supplementary Figure 2). However, no SNP is simultaneously associated with disease progression in both localized and advanced prostate cancers, suggesting that distinct biological pathways might be involved in different stages of the disease and the treatments adopted.
The associations of HFE rs9393682 and TUSC3 rs1378033 with prostate cancer progression was replicated across both cohorts, which would reduce false-positive findings in this study. In addition, our functional studies support these gene-disease associations. However, several limitations in the present study should be considered. First, we chose the TTP end point based on serum PSA, due to its biologic and clinical relevance. A rising PSA is often the first indication of the development of progressive disease and precipitates a change in therapy. Thus, the end point of PSA TTP most closely identifies the timing of treatment failure. Second, the modest sample size of both cohorts did not have optimal statistical power for discovering and replicating the association, so the observed P values did not reach a level that would avoid false-positives arising from multiple testing. Third, although we speculated that VDRE SNPs might affect VDR binding to genomic sequences, the specific mechanism was not Oncotarget 74126 www.impactjournals.com/oncotarget clear. However, the two susceptibility genes identified in this study, HFE and TUSC3, can indeed be regulated by vitamin D. Finally, our findings in this homogeneous Chinese Han population might not be applicable to other ethnic groups. Further functional characterizations and studies with larger patient cohorts are required to validate our findings.
This study shows that HFE rs9393682 and TUSC3 rs1378033 influence TTP in patients with prostate cancer. Silencing TUSC3 promotes prostate cancer cell proliferation and migration, and its expression is decreased in advanced stage cancer tissues, as well as in patients with poor prognosis. These results suggest that TUSC3 may function as a candidate susceptibility gene, and is a promising target for prostate cancer progression.

Patient recruitment and data collection
The study included 926 patients with prostate cancer divided into two independent cohorts ( Table 1). The discovery cohort comprised 150 men with localized prostate cancer who underwent RP as initial therapy, and 365 men with advanced cancer on ADT from the National Taiwan University Hospital located in northern Taiwan, as described previously [15][16][17]. The replication cohort was composed of 171 patients with localized prostate cancer, and 240 patients with advanced disease from the Kaohsiung Medical University Hospital, E-Da Hospital, and Kaohsiung Veterans General Hospital, all located in southern Taiwan. Detailed clinical information was obtained from patients' medical records. The primary outcome variable was TTP. TTP after RP for localized prostate cancer was defined as two consecutive PSA measurements of more than 0.2 ng/mL at an interval of more than three months, and the first of the consecutive rises was considered as the time of progression [18]. TTP after ADT for advanced prostate cancer was defined as a serial rise in PSA, at least 2 rises in PSA (> 1 week apart), greater than the PSA nadir [19]. Initiation of secondary hormone treatment for rising PSA was also considered as a progression event. This study was approved by the Institutional Review Board of Kaohsiung Medical University Hospital. Written informed consent was obtained from each patient, and the study was carried out in accordance with approved guidelines.

SNP selection and genotyping
Since transcription factors are known to regulate different genes in different cellular contexts [20], we used a genome-wide cis-regulatory module prediction database, PReMod (genomequebec.mcgill.ca/PReMod) [21], to identify putative VDREs in the whole human genome instead of the chromatin immunoprecipitation data. The PReMod algorithm predicts that a total of 11,342 sites within the human genome are bound by the VDR (canonical VDR position weight matrix: M00444; consensus: GGGKNARNRRGGWSA) [22]. We identified SNPs within VDREs by comparing two hexameric halfsites of these putative VDREs with HapMap SNPs CHB (Han Chinese, Beijing, China) data in the UCSC table browser (NCBI35/hg17) [23,24]. SNPs with a minor allele frequency less of than 0.10 in the HapMap CHB population were excluded, thus leaving 68 SNPs in VDREs that were initially selected for analysis. The threshold of 0.10 was chosen because it was considered to be the lowest minor allele frequency for a SNP with a relative risk of at least two being detectable with a sample size of 500.
Genomic DNA was extracted from patients' peripheral blood using the QIAamp DNA Blood Mini Kit (Qiagen, Valencia, CA) and stored at −80 °C until the time of the study. Genotyping was performed as described previously [18] at the National Center for Genome Medicine, Academia Sinica, Taiwan, using Agena Bioscience iPLEX matrix-assisted laser desorption/ ionization time-of-flight mass-spectrometry technology. The average genotype call rate for these SNPs was 95.2%. Ten samples were blindly duplicated for quality control and the genotype concordance was 99.7%. Six SNPs that significantly deviated from the Hardy-Weinberg equilibrium (P < 0.005) were removed, leaving 62 SNPs for further statistical analysis.
qRT-PCR analysis, lentiviral transduction, cell proliferation and wound healing assays, and bioinformatics analysis PC-3 (CRL-1435) human prostate cancer cell line was purchased from ATCC (Manassas, VA) and maintained in the recommended culture media. The identity of the cell line was checked by Cell ID System and Promega GenePrint 10 System through short tandem repeat analysis (Mission Biotech, Taipei, Taiwan). A comprehensive methods section is available in Supplementary Methods.

Statistical analysis
Patient clinicopathologic characteristics were summarized as number and percentage of patients or median and interquartile range of values. The association between patient characteristics with TTP was assessed by the log-rank test or Cox regression analysis. Individual SNPs were first assessed by association with TTP using the log-rank test under dominant, recessive, and additive models because the function of the SNPs was unknown. Cox proportional hazards regression was then conducted on each SNP as an isolated covariate with adjustment for known prognostic factors. In localized prostate www.impactjournals.com/oncotarget cancer, multiple explanatory variables included age, PSA at diagnosis, pathologic Gleason score, and stage, as previously defined [25]. In advanced prostate cancer, multiple explanatory variables included age, PSA at ADT initiation, biopsy Gleason score, clinical stage, PSA nadir, and treatment modality, as previously defined [26]. Heterogeneity between cohorts was evaluated by Cochran's χ 2 -based Q statistical test. If the results of the Q test were significant, a random-effects model was used to accommodate the diversity; otherwise, the combined HR was estimated using the fixed-effects model. The Statistical Package for the Social Sciences software version 22.0.0 (IBM, Armonk, NY) was used for statistical analyses. A two-sided P value of < 0.05 was considered statistically significant.