Variants of SCARB1 and VDR Involved in Complex Genetic Interactions May Be Implicated in the Genetic Susceptibility to Clear Cell Renal Cell Carcinoma

The current data are still inconclusive in terms of a genetic component involved in the susceptibility to renal cell carcinoma. Our aim was to evaluate 40 selected candidate polymorphisms for potential association with clear cell renal cell carcinoma (ccRCC) based on independent group of 167 patients and 200 healthy controls. The obtained data were searched for independent effects of particular polymorphisms as well as haplotypes and genetic interactions. Association testing implied position rs4765623 in the SCARB1 gene (OR = 1.688, 95% CI: 1.104–2.582, P = 0.016) and a haplotype in VDR comprising positions rs739837, rs731236, rs7975232, and rs1544410 (P = 0.012) to be the risk factors in the studied population. The study detected several epistatic effects contributing to the genetic susceptibility to ccRCC. Variation in GNAS1 was implicated in a strong synergistic interaction with BIRC5. This effect was part of a model suggested by multifactor dimensionality reduction method including also a synergy between GNAS1 and SCARB1 (P = 0.036). Significance of GNAS1-SCARB1 interaction was further confirmed by logistic regression (P = 0.041), which also indicated involvement of SCARB1 in additional interaction with EPAS1 (P = 0.008) as well as revealing interactions between GNAS1 and EPAS1 (P = 0.016), GNAS1 and MC1R (P = 0.031), GNAS1 and VDR (P = 0.032), and MC1R and VDR (P = 0.035).


Introduction
Kidney cancer is according to the International Agency for Research on Cancer the third fastest increasing carcinoma with the 30% growth noted in years 1993-2009. Clear cell renal cell carcinoma (ccRCC) is the most often type of kidney cancer in humans [1,2]. Several factors including obesity, hypertension, and smoking have been shown to increase risk of developing renal cell carcinoma [3,4]. It has also been noted that this cancer is twice more frequent in males than in females [5]. So far only few genes have been associated with increased risk of RCC. The data obtained from three genome wide association (GWA) studies performed on European ancestry populations are inconsistent. First GWAS revealed three loci located on 2p21, 11q13.3, and 12q24.31 to be associated with RCC [6]. The gene EPAS1 located in the region 2p21 has been previously suggested as a candidate gene in RCC due to its overexpression in this carcinoma [7]. The region 12q24.31 contains the SCARB1 gene that encodes receptor involved in uptake of high-density lipoprotein cholesterol. Fine mapping analysis of a 120 kb area near EPAS1 has confirmed association of 2p21 region with RCC but stronger signal was noted outside EPAS1 [8]. The second GWAS has implicated the ITPR2 gene on 12p.11.23 as a novel susceptibility locus for RCC while no other loci reached genome-wide significance in this extended population sample [9]. The third GWA study identified additional association with RCC for ZEB2 gene encoding zinc finger E-box-binding 2 BioMed Research International homeobox 2 [10]. The recent GWAS conducted on samples from African Americans has revealed the region 11q13. 3 containing no genes to be associated with RCC, confirming the outcome obtained before for Europeans [6,11]. Other loci implicated in RCC include GNAS1 showing a prognostic value [12] and several genes which have been associated with RCC risk. These latter ones include two genes from the vitamin D pathway, namely, VDR and RXRA which encode proteins forming a complex involved in vitamin D dependent transcription regulation [13] and polymorphism in the promoter of the VEGFA gene [14]. Two studies have reported associations with RCC in Chinese population showing polymorphism in the promoter of the BIRC5 gene [15] and variation in the XRCC6 locus, to be the susceptibility candidates for renal cell carcinoma [16]. One study indicated insertion/deletion polymorphism in the promoter region of CASP8 encoding a key regulator of apoptosis, caspase-8, to be a factor in RCC susceptibility and metastasis [17]. To further investigate the problem of genetic susceptibility to RCC we analysed 35 candidate polymorphisms within genes previously implicated in RCC development and prognosis. Additionally, 5 SNPs located in novel gene ZC3H12A encoding MCP-1-induced protein (MCPIP1) involved in a regulation of inflammatory reaction were included in the analyses. MCPIP1 acts as an RNAse regulating stability of some mRNA coding for proinflammatory cytokines and also regulates activity of transcription factors, which are key regulators in carcinogenesis (e.g., NF-B) [18][19][20][21]. Therefore, we hypothesize that variation in ZC3H12A may be associated with the risk of ccRCC development. Performed analyses included examination of main effects associated with all the selected 40 polymorphisms and the reconstructed haplotypes in EPAS1, RXRA, and VDR as well as a thorough search for epistatic effects using multifactor dimensionality reduction (MDR) and logistic regression methods.

