The Role of Polymorphisms in Vitamin D-Related Genes in Response to Vitamin D Supplementation

Background. Vitamin D deficiency represents a major healthcare problem. Vitamin D status is influenced by genetic and environmental determinants. Several observational studies have evaluated the association of single-nucleotide polymorphisms (SNPs) in vitamin D-related genes and vitamin D levels. Nevertheless, little is known about the role of these SNPs in the response to vitamin D supplementation. We conducted an interventional study to define the association between SNPs in vitamin D-related genes and the response to vitamin D supplementation in 100 self-reported healthy women of Arab ancestry for the majority. Methods. A total of 100 healthy female subjects received a weekly oral dose of 50,000 IU vitamin D for 12 weeks. Serum vitamin D concentration and metabolic profiles were measured at baseline and 12 weeks post-vitamin D supplementation. The genotypes of 37 SNPs selected from previously reported vitamin D-related genes have been assessed by Fluidigm genotyping assay. Results. Rs731236 (VDR gene) and rs7116978 (CYP2R1 gene) showed a significant association with vitamin D status. The rs731236 GG genotype and the rs7116978 CC genotype were associated with a “vitamin D sufficiency” state. Rs731236 GG and rs7116978 CC genotypes showed a higher response to vitamin D supplementation. Transcription factor binding site prediction analysis showed altered binding sites for transcription factors according to the different rs7116978 alleles. Interestingly, the 37 SNPs previously established to play a role in vitamin D-related pathways explained very little of the response to vitamin D supplementation in our cohort, suggesting the existence of alternative loci whose number and effect size need to be investigated in future studies. Conclusion. In this paper, we present novel data on vitamin D-related SNPs and response to vitamin D supplementation demonstrating the feasibility of applying functional genomic approaches in interventional studies to assess individual-level responses to vitamin D supplementation.


Introduction
Vitamin D plays an important role in the endocrine system, and it takes part in several biological processes such as blood pressure regulation, calcium and phosphate homeostasis, nerve conduction, skeletal development, erythropoiesis, and so forth. [1][2][3]. The active form of vitamin D, 1,25-dihydroxyvitamin D [1,25(OH)2D], regulates the expression of vitamin D-related genes involved in calcium transport and bone matrix protein [4,5]. Nevertheless, vitamin D deficiency has been well documented worldwide [6]. Several factors have been shown to contribute to vitamin D Nutrients 2020, 12, 2608 3 of 16 To classify the vitamin D levels S/D/I (sufficient, deficient, insufficient), pre-and postsupplementation, we followed the Endocrine Society criteria, where vitamin D deficiency is defined as serum 25(OH)D levels below 20 ng/mL (50 nmol/liter), vitamin D insufficiency is defined as 25(OH)D levels between 21 and 29 ng/mL (52.5-72.5 nmol/liter), and sufficiency is defined as serum 25(OH)D levels ≥30 ng/mL [41]. At the end of the intervention, participants were classified as either responder (R) to vitamin D supplementation (those who achieved serum levels of 25(OH) D above 20 ng/mL) or nonresponder (NR) (those whose serum levels of 25(OH) D remained <20 ng/mL) [19,42]. The biochemical measurements were performed at the Biomedical Labs in QU, while the molecular profiling was performed at Sidra Medicine.

