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

Authors: Ewelina Pośpiech; Janusz Ligęza; Wacław Wilk; Aniela Gołas; Janusz Jaszczyński; Andrzej Stelmach; Janusz Ryś; Aleksandra Blecharczyk; Anna Wojas-Pelc; Jolanta Jura; Wojciech Branicki
Journal: BioMed Research International
 (2015)

---

## Abstract

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


---

## Body

## 1. 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 geneEPAS1 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 theSCARB1 gene that encodes receptor involved in uptake of high-density lipoprotein cholesterol. Fine mapping analysis of a 120 kb area nearEPAS1 has confirmed association of 2p21 region with RCC but stronger signal was noted outsideEPAS1 [8]. The second GWAS has implicated theITPR2 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 forZEB2 gene encoding zinc finger E-box-binding 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 includeGNAS1 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 andRXRA which encode proteins forming a complex involved in vitamin D dependent transcription regulation [13] and polymorphism in the promoter of theVEGFA gene [14]. Two studies have reported associations with RCC in Chinese population showing polymorphism in the promoter of theBIRC5 gene [15] and variation in theXRCC6 locus, to be the susceptibility candidates for renal cell carcinoma [16]. One study indicated insertion/deletion polymorphism in the promoter region ofCASP8 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 geneZC3H12A 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–21]. Therefore, we hypothesize that variation inZC3H12Amay 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 inEPAS1,RXRA, andVDR as well as a thorough search for epistatic effects using multifactor dimensionality reduction (MDR) and logistic regression methods.


## 2. Materials and Methods

### 2.1. 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–11] and seven SNPs associated with RCC in the selected paramount candidate gene approach association studies [12–17] were included in the research. Additionally, nine most often investigated SNPs in theVDR gene were examined due to a suggested association of this locus with RCC [13]. Eleven most relevant SNPs inMC1R, 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 theZC3H12A gene, which is involved in regulation of inflammatory reaction [20]. Details on the selected forty polymorphisms are given in Table 1.Table 1
Data for the studied 40 polymorphisms.
Locus
Gene
Chromosome
Position (GRCh38)
Allele variants
MAF∗
Reference
rs34796867
ZC3H12A
1
37475773
C/T
T = 0.000
[20]
rs113322875
ZC3H12A
1
37477383
G/A
A = 0.000
[20]
rs34031609
ZC3H12A
1
37482928
G/A
A = 0.000
[20]
rs113655247
ZC3H12A
1
37483276
G/C
C = 0.000
[20]
rs17849897
ZC3H12A
1
37483451
C/T
T = 0.014
[20]
rs12105918
ZEB2
2
144450626
G/A
G = 0.051
[10]
rs11894252
EPAS1
2
46306237
C/T
T = 0.356
[6]
rs7579899
EPAS1
2
46310465
C/T
T = 0.357
[6]
rs9679290
EPAS1
2
46330505
G/C
C = 0.451
[6]
rs4953346
EPAS1
2
46331069
C/A
C = 0.452
[6]
rs3834129
CASP8
2
201232809:14
CTTACT/—
Del = 0.450
[17]
rs833061
VEGFA
6
43769749
C/T
C = 0.484
[14]
rs748964
RXRA
9
134442243
G/C
G = 0.144
[13]
rs3118523
RXRA
9
134443675
G/A
G = 0.202
[13]
rs7105934
11q13.3
11
69424973
C/T
T = 0.040
[6]
rs1049380
ITPR2
12
26336611
G/T
G = 0.436
[9]
rs739837
VDR
12
47844438
C/A
A = 0.471
[13]
rs731236
VDR
12
47844974
G/A
G = 0.354
[13]
rs7975232
VDR
12
47845054
G/T
T = 0.469
[13]
rs1544410
VDR
12
47846052
C/T
T = 0.361
[13]
rs2228570
VDR
12
47879112
G/A
A = 0.456
[13]
rs2238136
VDR
12
47883930
G/A
A = 0.365
[13]
rs4516035
VDR
12
47906043
C/T
C = 0.458
[13]
rs7139166
VDR
12
47906551
G/C
G = 0.454
[13]
rs11568820
VDR
12
47908762
C/T
T = 0.128
[13]
rs4765623
SCARB1
12
124836304
G/A
A = 0.361
[6]
N29insA
MC1R
16
89919345
C/A
A = 0.000
[2]
rs1805005
MC1R
16
89919436
G/T
T = 0.078
[2]
rs1805006
MC1R
16
89919510
C/A
A = 0.001
[2]
rs2228479
MC1R
16
89919532
G/A
A = 0.108
[2]
rs11547464
MC1R
16
89919683
G/A
A = 0.007
[2]
rs1805007
MC1R
16
89919709
C/T
T = 0.040
[2]
Y152OCH
MC1R
16
89919714
C/A
A = 0.000
[2]
rs1110400
MC1R
16
89919722
C/T
C = 0.026
[2]
rs1805008
MC1R
16
89919736
C/T
T = 0.068
[2]
rs885479
MC1R
16
89919746
C/T
T = 0.035
[2]
rs1805009
MC1R
16
89920138
G/C
C = 0.001
[2]
rs9904341
BIRC5
17
78214286
G/C
C = 0.343
[15]
rs7121
GNAS1
20
58903752
C/T
T = 0.465
[12]
rs132770
XRCC6
22
41621260
G/A
A = 0.153
[16]
∗MAF  =  minor allele frequency.