Study Population and Candidate Genetic
Loci. This research was approved by the bioethical commission of the Regional Research Council in Kraków, number 68/KBL/OIL/ 2011. All the subjects gave an informed consent. The study sample comprised 167 ccRCC patients treated with a surgery in the Centre of Oncology in Kraków and 200 healthy individuals frequency matched by sex and age. Tissues collected from ccRCC patients were subjected to DNA extraction using silica-based method with GeneJET Genomic DNA Purification Kit (Thermo Fisher Scientific Inc., Waltham, MA) using the manufacturer's protocol. DNA from the buccal swabs collected from the control group was extracted by silica adsorption using the NucleoSpin Tissue extraction kit (Macherey-Nagel GmbH & Co. KG, Germany) as previously described [22]. Forty candidate polymorphisms including eight SNPs previously associated with renal cell carcinoma through the 4 genome wide association studies [6,[9][10][11] and seven SNPs associated with RCC in the selected paramount candidate gene approach association studies [12][13][14][15][16][17] were included in the research. Additionally, nine most often investigated SNPs in the VDR gene were examined due to a suggested association of this locus with RCC [13]. Eleven most relevant SNPs in MC1R, which is the important melanoma associated gene, were included due to a postulated increased risk of RCC in melanoma patients [2]. The remaining 5 SNPs are located in the ZC3H12A gene, which is involved in regulation of inflammatory reaction [20]. Details on the selected forty polymorphisms are given in Table 1.