Anthropometric Measures
The height was measured to the nearest 0.1 cm without shoes, using a Seca GmbH & Co. kg stadiometer. WC (waist circumference, cm) was measured midway between the lowest rib and the superior border of the iliac crest, just above the navel on standing subjects using an inelastic tape. Weight to the nearest 0.1 kg was measured in light clothing without shoes using a Seca GmbH & Co. kg weight scale. BMI (Body Mass Index) was calculated as weight (kg) divided by the squared height (m 2 ).
SNP-specific assays were designed and ordered through the Fluidigm D3 Assay Design tool (https://d3.fluidigm.com/account/login). SNP-related sequences are reported in Supplementary Table S1. The genotypes of the polymorphisms listed above were determined by polymerase chain reaction (PCR). In the first step, two preamplification primers (locus-specific primer (LSP) and specific target amplification (STA) primer) were used to amplify the target region containing the SNP. The DNA sequences of the SNPs of interest were preamplified simultaneously in one multiplex PCR, for each sample separately, on a Veriti Thermal Cycler (Applied Biosystems, Foster City, CA, USA), with the following conditions: hold at 95 • C for 15 min, 14 cycles at 95 • C for 15 s, and 60 • C for 4 min. A second amplification was performed on the Fluidigm 96.96 Dynamic Array.
Assay mixture was prepared by mixing 3 µL of each allele-specific primer (ASP), 8 µL of each locus-specific primer (LSP), and 29 µL DNA hydration buffer. One µl of each assay mix was combined with 2.5 µL of 2 × Assay Loading Reagent and 1.5 µL of nuclease-free water. Sample mixture was prepared by combining 2.5 µL of 1:50 diluted samples from the STA reaction and 3 µL Biotium Fast Probe Master Mix. The BioMark HD dynamic array was first primed with control line fluid, and then loaded with the samples and assay mixtures via the appropriate inlets using an IFC (integrated fluidic circuit) controller. The array chip was placed in the BioMark HD Instrument. PCR was carried out using the following cycling conditions: 50 • C for 2 min, 70 • C for 30 min, 25 • C for 10 min, and 95 • C for 5 min, followed by four touchdown cycles (95 • C for 15 s, from 64 • C to 61 • C for 45 s, 72 • C for 15 s) and 34 additional cycles (95 • C for 15 s, 60 • C for 45 s, 72 • C for 15 s). Each PCR reaction used distilled water instead of DNA as negative control. Results were plotted on a two-dimensional scatter plot of the major versus the minor allele using the BioMark SNP Genotyping Analysis software version 2.1.1.
Genotyping calls were assessed based on the allele discrimination plots and manually reviewed by looking at the single amplification plots. Genotyping calls were exported as a CSV file and processed for secondary analysis.

Data Analysis
Chi-square was used to test association and Hardy-Weinberg equilibrium (HWE) for each SNP. SNPs that did not pass the HWE criteria have been discarded from further analysis. All statistical assessments were two-sided and considered to be significant when p-value was <0.05. Pairwise comparisons and Bonferroni corrections were performed where relevant.
SNP genotypes were categorized as either major, intermediate (Inter), or minor based on the frequency. The genotype category for each SNP of each donor was used to build a color grid (pheatmap, R: A Language and Environment for Statistical Computing) [45].

Results
All characteristics of the participants are shown in Table 1. The mean age was 21 years. The mean BMI, height, waist circumference (WC), and weight were 24 kg/m 2 , 159 cm, 33 inches, and 62 kg, respectively. The participants had a mean level of ALT, AST, BUN, calcium, and creatinine in Phase 1 of 9.8 U/L, 15.1 U/L, 9.8 mg/dl, 9.2 mg/dl, and 0.46 mg/dL, respectively. The participants had a mean level of ALT, AST, BUN, calcium, and creatinine in Phase 2 of 12.9 U/L, 16.3 U/L, 12.7 mg/dL, 11.3 mg/dl, and 0.67 mg/dL, respectively. The mean level of the active form of vitamin D was 11.01 in Phase 1 and 34.3 in Phase 2. Table 1 also lists the 95% confidence interval (CI) of the mean and the p-value associated with the hypothesis of null distribution. The allele frequencies of the 37 SNPs assessed in this study are reported in Table 2. Four deviations from HWE were detected (SNP IDs: rs4516035, rs16846876, rs2298849, and rs2298850). The SNPs that did not meet HWE were removed from further analyses. One SNP (rs11234027) showed only one allele in the current cohort and was also removed from further analyses. We compared the allele distribution of the 37 SNPs with that of the populations included in the 1000 Genomes Project (1000G). Several SNPs showed differences in frequencies of the Ref (reference) allele in our cohort as compared to the data from the 1000G project. However, we cannot exclude that those differences are due to a small cohort size, or perhaps due to gender specificity.
Haplotype analysis was performed by chromosome (Chr.4 for GC gene, Chr. 11 for CYP2R1 and DHCR7 genes, Chr. 20 for CYP24A1 gene, Chr. 12 for CYP27B1 and VDR genes); no significant linkage disequilibrium was detected, probably due to the small cohort size (Supplementary Figure S1).  We also assessed the association of each SNP with the response to vitamin D supplementation (i.e., R/NR; responder, nonresponder) and pretreatment vitamin D status (i.e., S/D/I; sufficient, deficient, insufficient) classification. Out of all the SNPs assessed, rs731236 (VDR gene) and rs7116978 (CYP2R1 gene) showed a significant association with S/D/I classification (two-tailed chi-square test, p = 0.0336 and p = 0.0163, respectively, shown in bold in Table 3). Figure 1 shows the percentage of the genotypes of the two SNPs mentioned above and according to the S/D/I classification. The rs731236 GG genotype and the rs7116978 CC genotype were associated with a "vitamin D sufficiency" state. p = 0.0336 and p = 0.0163, respectively, shown in bold in Table 3). Figure 1 shows the percentage of the genotypes of the two SNPs mentioned above and according to the S/D/I classification. The rs731236 GG genotype and the rs7116978 CC genotype were associated with a "vitamin D sufficiency" state. To better understand the association of rs731236 and rs7116978 SNPs with the vitamin D sufficiency state, we performed pairwise comparisons by creating 3 × 2 contingency tables. For rs731236, only the contingency table of AA, AG, and GG genotypes versus I (Insufficient) and S (Sufficient) classification produced a significant p-value (chi-square = 8.868; p-value = 0.0119). For rs7116978, the two contingency tables of CC, CT, and TT genotypes versus D (Deficient) and S (Sufficient) classification and of CC, CT, and TT genotypes versus I (Insufficient) and S (Sufficient) classification produced a significant p-value (CC,CT,TT vs. D,S: chi-square = 7.272, p-value = 0.0264, CC,CT,TT vs. I,S: chi-square = 6.802, p-value = 0.0333). Nevertheless, after applying Bonferroni correction, only the pairwise comparison rs731236 AA, AG, and GG genotypes versus I and S classification remained significant (p-value = 0.0357). We cannot exclude this to be due to the small cohort size. To better understand the association of rs731236 and rs7116978 SNPs with the vitamin D sufficiency state, we performed pairwise comparisons by creating 3 × 2 contingency tables. For rs731236, only the contingency table of AA, AG, and GG genotypes versus I (Insufficient) and S (Sufficient) classification produced a significant p-value (chi-square = 8.868; p-value = 0.0119). For rs7116978, the two contingency tables of CC, CT, and TT genotypes versus D (Deficient) and S (Sufficient) classification and of CC, CT, and TT genotypes versus I (Insufficient) and S (Sufficient) classification produced a significant p-value (CC,CT,TT vs. D,S: chi-square = 7.272, p-value = 0.0264, CC,CT,TT vs. I,S: chi-square = 6.802, p-value = 0.0333). Nevertheless, after applying Bonferroni correction, only the pairwise comparison rs731236 AA, AG, and GG genotypes versus I and S classification remained significant (p-value = 0.0357). We cannot exclude this to be due to the small cohort size.
The increase of serum 25(OH) vitamin D level (delta vitamin D values, difference between Phase 1 and Phase 2) is also shown according to the different genotypes. For SNP rs731236, genotype GG showed a higher increase as compared to AG and AA genotypes. As SNP rs7116978 is intronic, we have assessed whether the different alleles may alter the binding sites for transcription factors. The results of these analyses are reported in Supplementary Figure S2. Figure 2 shows the allelic distribution of rs731236 and rs7116978 among the populations included in 1000G and the cohort included in the current study. For rs731236, the current study displayed a lower frequency of the G allele as compared to the African, European, South Asian, and American cohorts. Allelic frequencies of rs7116978 were more similar across the populations included in 1000G. However, the current study showed a higher frequency of rs7116978 T allele as compared to all the other populations assessed. We cannot exclude the possibility that these differences are related to gender specificity; those differences may also not be represented by a small cohort as the one employed in this study. We further sought to test whether the combination of different SNPs may predict R versus NR classification. Figure 3 shows heatmaps of SNPs categorized according to the genes they belong to. The heatmaps shown display the hierarchical clustering (complete-linkage clustering based on Euclidean distance) of SNPs and donors based on genotype. Horizontally, at the top of the heatmaps, we added color-coded labels that correspond to the genes. The hierarchical clustering of the donors (vertical hierarchical tree, left side of heatmaps) showed that groups formed based on the genotype of certain SNPs. No clear grouping of the participants based on their response to vitamin D supplementation was found in any of the genes assessed, suggesting that the genetic component alone does not predict responsiveness to vitamin D in our cohort. We further sought to test whether the combination of different SNPs may predict R versus NR classification. Figure 3 shows heatmaps of SNPs categorized according to the genes they belong to. The heatmaps shown display the hierarchical clustering (complete-linkage clustering based on Euclidean distance) of SNPs and donors based on genotype. Horizontally, at the top of the heatmaps, we added color-coded labels that correspond to the genes. The hierarchical clustering of the donors (vertical hierarchical tree, left side of heatmaps) showed that groups formed based on the genotype of certain SNPs. No clear grouping of the participants based on their response to vitamin D supplementation was found in any of the genes assessed, suggesting that the genetic component alone does not predict responsiveness to vitamin D in our cohort.