### 2.2. Genotyping

Polymorphisms were genotyped in 5 multiplex and 1 singleplex PCR reactions. SNPs in multiplexes were amplified using Qiagen Multiplex PCR kit (Qiagen, Hilden, Germany) and further analysed with single base extension (SBE) reactions using SNaPshot Multiplex Kit (Applied Biosystems, Foster City, CA). Multiplex 1 included 7 SNP positions (rs11894252, rs7579899, rs833061, rs3118523, rs748964, rs7105934, and rs4765623), multiplex 2 consisted of 6 SNPs (rs9679290, rs4953346, rs1049380, rs9904341, rs7121, and rs132770), and multiplex 3 comprised further 6 SNPs (rs34796867, rs113322875, rs34031609, rs113655247, rs17849897, and rs12105918). Sequences and concentrations of PCR and SBE primers are shown in Supplementary Tables  1 and  2, respectively in Supplementary Material available online athttp://dx.doi.org/10.1155/2015/860405. Amplifications were performed in 5 μL reaction volume consisting of 2.5 μL Qiagen Multiplex PCR mixture, 0.5 μL of Q solution, 0.5 μL of primer premix, and 1.5 μL of template DNA. The PCR reactions were carried out with the following temperature profile: 95°C/15 min, [94°/30 s, 58°C/90 s, 72°C/90 s] × 32, 72°C/10 min. The amplification products were cleaned with a mixture (1 : 1) of Exonuclease I (Exo I) and FastAP Thermosensitive Alkaline Phosphatase (FastAP) enzymes (Thermo Fisher Scientific Inc., Waltham, MA). The minisequencing reactions consisted of 0.5 μL of SNaPshot mix, 0.5 μL of SBE primer premix, 1 μL of purified PCR product, and 3 μL of DNase free water. The SBE temperature profile was as follows: [96°C/10 s, 50°C/5 s, 60°C/30 s] × 26. Products of the extension reactions were purified with 1 μL of FastAP Thermosensitive Alkaline Phosphatase enzyme and analyzed with ABI 3100 Avant Genetic Analyser (Applied Biosystems, Foster City, CA). The remaining 11 polymorphisms inMC1R (N29insA, rs1805005, rs1805006, rs2228479, rs11547464, rs1805007, Y152OCH, rs1110400, rs1805008, rs885479, and rs1805009) and 9 SNPs from theVDR gene (rs739837, rs731236, rs7975232, rs1544410, rs2228570, rs2238136, rs4516035, rs7139166, and rs11568820) were analysed in additional two multiplexes using protocols described elsewhere [22, 23]. The remaining insertion/deletion polymorphism inCASP8 was analysed in a single PCR reaction using HotStarTaq Master Mix Kit (Qiagen, Hilden, Germany) and following primer sequences 5′-6-FAM-AAACTTCTCCCATGGCCTCT-3′ and 5′-TATGAATGAGCCGAGGAAGG3′. Amplifications were carried out in 5 μL reaction volume consisting of 2.5 μL HotStarTaq Master Mix, 0.5 μL of primer premix (concentration of 0.125 μM), 0.5 μL of distilled water, and 1.5 μL of template DNA. The following temperature profile was applied: 95°C/15 min, [94°/30 s, 60°C/45 s, 72°C/60 s] × 30, 72°C/10 min. PCR products were directly analyzed by capillary electrophoresis using ABI 3100 Avant Genetic Analyser (Applied Biosystems, Foster City, CA).


### 2.3. Statistical Analyses

#### 2.3.1. Population Analyses and Haplotype Reconstruction