Association Testing.
Univariate associations between single selected polymorphisms and ccRCC status were assessed through binary logistic regression with PASW statistics version 21 computer software (SPSS Inc., Chicago, IL, USA). Odds ratios (ORs) with corresponding 95% confidence intervals (CIs) and values were assigned. Additive, recessive, and dominant models of allele categorization were tested. In case of the MC1R gene, polymorphisms were divided in two groups according to their known influence on pigmentation phenotype. Strong variants marked with " " (N29insA, rs1805006, rs1805007, rs11547464, rs1805008, rs1805009, and Y152OCH) have been considered as highpenetrance since they significantly diminish the receptor performance and have strong phenotypic effect (OR values of ∼50 or higher). Weak variants marked with " " are lowpenetrance and have weaker effect on the receptor performance and therefore weaker phenotypic effect (OR below ∼ 10). Such a categorization has been applied in several previous studies (e.g., [27]). For " " variants, three states were considered according to the existence of major function mutations (0 = no " " variant carriers, 1 = one " " variant carriers, and 2 = two " " variant carriers) and the same approach was used in case of " " variants. Binary logistic regression was also used to evaluate the association of the inferred haplotypes (with frequency exceeding 0.5%) under the assumption of additive inheritance mode. Association analyses were conducted considering value of <0.05 level as statistically significant; however associations were also tested with false discovery rate (FDR) correction for multiple comparisons (http://www.seu.ac.lk/cedpl/student download/Benja-miniHochberg.xlsx). Finally, multivariate logistic regression was applied including simultaneous analysis of all the polymorphisms and haplotypes as well as disclosed interactions.

Epistasis Examination.
Statistical locus-locus interactions were examined using two methods, that is, multifactor dimensionality reduction (MDR v2.0 beta 8.1 http://www.epistasis.org/) and binary logistic regression (PASW statistics v.21). MDR approach enables evaluation of epistatic effects through assignment of all genotype combinations between interacting loci to one of the two groups: "high-risk" or "low-risk" of disease development based on the ratio of cases to controls and the analysis is performed as described elsewhere [28]. The results of MDR analysis were interpreted using entropy-based approach which is a nonparametric method for estimation of the benefit in information gain from interaction between two analyzed attributes [29]. In addition to MDR analysis, the parametric binary logistic regression method was applied. Pairwise interactions were tested by introducing interaction terms into the two-factor logistic regression models using forward selection strategy. ORs with 95% CIs and respective values were estimated for interactions disclosed with logistic regression.

Meta-Analysis.
Meta-analysis for rs4765623 in SCARB1 was performed with MetaEasy v1.0.5 [30] and involved data from our population and data for three different sample sets (IARC, NCI, Replication) from the study of Purdue et al. (2011) [6]. OR value for the combined data was calculated with the fixed effects model and value was estimated based on 95% CI with the method described in Altman and Bland [31]. Heterogeneity was investigated using Cochrane Q value and 2 statistics which assessed the consistency in meta-analysis studies. 2 ranges from 0% to 100% where 0% indicates no heterogeneity whereas increasing value of 2 indicates increasing level of heterogeneity.

Population Analyses and Haplotype Inference.
No significant departures from the Hardy-Weinberg equilibrium were noted after Bonferroni's correction for multiple testing ( > 0.0013) for the examined loci. Results of LD analysis among DNA variants investigated in EPAS1, RXRA, and VDR are presented in Supplementary Figure 1. Two haplotype blocks (rs11894252-rs7579899; rs9679290-rs4953346) with strong LD ( 2 > 0.9) were found in the EPAS1 gene, moderate linkage disequilibrium was revealed for two SNPs rs748964-rs3118523 ( 2 = 0.65) in the RXRA gene whereas in the VDR gene strong linkage for two SNPs rs4516035-rs7139166 ( 2 > 0.9) and 4 SNPs-haplotype block comprising rs739837-rs731236-rs7975232-rs1544410 ( 2 > 0.5) were detected. From 7 haplotypes reconstructed with PHASE program for all SNPs in EPAS1, 4 were very frequent. In case of RXRA 3 frequent and 1 rare haplotypes were reconstructed. In the block number 1 in VDR, comprising 4 SNPs, 10 haplotypes were reconstructed but only 6 were found to have frequency above 0.5%. In case of block number 3 in VDR 11 haplotypes were reconstructed with 7 exceeding the 0.5% threshold. The reconstructed haplotypes and their frequencies are listed in Supplementary Table 3.

Association and Interaction Analyses.
Samples collected from patients and controls were frequency matched by sex (58.1% of males in the ccRCC group and 51% in the control group) and age (the mean age of the participants was 64 years in the ccRCC group and 62 years in the control group). Univariate analyses performed for the particular 40 polymorphisms under study revealed rs4765623 in the SCARB1 gene to be associated with ccRCC in our population but only considering dominant mode of inheritance ( = 0.016) and this result did not pass FDR procedure. According to logistic An association testing performed for the reconstructed haplotypes revealed that CAGT haplotype in VDR gene (comprising rs739837-rs731236-rs7975232-rs1544410) observed with frequency 0.68% might be a risk factor for ccRCC ( = 0.012). This result was also found to be insignificant after applying the FDR procedure. Notably, however, CAGT haplotype was observed in 5 among 167 cases of ccRCC but was completely absent in the control group of 200 individuals. Haplotypes reconstructed in the EPAS1 and RXRA did not reveal any association with ccRCC under the significance threshold of 0.05.
The study revealed several epistatic effects involving the abovementioned DNA variants in SCARB1 and VDR, which confirmed suggestive association trends with ccRCC. According to MDR analysis, ccRCC is best explained by the four-factor model comprising rs4765623 in SCARB1, rs4953346 in EPAS1, rs9904341 in BIRC5, and rs7121 in GNAS1 with balanced accuracy (BA) of 0.5828, cross-validation consistency (CVC) of maximum level 10/10, and permutation testing value of 0.036 ( Table 2). Analysis of interaction dendrogram provided by MDR confirmed implication of these factors in epistatic effects showing strong positive (synergistic) interactions between rs7121 in GNAS1, rs9904341 in BIRC5, and rs4765623 in SCARB1 (Figure 1). Entropy-based interaction graph confirmed strong positive effects between rs7121 and rs9904341 and between rs7121 and rs4765623. According to the entropy-based analysis the largest main effect is attributed to rs4765623 in SCARB1. This position removes 1.28% of entropy, what is understood to remove 1.28% of "uncertainty" in prediction of ccRCC status. The position rs7121 in GNAS1 eliminates 0.43% of entropy, whereas the interaction between these two factors removes additional 0.55% of entropy which is not removed by any of these factors treated individually. In case of the second interaction revealed, between rs7121 in GNAS1 and rs9904341 in BIRC5 positive value of entropy brought by interaction equaled 0.62% which indicates additional benefit in information gain when considering epistasis between these factors ( Figure 2).
Finally, multivariate logistic regression involving simultaneous analysis of all the studied polymorphisms, haplotypes, and gene-gene interactions revealed the following factors to be relevant in susceptibility to ccRCC: CAGT haplotype in   VDR ( = 0.005), interaction between rs4765623 in SCARB1 and rs9679290 in EPAS1 ( = 0.007), interaction between variants " " in MC1R and rs7975232 in VDR ( = 0.010), and interaction between rs7121 in GNAS1 and rs2228570 in VDR ( = 0.041) ( Table 4). This final logistic regression model has the significance of 2 = 6.98 × 10 −6 and explains 9.2% of the total risk in ccRCC development.

Discussion
The picture arising from the available genetic data for renal cell carcinoma is ambiguous with no single polymorphism showing a record of repeatable and reliable association with RCC in independent studies. Here, we analyzed a group of 167 ccRCC patients and 200 healthy controls examining 40 candidate polymorphisms in 14 loci selected from the literature [2,6,9,10,[12][13][14][15][16][17][18][19].
The study revealed rs4765623 in the SCARB1 gene to show association trend with ccRCC in the studied population sample. The SCARB1 gene encodes the scavenger receptor class B member 1 which binds to the high-density lipoprotein (HDL) cholesterol and therefore mediates cellular cholesterol homeostasis. Cholesterol transfer towards liver for excretion is known as a protective mechanism against the development of atherosclerosis [32]. The role of SCARB1 in cancer biology is not well investigated but it has been recently suggested that SCARB1 receptor may be involved in multiple functions [33] and that cholesterol entry through SCARB1-HDL can activate signal transduction pathways that can regulate cellular proliferation and migration and therefore carcinogenesis [34]. Association of rs4765623 position in SCARB1 with RCC was first identified in GWAS study performed on European population by Purdue et al. in 2011 [6]. The revealed effect was weak (OR = 1.15) but significant ( = 2.6 × 10 −8 ). However, this association has not been confirmed in replication studies performed on Chinese population [35,36]. Our results confirmed association of the rs4765623 minor allele A ( = 36.1%) with susceptibility to ccRCC in European (Polish) population sample and can serve as an independent replication to further establish the association for rs4765623 for renal cell carcinoma risk. The effect identified in our study (OR = 1.7) was stronger than in the study of Purdue but this value was observed when dominant mode of minor allele classification was applied. Frequency of A risk allele reported in Purdue et al. was similar to our study (34%). Purdue et al. has reported significant heterogeneity for rs4765623 in SCARB1 ( = 0.03) among the three studied groups [6]. Meta-analysis including our data also suggests moderate interpopulation heterogeneity with ∼62% 2 value but the significance decreased from 0.03 to 0.048. The effect size of rs4765623 in SCARB1 for the combined data, including our population, was very similar to that obtained in Purdue et al.
Our study also revealed suggestive association trend of VDR haplotype comprised of four SNPs, namely, rs739837, rs731236, rs7975232, and rs1544410, with ccRCC risk. VDR gene encodes vitamin D receptor which is involved in vitamin D turnover. The active form of vitamin D [1,25(OH) 2 D 3 ] is considered to be a protective factor against cancer development due to its role in cell growth regulation, differentiation, programmed cell death, angiogenesis, and moderation of gene transcription [37,38]. Indeed, several functional studies have shown that the level of the active form of vitamin D is significantly lower in ccRCC patients than in healthy controls [39,40]. The active form of vitamin D binds to the vitamin D receptor which is predicted to regulate expression of over 2000 genes in humans [41]. Notably, variation in VDR gene has been implicated in many carcinomas (e.g., [13,[42][43][44][45][46][47]). However, few reports are available on association between VDR polymorphism and renal cell carcinoma and the obtained results are inconsistent [48,49]. Most of the previous studies have been performed on Asian populations [50,51] and recent meta-analyses showed that association of VDR polymorphism with ccRCC may be population specific [52,53]. Polymorphisms included in the haplotype associated with ccRCC in our study are located in the 3 region of VDR gene and are supposed to influence mRNA level [54]. Allele T of rs731236 (TaqI) polymorphism was associated with increased risk of RCC in a Japanese population sample [50] and then confirmed in a population of European descent [13]. However, several other studies have not confirmed that association (e.g., [51]). Our research revealed increased risk of RCC development associated with allele T in rs731236 when considered in the haplotype with three other SNPs. Interestingly, the identified risk haplotype CAGT was absent in a group of 200 healthy controls. Moreover, we analysed VDR variation in additional 300 controls and did not observe this haplotype ( = 0.008, data not shown). Association of a haplotype in the VDR gene with RCC was also noted by Karami et al. but in that study different region of the VDR gene was involved [13].
The discovered main association effects need to be interpreted as trends only because they did not pass FDR procedure. Small number of samples was an obvious limitation of this study that could prevent detection of weak genetic effects known to be typical for ccRCC. However, nonreplication may have also been caused by interpopulation heterogeneity [55] which was observed in meta-analysis for rs4765623 in SCARB1. Our sample size was sufficient to detect the true effects of OR = 1.803-1.941 (depending on minor allele frequency) with 80% probability for 19 out of 40 tested DNA variants. Among variants associated with RCC through GWA studies the strongest effect has been shown so far for rs7105934 (11q13.3) with OR equaling 1.41 in European samples [6] and 1.79 in African American samples [11]. It seems therefore that most of the known genetic RCC risk factors might have escaped detection in this study. However, besides main effects of the studied loci we also used statistical methods to test interactions between them. Notably, the SCARB1 locus was implicated by the MDR method in complex genegene interactions in the studied population sample. The significance of epistasis in the determination of complex traits is currently often emphasized and it has been suggested that the impact of genetic interactions may even outrank an independent main effect of single susceptibility loci [56]. Importantly, the applied MDR approach was developed as nonparametric and genetic model-free data mining strategy for epistasis identification in studies dealing with relatively small sample size [28,57]. It has been postulated that this method is less prone to I type error and thus more reliable than logistic regression. MDR showed that four-factor model comprising two synergistic interactions, the first between rs7121 in GNAS1 and rs9904341 in BIRC5 and the second between rs7121 in GNAS1 and rs4765623 in SCARB1, best explains susceptibility to ccRCC. The GNAS1 gene encodes a s subunit of heterotrimeric G protein that is required for activation of adenylyl cyclase and generation of cAMP and plays the key role in multiple signal transduction pathways, linked to proapoptotic processes in cancer cells. Modulated expression of this gene has been associated with various disorders (e.g., [58,59]). GNAS1 has been linked to RCC in two independent studies [12,60]. The BIRC5 gene encodes survivin, a protein preventing apoptotic cell death. Overexpression of this gene has a continuous literature record of its prognostic significance in many carcinomas (e.g., [61,62]). Increased expression of survivin has also been reported in RCC [63]. Recently, promoter mutation rs9904341 in BIRC5 has been associated with susceptibility to RCC in Asians [15]. There is evidence that this polymorphism regulates transcriptional activity increasing survivin level [64]. Notably, there are evidences that there is a direct dependency between GNAS1 and BIRC5 products. Subunit s activates STAT3 transcription factor inducing signal through the PKA kinase, JNK, and PI3 K [65]. In turn, STAT3 is important in activation of survivin encoded by BIRC5 gene (e.g., [66]). Therefore, synergistic interaction of GNAS1 and BIRC5 in susceptibility to ccRCC discovered in this study seems to be justified. In one of the recent studies utilizing ccRCC data provided by the Cancer Genome Atlas (TCGA), another member of the same group of inhibitors of apoptosis BIRC7 has been listed among genes differentially expressed in ccRCC [67]. The MDR method also revealed that rs7121 in GNAS1 is also implicated in synergistic interaction with rs4765623 in SCARB1. This interaction was further confirmed by logistic regression ( = 0.041). Interestingly, logistic regression also revealed additional and more statistically significant epistatic effect with SCARB1 as a component, that is, interaction between rs4765623 in SCARB1 and rs9679290 in EPAS1 ( = 0.008). Although a direct dependency between EPAS1 and SCARB1 or their protein products is unclear, it is worth noting that both these loci are related to the processes connected with angiogenesis. The SCARB1 gene, as it was mentioned before encodes receptor that binds to the high-density lipoprotein (HDL). It has been suggested that beyond its best known function in cholesterol mediation, HDL is also involved in many different activities like anti-inflammatory, antiapoptotic processes and a variety of endothelial behaviors and therefore angiogenesis (e.g., [68,69]). The second locus involved in the revealed interaction is EPAS1 which was associated with ccRCC in a large GWA study performed by Purdue et al. [11]. This gene encodes a hypoxia-inducible factor 2 (HIF2 ) which belongs to the transcription factors responsible for induction of genes controlled by oxygen. Under normal conditions, the level of HIFs is regulated by ubiquitinase complex, comprising a von Hippel-Lindau tumor suppressor protein (pVHL). In normoxia HIF factors are ubiquitinated and degraded. With a drop of oxygen level (e.g., in tumors), stabilization of HIF1 and HIF2 (encoded by EPAS1) occurs, which is what results in induction of transcription of many genes encoding proteins involved in angiogenesis [70]. Inactivation of VHL gene is observed in approximately 90% of patients with clear cell RCC and leads to accumulation of HIF1 and HIF2 proteins and in consequence to stimulation of expression of protooncogenes TGF and c-Met [7,71,72].
The GNAS1 gene was found in our study to be implicated in other epistatic effects with MC1R ( = 0.031), EPAS1 ( = 0.016), and VDR ( = 0.032). The MC1R gene encodes melanocortin 1 receptor which is involved in the regulation of melanin pigment synthesis which has also been found to increase the risk of developing melanoma (e.g., [44,73,74]). It is well established that patients with cancer are at higher risk to develop multiple cancers and some examples of melanoma and RCC coexistence have been reported (e.g., [75,76]). However, risk factors common for both cancers are only partially explained and include mutations in MITF gene [77]. Interestingly, product of MC1R gene acts through Gprotein. Therefore, there is a direct molecular dependence between MC1R and GNAS1 products. Importantly, it has been demonstrated that MC1R is expressed in the kidney cells and shown that treatment with MC1R agonists ameliorated kidney diseases in rats with passive Heymann nephritis [78]. Finally, logistic regression also indicated interaction between the MC1R and VDR genes ( = 0.035). It is not clear how this epistatic effect could affect ccRCC development, but interaction between these two genes has been suggested to influence human pigmentation [23,79]. Interestingly, in the recent TCGA study another gene involved in melanogenesis and increased risk of melanoma that is the TYRP1 locus [74] has been selected as differentially expressed in ccRCC [67]. The overlap between risk factors in melanoma and ccRCC seems to be interesting and worth further investigation.
In conclusion, position rs4765623 in SCARB1 and haplotype in VDR showed suggestive association trends with ccRCC susceptibility in the studied Polish population. Moreover using MDR and logistic regression methods, a complex network of interactions involving six genes previously implicated in RCC (SCARB1, GNAS1, BIRC5, EPAS1, VDR, and MC1R) was discovered. Due to a relatively small sample number these epistatic effects should be further studied and confirmed on a larger cohort. The risk haplotype in VDR and the gene-gene interactions SCARB1-EPAS1, GNAS1-VDR, and MC1R-VDR were included in the final logistic regression model explaining 9.2% of the total risk in ccRCC development.

Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.