Discussion
Vitamin D is a fat-soluble vitamin that is produced when the skin is exposed to UVB radiation [46]. Vitamin D deficiency is associated with chronic liver [47] and kidney [48] diseases.
The actions of vitamin D are mediated by several proteins. VDR (vitamin D receptor) is a ligand-activated transcription factor that acts through vitamin D response elements located near the start sites of target genes to regulate gene expression [49]. GC (group-specific component) encodes for a vitamin D binding protein that plays an important role in the transport and metabolism of vitamin D, being the major plasma carrier for vitamin D and its metabolites [50]. The vitamin D pathway also involves a series of cytochrome P450-containing sterol hydroxylases that generate and degrade the active hormone serving as a ligand for the vitamin D receptor-mediated transcriptional gene expression [51]. Among these hydroxylases, CYP2R1 is the principal enzyme carrying the hydroxylation of vitamin D to 25-hydroxyvitamin D in the liver [52]; the 1a-hydroxylation of 25(OH)D in the kidney by CYP27B1 generates the fully active vitamin D metabolite [51,52]. Cytochrome P450 family 24 subfamily A member 1 (CYP24A1) gene encodes a mitochondrial monooxygenase which catalyzes the 24-hydroxylation of 1,25-dihydroxyvitamin D3 [51].
Recent studies have suggested that SNPs within the genes above may influence the level or activity of vitamin D, but studies exploring the association of SNPs in vitamin D-related genes with the response to supplementation of vitamin D are still lacking.
The present interventional study assessed the association between SNPs previously established to play a role in vitamin D-related genes and the responsiveness to vitamin D supplementation after the intervention. Our data showed that both rs731236 (VDR gene) and rs7116978 (CYP2R1 gene) have a significant association with the vitamin D status, where the rs731236 GG and the rs7116978 CC genotypes were associated with a "vitamin D sufficiency" state. In addition, rs731236 GG and rs7116978 CC genotypes were associated with a higher response to vitamin D supplementation. With the exception of rs731236 and rs7116978, the remaining SNPs previously established to play a role in vitamin D-related pathways explained very little the response to vitamin D supplementation in our cohort, suggesting the existence of alternative loci whose number and effect size warrant future studies.
Rs7116978 is located in an intronic region. Whereas the recognition of functional variants in the coding region is relatively simple, detecting changes in the noncoding region is more challenging due to the lack of a clear connection between nucleotide differences and regulatory functions. Some introns have been shown to contain transcription factor binding sites [53,54] and to influence gene expression by several known and unknown mechanisms such as intron-mediated enhancement [55,56]. In this study, we speculated that the two alleles of rs7116978 may change potential binding sites for different transcription factors. We have performed TFBS (transcription factor binding site) prediction analysis in silico to verify this hypothesis. The two alleles were predicted to modify transcription binding sites of several transcription factors, including transcription factor 4 (TCF4) and the human glucocorticoid receptor (GR) (Supplementary Figure S2). Several studies have shown the link between these transcription factors with vitamin D metabolism. Beildeck et al. found that 1,25(OH)2D3 induces the expression of TCF4 in several human cell lines via a VDR-dependent mechanism [57]. As per the proposed model, TCF4 upregulation would enhance the repression of β-catenin/TCF target genes. With regards to the human glucocorticoid receptor (GR), cross-sectional studies have shown that the chronic use of glucocorticoids is associated with low levels of 25(OH)D [58]. It appears that the glucocorticoids can upregulate the renal expression of CYP24A1 which in turn catabolizes 25(OH)D and 1,25(OH)2D to water-soluble inactive agents, thus mediating deficiency of vitamin D [58,59]. Another link between the use of glucocorticoids and vitamin D is that 1,25(OH)2D/VDR and glucocorticoids/glucocorticoid receptor (GR) intracellular signaling pathways cross-talk so that the increased levels of vitamin D may enhance the responsiveness of certain target cells to glucocorticoids [60][61][62]. Glucocorticoids with their cognate receptors translocate from the cell cytoplasm to the nucleus where they bind to glucocorticoid response element (GRE) to regulate gene transcription. As VDR and GR share some coactivators, VDR may promote individual gene transcription induced by GR; additionally, it has been reported that vitamin D may upregulate the binding of GR to GRE [60]. These may result in vitamin D enhancing the cellular responsiveness to glucocorticoids. Functional studies are warranted to understand the importance of rs7116978 polymorphism as a potential regulator of gene expression.
The vitamin D molecular pathway is a complex process that varies across individual genetic profiles and according to their health status. As this current topic is still in its infancy, additional studies are warranted to elucidate how genetic variation contributes to vitamin D supplementation outcomes.

Conclusions
In conclusion, the present interventional study shows novel data about the association of vitamin D-related SNPs with responsiveness to vitamin D supplementation in individuals of Arab ancestry. Our findings should be validated in additional studies employing a cohort of a bigger size.