Genetic data obtained for the 40 selected polymorphisms were tested for deviations from Hardy-Weinberg equilibrium. Hardy-Weinberg analysis was performed with Arlequin version 3.1 software (http://cmpg.unibe.ch/software/arlequin3/). Haploview version 4.2 software was used to test the degree of linkage disequilibrium among DNA variants inEPAS1,RXRA, andVDR (http://www.broadinstitute.org/). Haplotype reconstruction and frequency estimation were done for SNPs located inEPAS1,RXRA, andVDR using PHASE version 2.1 computer software (http://stephenslab.uchicago.edu/software.html#phase). Haplotype reconstruction and LD analysis did not include polymorphisms in theMC1R gene because they all may individually affect receptor performance [24] and polymorphisms inZC3H12a as besides one SNP (rs17849897) they turned out to be monomorphic in the tested population. In the case ofVDR gene, polymorphisms were divided into 3 blocks following the obtained results and literature suggestions [25, 26]. Positions rs739837, rs731236, rs7975232, and rs1544410 were included in the first block. The second block was formed by a single rs2228570 position whereas in the third block the remaining rs2238136, rs4516035, rs713966, and rs11568820 positions were considered. Minimal ORs detectable with the power of at least 80% were calculated with Power and Sample Size Program (PS Program) version 3.1.2 (http://biostat.mc.vanderbilt.edu/wiki/Main/PowerSampleSize).


#### 2.3.2. 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) andP values were assigned. Additive, recessive, and dominant models of allele categorization were tested. In case of theMC1R gene, polymorphisms were divided in two groups according to their known influence on pigmentation phenotype. Strong variants marked with “R” (N29insA, rs1805006, rs1805007, rs11547464, rs1805008, rs1805009, and Y152OCH) have been considered as high-penetrance since they significantly diminish the receptor performance and have strong phenotypic effect (OR values of ~50 or higher). Weak variants marked with “r” are low-penetrance 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 “R” variants, three states were considered according to the existence of major function mutations (0 = no “R” variant carriers, 1 = one “R” variant carriers, and 2 = two “R” variant carriers) and the same approach was used in case of “r” 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 P 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/BenjaminiHochberg.xlsx). Finally, multivariate logistic regression was applied including simultaneous analysis of all the polymorphisms and haplotypes as well as disclosed interactions.


#### 2.3.3. Epistasis Examination

Statistical locus-locus interactions were examined using two methods, that is, multifactor dimensionality reduction (MDR v2.0 beta 8.1http://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 P values were estimated for interactions disclosed with logistic regression.


#### 2.3.4. Meta-Analysis

Meta-analysis for rs4765623 inSCARB1 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 P value was estimated based on 95% CI with the method described in Altman and Bland [31]. Heterogeneity was investigated using Cochrane Q P value and I
2 statistics which assessed the consistency in meta-analysis studies. I
2 ranges from 0% to 100% where 0% indicates no heterogeneity whereas increasing value of I
2 indicates increasing level of heterogeneity.


## 2.1. 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–11] and seven SNPs associated with RCC in the selected paramount candidate gene approach association studies [12–17] were included in the research. Additionally, nine most often investigated SNPs in theVDR gene were examined due to a suggested association of this locus with RCC [13]. Eleven most relevant SNPs inMC1R, 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 theZC3H12A gene, which is involved in regulation of inflammatory reaction [20]. Details on the selected forty polymorphisms are given in Table 1.Table 1
Data for the studied 40 polymorphisms.
Locus
Gene
Chromosome
Position (GRCh38)
Allele variants
MAF∗
Reference
rs34796867
ZC3H12A
1
37475773
C/T
T = 0.000
[20]
rs113322875
ZC3H12A
1
37477383
G/A
A = 0.000
[20]
rs34031609
ZC3H12A
1
37482928
G/A
A = 0.000
[20]
rs113655247
ZC3H12A
1
37483276
G/C
C = 0.000
[20]
rs17849897
ZC3H12A
1
37483451
C/T
T = 0.014
[20]
rs12105918
ZEB2
2
144450626
G/A
G = 0.051
[10]
rs11894252
EPAS1
2
46306237
C/T
T = 0.356
[6]
rs7579899
EPAS1
2
46310465
C/T
T = 0.357
[6]
rs9679290
EPAS1
2
46330505
G/C
C = 0.451
[6]
rs4953346
EPAS1
2
46331069
C/A
C = 0.452
[6]
rs3834129
CASP8
2
201232809:14
CTTACT/—
Del = 0.450
[17]
rs833061
VEGFA
6
43769749
C/T
C = 0.484
[14]
rs748964
RXRA
9
134442243
G/C
G = 0.144
[13]
rs3118523
RXRA
9
134443675
G/A
G = 0.202
[13]
rs7105934
11q13.3
11
69424973
C/T
T = 0.040
[6]
rs1049380
ITPR2
12
26336611
G/T
G = 0.436
[9]
rs739837
VDR
12
47844438
C/A
A = 0.471
[13]
rs731236
VDR
12
47844974
G/A
G = 0.354
[13]
rs7975232
VDR
12
47845054
G/T
T = 0.469
[13]
rs1544410
VDR
12
47846052
C/T
T = 0.361
[13]
rs2228570
VDR
12
47879112
G/A
A = 0.456
[13]
rs2238136
VDR
12
47883930
G/A
A = 0.365
[13]
rs4516035
VDR
12
47906043
C/T
C = 0.458
[13]
rs7139166
VDR
12
47906551
G/C
G = 0.454
[13]
rs11568820
VDR
12
47908762
C/T
T = 0.128
[13]
rs4765623
SCARB1
12
124836304
G/A
A = 0.361
[6]
N29insA
MC1R
16
89919345
C/A
A = 0.000
[2]
rs1805005
MC1R
16
89919436
G/T
T = 0.078
[2]
rs1805006
MC1R
16
89919510
C/A
A = 0.001
[2]
rs2228479
MC1R
16
89919532
G/A
A = 0.108
[2]
rs11547464
MC1R
16
89919683
G/A
A = 0.007
[2]
rs1805007
MC1R
16
89919709
C/T
T = 0.040
[2]
Y152OCH
MC1R
16
89919714
C/A
A = 0.000
[2]
rs1110400
MC1R
16
89919722
C/T
C = 0.026
[2]
rs1805008
MC1R
16
89919736
C/T
T = 0.068
[2]
rs885479
MC1R
16
89919746
C/T
T = 0.035
[2]
rs1805009
MC1R
16
89920138
G/C
C = 0.001
[2]
rs9904341
BIRC5
17
78214286
G/C
C = 0.343
[15]
rs7121
GNAS1
20
58903752
C/T
T = 0.465
[12]
rs132770
XRCC6
22
41621260
G/A
A = 0.153
[16]
∗MAF  =  minor allele frequency.


## 2.2. Genotyping

Polymorphisms were genotyped in 5 multiplex and 1 singleplex PCR reactions. SNPs in multiplexes were amplified using Qiagen Multiplex PCR kit (Qiagen, Hilden, Germany) and further analysed with single base extension (SBE) reactions using SNaPshot Multiplex Kit (Applied Biosystems, Foster City, CA). Multiplex 1 included 7 SNP positions (rs11894252, rs7579899, rs833061, rs3118523, rs748964, rs7105934, and rs4765623), multiplex 2 consisted of 6 SNPs (rs9679290, rs4953346, rs1049380, rs9904341, rs7121, and rs132770), and multiplex 3 comprised further 6 SNPs (rs34796867, rs113322875, rs34031609, rs113655247, rs17849897, and rs12105918). Sequences and concentrations of PCR and SBE primers are shown in Supplementary Tables  1 and  2, respectively in Supplementary Material available online athttp://dx.doi.org/10.1155/2015/860405. Amplifications were performed in 5 μL reaction volume consisting of 2.5 μL Qiagen Multiplex PCR mixture, 0.5 μL of Q solution, 0.5 μL of primer premix, and 1.5 μL of template DNA. The PCR reactions were carried out with the following temperature profile: 95°C/15 min, [94°/30 s, 58°C/90 s, 72°C/90 s] × 32, 72°C/10 min. The amplification products were cleaned with a mixture (1 : 1) of Exonuclease I (Exo I) and FastAP Thermosensitive Alkaline Phosphatase (FastAP) enzymes (Thermo Fisher Scientific Inc., Waltham, MA). The minisequencing reactions consisted of 0.5 μL of SNaPshot mix, 0.5 μL of SBE primer premix, 1 μL of purified PCR product, and 3 μL of DNase free water. The SBE temperature profile was as follows: [96°C/10 s, 50°C/5 s, 60°C/30 s] × 26. Products of the extension reactions were purified with 1 μL of FastAP Thermosensitive Alkaline Phosphatase enzyme and analyzed with ABI 3100 Avant Genetic Analyser (Applied Biosystems, Foster City, CA). The remaining 11 polymorphisms inMC1R (N29insA, rs1805005, rs1805006, rs2228479, rs11547464, rs1805007, Y152OCH, rs1110400, rs1805008, rs885479, and rs1805009) and 9 SNPs from theVDR gene (rs739837, rs731236, rs7975232, rs1544410, rs2228570, rs2238136, rs4516035, rs7139166, and rs11568820) were analysed in additional two multiplexes using protocols described elsewhere [22, 23]. The remaining insertion/deletion polymorphism inCASP8 was analysed in a single PCR reaction using HotStarTaq Master Mix Kit (Qiagen, Hilden, Germany) and following primer sequences 5′-6-FAM-AAACTTCTCCCATGGCCTCT-3′ and 5′-TATGAATGAGCCGAGGAAGG3′. Amplifications were carried out in 5 μL reaction volume consisting of 2.5 μL HotStarTaq Master Mix, 0.5 μL of primer premix (concentration of 0.125 μM), 0.5 μL of distilled water, and 1.5 μL of template DNA. The following temperature profile was applied: 95°C/15 min, [94°/30 s, 60°C/45 s, 72°C/60 s] × 30, 72°C/10 min. PCR products were directly analyzed by capillary electrophoresis using ABI 3100 Avant Genetic Analyser (Applied Biosystems, Foster City, CA).


## 2.3. Statistical Analyses

### 2.3.1. Population Analyses and Haplotype Reconstruction

Genetic data obtained for the 40 selected polymorphisms were tested for deviations from Hardy-Weinberg equilibrium. Hardy-Weinberg analysis was performed with Arlequin version 3.1 software (http://cmpg.unibe.ch/software/arlequin3/). Haploview version 4.2 software was used to test the degree of linkage disequilibrium among DNA variants inEPAS1,RXRA, andVDR (http://www.broadinstitute.org/). Haplotype reconstruction and frequency estimation were done for SNPs located inEPAS1,RXRA, andVDR using PHASE version 2.1 computer software (http://stephenslab.uchicago.edu/software.html#phase). Haplotype reconstruction and LD analysis did not include polymorphisms in theMC1R gene because they all may individually affect receptor performance [24] and polymorphisms inZC3H12a as besides one SNP (rs17849897) they turned out to be monomorphic in the tested population. In the case ofVDR gene, polymorphisms were divided into 3 blocks following the obtained results and literature suggestions [25, 26]. Positions rs739837, rs731236, rs7975232, and rs1544410 were included in the first block. The second block was formed by a single rs2228570 position whereas in the third block the remaining rs2238136, rs4516035, rs713966, and rs11568820 positions were considered. Minimal ORs detectable with the power of at least 80% were calculated with Power and Sample Size Program (PS Program) version 3.1.2 (http://biostat.mc.vanderbilt.edu/wiki/Main/PowerSampleSize).


### 2.3.2. 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) andP values were assigned. Additive, recessive, and dominant models of allele categorization were tested. In case of theMC1R gene, polymorphisms were divided in two groups according to their known influence on pigmentation phenotype. Strong variants marked with “R” (N29insA, rs1805006, rs1805007, rs11547464, rs1805008, rs1805009, and Y152OCH) have been considered as high-penetrance since they significantly diminish the receptor performance and have strong phenotypic effect (OR values of ~50 or higher). Weak variants marked with “r” are low-penetrance 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 “R” variants, three states were considered according to the existence of major function mutations (0 = no “R” variant carriers, 1 = one “R” variant carriers, and 2 = two “R” variant carriers) and the same approach was used in case of “r” 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 P 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/BenjaminiHochberg.xlsx). Finally, multivariate logistic regression was applied including simultaneous analysis of all the polymorphisms and haplotypes as well as disclosed interactions.


### 2.3.3. Epistasis Examination

Statistical locus-locus interactions were examined using two methods, that is, multifactor dimensionality reduction (MDR v2.0 beta 8.1http://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 P values were estimated for interactions disclosed with logistic regression.


### 2.3.4. Meta-Analysis

Meta-analysis for rs4765623 inSCARB1 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 P value was estimated based on 95% CI with the method described in Altman and Bland [31]. Heterogeneity was investigated using Cochrane Q P value and I
2 statistics which assessed the consistency in meta-analysis studies. I
2 ranges from 0% to 100% where 0% indicates no heterogeneity whereas increasing value of I
2 indicates increasing level of heterogeneity.


## 2.3.1. Population Analyses and Haplotype Reconstruction

Genetic data obtained for the 40 selected polymorphisms were tested for deviations from Hardy-Weinberg equilibrium. Hardy-Weinberg analysis was performed with Arlequin version 3.1 software (http://cmpg.unibe.ch/software/arlequin3/). Haploview version 4.2 software was used to test the degree of linkage disequilibrium among DNA variants inEPAS1,RXRA, andVDR (http://www.broadinstitute.org/). Haplotype reconstruction and frequency estimation were done for SNPs located inEPAS1,RXRA, andVDR using PHASE version 2.1 computer software (http://stephenslab.uchicago.edu/software.html#phase). Haplotype reconstruction and LD analysis did not include polymorphisms in theMC1R gene because they all may individually affect receptor performance [24] and polymorphisms inZC3H12a as besides one SNP (rs17849897) they turned out to be monomorphic in the tested population. In the case ofVDR gene, polymorphisms were divided into 3 blocks following the obtained results and literature suggestions [25, 26]. Positions rs739837, rs731236, rs7975232, and rs1544410 were included in the first block. The second block was formed by a single rs2228570 position whereas in the third block the remaining rs2238136, rs4516035, rs713966, and rs11568820 positions were considered. Minimal ORs detectable with the power of at least 80% were calculated with Power and Sample Size Program (PS Program) version 3.1.2 (http://biostat.mc.vanderbilt.edu/wiki/Main/PowerSampleSize).


## 2.3.2. 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) andP values were assigned. Additive, recessive, and dominant models of allele categorization were tested. In case of theMC1R gene, polymorphisms were divided in two groups according to their known influence on pigmentation phenotype. Strong variants marked with “R” (N29insA, rs1805006, rs1805007, rs11547464, rs1805008, rs1805009, and Y152OCH) have been considered as high-penetrance since they significantly diminish the receptor performance and have strong phenotypic effect (OR values of ~50 or higher). Weak variants marked with “r” are low-penetrance 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 “R” variants, three states were considered according to the existence of major function mutations (0 = no “R” variant carriers, 1 = one “R” variant carriers, and 2 = two “R” variant carriers) and the same approach was used in case of “r” 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 P 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/BenjaminiHochberg.xlsx). Finally, multivariate logistic regression was applied including simultaneous analysis of all the polymorphisms and haplotypes as well as disclosed interactions.


## 2.3.3. Epistasis Examination

Statistical locus-locus interactions were examined using two methods, that is, multifactor dimensionality reduction (MDR v2.0 beta 8.1http://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 P values were estimated for interactions disclosed with logistic regression.


## 2.3.4. Meta-Analysis

Meta-analysis for rs4765623 inSCARB1 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 P value was estimated based on 95% CI with the method described in Altman and Bland [31]. Heterogeneity was investigated using Cochrane Q P value and I
2 statistics which assessed the consistency in meta-analysis studies. I
2 ranges from 0% to 100% where 0% indicates no heterogeneity whereas increasing value of I
2 indicates increasing level of heterogeneity.


## 3. Results

### 3.1. Population Analyses and Haplotype Inference

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


### 3.2. 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 theSCARB1 gene to be associated with ccRCC in our population but only considering dominant mode of inheritance (P
=
0.016) and this result did not pass FDR procedure. According to logistic regression analysis individuals with at least one minor A allele have 1.7 higher odd of ccRCC development comparing to individuals with GG genotype (OR = 1.7, 95% CI: 1.104–2.582). Additive mode of allele inheritance did not reach the threshold of 0.05 significance (P
=
0.075). Except for rs4765623 inSCARB1 no other SNP showed main effect in our study. We also could not verify significance of 4 SNPs in theZC3H12a gene which were found to be monomorphic in this sample set. The only polymorphic position inZC3H12a (rs17849897) with the minor allele frequency of 1.4% turned out not to be associated with ccRCC. Results of univariate association analyses are provided in Supplementary Table 4.An association testing performed for the reconstructed haplotypes revealed that CAGT haplotype inVDR gene (comprising rs739837-rs731236-rs7975232-rs1544410) observed with frequency 0.68% might be a risk factor for ccRCC (P
=
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 theEPAS1 andRXRA 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 inSCARB1 andVDR, which confirmed suggestive association trends with ccRCC. According to MDR analysis, ccRCC is best explained by the four-factor model comprising rs4765623 inSCARB1, rs4953346 inEPAS1, rs9904341 inBIRC5, and rs7121 inGNAS1 with balanced accuracy (BA) of 0.5828, cross-validation consistency (CVC) of maximum level 10/10, and permutation testing P 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 inGNAS1, rs9904341 inBIRC5, and rs4765623 inSCARB1 (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 inSCARB1. This position removes 1.28% of entropy, what is understood to remove 1.28% of “uncertainty” in prediction of ccRCC status. The position rs7121 inGNAS1 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 inGNAS1 and rs9904341 inBIRC5 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).Table 2
Results of multifactor dimensionality reduction analysis performed for ccRCC status.
Multifactor dimensionality reduction analysis
Best candidate model (for up to five factor combinations)
BA∗
CVC∗
∗
P-value∗
∗
∗
rs4765623
0.5469
10/10
0.244
rs4765623; rs4953346
0.5125
4/10
0.658
rs4765623; rs4953346; rs9904341
0.5319
5/10
0.432
rs4765623; rs4953346; rs9904341; rs7121
0.5828
10/10
0.036
rs4765623; rs4953346; rs9904341; rs7121; rs132770
0.5598
10/10
0.140
Best model is marked with bold.
∗BA: balanced accuracy.
∗
∗CVC: cross-validation consistency.
∗
∗
∗Permutation testing P value.Figure 1
Interaction dendrogram provided with MDR analysis. Dendrogram interaction graphs are built using hierarchical cluster analysis and illustrate the presence, strength, and nature of epistatic effects. The more the line connecting two factors is moved to the right side, the stronger the interaction effect is. The red line between factors indicates a high degree of synergy whereas the orange line indicates a lesser degree of synergy.Figure 2
Entropy-based interaction radial graph provided with MDR analysis. Entropy values in the cells of individual factors indicate the main independent effects whereas the entropy values marked on the lines connecting two factors represent the effect of interaction.Pairwise interaction testing performed with logistic regression confirmed significance of rs4765623-rs7121 interaction (P
=
0.041) and revealed 5 additional epistatic effects (Table 3). The highest statistical significance was noted for interaction between rs4765623 inSCARB1 and rs9679290 inEPAS1 (P
=
0.008). Statistical significance was also noted for interaction between rs7121 inGNAS1 and rs9679290 inEPAS1 (P
=
0.016). The remaining interactions of potential influence on ccRCC includeGNAS1–MC1R (P
=
0.031),GNAS1–VDR (P
=
0.032), andMC1R-VDR (P
=
0.035). Overall, 6 genes (SCARB1,GNAS1,EPAS1,BIRC5,MC1R, andVDR) were found to be implicated in interaction effects of potential influence on ccRCC development (Table 3).Table 3
Interaction effects of statistical significance in ccRCC development revealed with binary logistic regression analysis.
Interacting loci
Genes
OR (95% CI) for interaction model
P value for interaction model
rs4765623∗rs9679290
SCARB1
∗
EPAS1
1.371 (1.087–1.729)
0.008
rs7121∗rs9679290
GNAS1
∗
EPAS1
1.284 (1.047–1.575)
0.016
rs7121∗rs1805008
GNAS1
∗
MC1R
2.601 (1.089–6.212)
0.031
rs7121∗rs2228570
GNAS1
∗
VDR
1.240 (1.018–1.510)
0.032
MC1R “r”∗rs7975232
MC1R
∗
VDR
1.341 (1.020–1.763)
0.035
rs7121∗rs4765623
GNAS1
∗
SCARB1
1.639 (1.020–2.639)
0.041Finally, 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 inVDR (P
=
0.005), interaction between rs4765623 inSCARB1 and rs9679290 inEPAS1 (P
=
0.007), interaction between variants “r” inMC1R and rs7975232 inVDR (P
=
0.010), and interaction between rs7121 inGNAS1 and rs2228570 inVDR (P
=
0.041) (Table 4). This final logistic regression model has the significance of χ
2  P
=
6.98
×
10
-
6 and explains 9.2% of the total risk in ccRCC development.Table 4
Results of multivariate binary logistic regression analysis concerning all SNPs under study, haplotypes, and disclosed interactions.
Multivariate binary logistic regression analysis
Variables in the model
Genes
OR (95% CI)
P value of variable in the model
Haplotype CAGT  (rs739837-rs731236-rs7975232-rs1544410)
VDR
cs
0.005
rs4765623∗rs9679290
SCARB1
∗
EPAS1
1.389 (1.094–1.763)
0.007
MC1R “r”∗rs7975232
MC1R
∗
VDR
1.455 (1.095–1.933)
0.010
rs7121∗rs2228570
GNAS1
∗
VDR
1.238 (1.009–1.518)
0.041


### 3.3. Results of Meta-Analysis

Meta-analysis of the pooled data from Purdue et al. [6] and the results of our study for rs4765623 inSCARB1 indicated that moderate heterogeneity may be present among the tested population samples with I
2
=
62.10% and Cochrane Q  P value = 0.048. Effect size of rs4765623 inSCARB1 calculated for the combined data was OR = 1.145, 95% CI: 1.093–1.200 with P
=
2.483
×
10
-
8.


## 3.1. Population Analyses and Haplotype Inference

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


## 3.2. 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 theSCARB1 gene to be associated with ccRCC in our population but only considering dominant mode of inheritance (P
=
0.016) and this result did not pass FDR procedure. According to logistic regression analysis individuals with at least one minor A allele have 1.7 higher odd of ccRCC development comparing to individuals with GG genotype (OR = 1.7, 95% CI: 1.104–2.582). Additive mode of allele inheritance did not reach the threshold of 0.05 significance (P
=
0.075). Except for rs4765623 inSCARB1 no other SNP showed main effect in our study. We also could not verify significance of 4 SNPs in theZC3H12a gene which were found to be monomorphic in this sample set. The only polymorphic position inZC3H12a (rs17849897) with the minor allele frequency of 1.4% turned out not to be associated with ccRCC. Results of univariate association analyses are provided in Supplementary Table 4.An association testing performed for the reconstructed haplotypes revealed that CAGT haplotype inVDR gene (comprising rs739837-rs731236-rs7975232-rs1544410) observed with frequency 0.68% might be a risk factor for ccRCC (P
=
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 theEPAS1 andRXRA 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 inSCARB1 andVDR, which confirmed suggestive association trends with ccRCC. According to MDR analysis, ccRCC is best explained by the four-factor model comprising rs4765623 inSCARB1, rs4953346 inEPAS1, rs9904341 inBIRC5, and rs7121 inGNAS1 with balanced accuracy (BA) of 0.5828, cross-validation consistency (CVC) of maximum level 10/10, and permutation testing P 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 inGNAS1, rs9904341 inBIRC5, and rs4765623 inSCARB1 (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 inSCARB1. This position removes 1.28% of entropy, what is understood to remove 1.28% of “uncertainty” in prediction of ccRCC status. The position rs7121 inGNAS1 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 inGNAS1 and rs9904341 inBIRC5 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).Table 2
Results of multifactor dimensionality reduction analysis performed for ccRCC status.
Multifactor dimensionality reduction analysis
Best candidate model (for up to five factor combinations)
BA∗
CVC∗
∗
P-value∗
∗
∗
rs4765623
0.5469
10/10
0.244
rs4765623; rs4953346
0.5125
4/10
0.658
rs4765623; rs4953346; rs9904341
0.5319
5/10
0.432
rs4765623; rs4953346; rs9904341; rs7121
0.5828
10/10
0.036
rs4765623; rs4953346; rs9904341; rs7121; rs132770
0.5598
10/10
0.140
Best model is marked with bold.
∗BA: balanced accuracy.
∗
∗CVC: cross-validation consistency.
∗
∗
∗Permutation testing P value.Figure 1
Interaction dendrogram provided with MDR analysis. Dendrogram interaction graphs are built using hierarchical cluster analysis and illustrate the presence, strength, and nature of epistatic effects. The more the line connecting two factors is moved to the right side, the stronger the interaction effect is. The red line between factors indicates a high degree of synergy whereas the orange line indicates a lesser degree of synergy.Figure 2
Entropy-based interaction radial graph provided with MDR analysis. Entropy values in the cells of individual factors indicate the main independent effects whereas the entropy values marked on the lines connecting two factors represent the effect of interaction.Pairwise interaction testing performed with logistic regression confirmed significance of rs4765623-rs7121 interaction (P
=
0.041) and revealed 5 additional epistatic effects (Table 3). The highest statistical significance was noted for interaction between rs4765623 inSCARB1 and rs9679290 inEPAS1 (P
=
0.008). Statistical significance was also noted for interaction between rs7121 inGNAS1 and rs9679290 inEPAS1 (P
=
0.016). The remaining interactions of potential influence on ccRCC includeGNAS1–MC1R (P
=
0.031),GNAS1–VDR (P
=
0.032), andMC1R-VDR (P
=
0.035). Overall, 6 genes (SCARB1,GNAS1,EPAS1,BIRC5,MC1R, andVDR) were found to be implicated in interaction effects of potential influence on ccRCC development (Table 3).Table 3
Interaction effects of statistical significance in ccRCC development revealed with binary logistic regression analysis.
Interacting loci
Genes
OR (95% CI) for interaction model
P value for interaction model
rs4765623∗rs9679290
SCARB1
∗
EPAS1
1.371 (1.087–1.729)
0.008
rs7121∗rs9679290
GNAS1
∗
EPAS1
1.284 (1.047–1.575)
0.016
rs7121∗rs1805008
GNAS1
∗
MC1R
2.601 (1.089–6.212)
0.031
rs7121∗rs2228570
GNAS1
∗
VDR
1.240 (1.018–1.510)
0.032
MC1R “r”∗rs7975232
MC1R
∗
VDR
1.341 (1.020–1.763)
0.035
rs7121∗rs4765623
GNAS1
∗
SCARB1
1.639 (1.020–2.639)
0.041Finally, 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 inVDR (P
=
0.005), interaction between rs4765623 inSCARB1 and rs9679290 inEPAS1 (P
=
0.007), interaction between variants “r” inMC1R and rs7975232 inVDR (P
=
0.010), and interaction between rs7121 inGNAS1 and rs2228570 inVDR (P
=
0.041) (Table 4). This final logistic regression model has the significance of χ
2  P
=
6.98
×
10
-
6 and explains 9.2% of the total risk in ccRCC development.Table 4
Results of multivariate binary logistic regression analysis concerning all SNPs under study, haplotypes, and disclosed interactions.
Multivariate binary logistic regression analysis
Variables in the model
Genes
OR (95% CI)
P value of variable in the model
Haplotype CAGT  (rs739837-rs731236-rs7975232-rs1544410)
VDR
cs
0.005
rs4765623∗rs9679290
SCARB1
∗
EPAS1
1.389 (1.094–1.763)
0.007
MC1R “r”∗rs7975232
MC1R
∗
VDR
1.455 (1.095–1.933)
0.010
rs7121∗rs2228570
GNAS1
∗
VDR
1.238 (1.009–1.518)
0.041


## 3.3. Results of Meta-Analysis

Meta-analysis of the pooled data from Purdue et al. [6] and the results of our study for rs4765623 inSCARB1 indicated that moderate heterogeneity may be present among the tested population samples with I
2
=
62.10% and Cochrane Q  P value = 0.048. Effect size of rs4765623 inSCARB1 calculated for the combined data was OR = 1.145, 95% CI: 1.093–1.200 with P
=
2.483
×
10
-
8.


## 4. 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–19].The study revealed rs4765623 in theSCARB1 gene to show association trend with ccRCC in the studied population sample. TheSCARB1 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 ofSCARB1 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 inSCARB1 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 (P
=
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 (f
=
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 inSCARB1 (P
=
0.03) among the three studied groups [6]. Meta-analysis including our data also suggests moderate interpopulation heterogeneity with ~62% I
2 value but the significance decreased from 0.03 to 0.048. The effect size of rs4765623 inSCARB1 for the combined data, including our population, was very similar to that obtained in Purdue et al. (OR = 1.145, P
=
2.483
×
10
-
8).Our study also revealed suggestive association trend ofVDR 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)2D3] 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 inVDR gene has been implicated in many carcinomas (e.g., [13, 42–47]). However, few reports are available on association betweenVDR 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 ofVDR 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 ofVDR 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 analysedVDR variation in additional 300 controls and did not observe this haplotype (P
=
0.008, data not shown). Association of a haplotype in theVDR gene with RCC was also noted by Karami et al. but in that study different region of theVDR 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 inSCARB1. 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, theSCARB1 locus was implicated by the MDR method in complex gene-gene 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 inGNAS1 and rs9904341 inBIRC5 and the second between rs7121 inGNAS1 and rs4765623 inSCARB1, best explains susceptibility to ccRCC. TheGNAS1 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]. TheBIRC5 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 inBIRC5 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 betweenGNAS1 andBIRC5 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 byBIRC5 gene (e.g., [66]). Therefore, synergistic interaction ofGNAS1 andBIRC5 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 apoptosisBIRC7 has been listed among genes differentially expressed in ccRCC [67]. The MDR method also revealed that rs7121 inGNAS1 is also implicated in synergistic interaction with rs4765623 inSCARB1. This interaction was further confirmed by logistic regression (P
=
0.041). Interestingly, logistic regression also revealed additional and more statistically significant epistatic effect withSCARB1 as a component, that is, interaction between rs4765623 inSCARB1 and rs9679290 inEPAS1 (P
=
0.008). Although a direct dependency betweenEPAS1 andSCARB1 or their protein products is unclear, it is worth noting that both these loci are related to the processes connected with angiogenesis. TheSCARB1 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 isEPAS1 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 byEPAS1) occurs, which is what results in induction of transcription of many genes encoding proteins involved in angiogenesis [70]. Inactivation ofVHL 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 protooncogenesTGFα andc-Met [7, 71, 72].TheGNAS1 gene was found in our study to be implicated in other epistatic effects withMC1R (P
=
0.031),EPAS1 (P
=
0.016), andVDR (P
=
0.032). TheMC1R 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 inMITFgene [77]. Interestingly, product ofMC1R gene acts through G-protein. Therefore, there is a direct molecular dependence betweenMC1R andGNAS1 products. Importantly, it has been demonstrated thatMC1R 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 theMC1R andVDR genes (P
=
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 theTYRP1 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 inSCARB1 and haplotype inVDR 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, andMC1R) 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 inVDR and the gene-gene interactionsSCARB1-EPAS1,GNAS1-VDR, andMC1R-VDR were included in the final logistic regression model explaining 9.2% of the total risk in ccRCC development.



---

*Source: 860405-2015-04-06.xml*