Genome-wide association study of vitamin D concentrations and bone mineral density in the African American-Diabetes Heart Study

Relative to European Americans, African Americans have lower 25-hydroxyvitamin D (25OHD) and vitamin D binding protein (VDBP) concentrations, higher 1,25-dihydroxyvitamin D (1,25(OH)2D3) concentrations and bone mineral density (BMD), and paradoxically reduced burdens of calcified atherosclerotic plaque (subclinical atherosclerosis). To identify genetic factors contributing to vitamin D and BMD measures, association analysis of >14M variants was conducted in a maximum of 697 African American-Diabetes Heart Study participants with type 2 diabetes (T2D). The most significant association signals were detected for VDBP on chromosome 4; variants rs7041 (β = 0.44, SE = 0.019, P = 9.4x10-86) and rs4588 (β = 0.17, SE = 0.021, P = 3.5x10-08) in the group-specific component (vitamin D binding protein) gene (GC). These variants were found to be independently associated. In addition, rs7041 was also associated with bioavailable vitamin D (BAVD; β = 0.16, SE = 0.02, P = 3.3x10-19). Six rare variants were significantly associated with 25OHD, including a non-synonymous variant in HSPG2 (rs116788687; β = -1.07, SE = 0.17, P = 2.2x10-10) and an intronic variant in TNIK (rs143555701; β = -1.01, SE = 0.18, P = 9.0x10-10), both biologically related to bone development. Variants associated with 25OHD failed to replicate in African Americans from the Insulin Resistance Atherosclerosis Family Study (IRASFS). Evaluation of vitamin D metabolism and bone mineral density phenotypes in an African American population enriched for T2D could provide insight into ethnic specific differences in vitamin D metabolism and bone mineral density.


Introduction
T2D and osteoporosis are diseases of aging which contribute to increased fracture risk [1]. Related to vitamin D metabolism, serum concentrations of 25-hydroxyvitamin D (25OHD), 1,25-dihydroxyvitamin D (1,25(OH) 2 D 3 ), vitamin D binding protein (VDBP), and intact parathyroid hormone (iPTH) differ markedly between individuals of African and European ancestry [2][3][4][5][6]; however, bioavailable vitamin D (BAVD) levels appear to be similar. In concert with differences in vitamin D metabolism, African (versus European) ancestry is associated with higher bone mineral density (BMD) and lower prevalence of osteoporosis [7][8][9]. The effects of 25OHD, 1,25(OH) 2 D 3 and iPTH on their target organs/cells may also differ based on race/ethnicity. For example, skeletal resistance to the effects of iPTH is more pronounced in individuals who possess recent African ancestry [10].
Beyond population ancestry-based differences in bone mineral density and vitamin D metabolism, calcified atherosclerotic plaque [11] and calcium-containing kidney stones [12] also develop significantly less often in African Americans than European Americans. Inverse relationships exist between the development and progression of calcified atherosclerotic plaque (e.g., subclinical atherosclerosis) with bone mineralization; this relationship is independent from ethnicity [13]. Literature evidence supports higher risks for development of subclinical atherosclerosis and osteoporosis in populations of European ancestry [14].
The African American-Diabetes Heart Study (AA-DHS) was designed to investigate biologic relationships between vitamin D and its metabolites, subclinical atherosclerosis, and skeletal health. Unique aspects of this type 2 diabetes (T2D)-affected African American cohort is their generally preserved kidney function, low levels of albuminuria, and well-treated blood pressures, serum lipids, and glucose levels suggesting adequate access to healthcare. Herein, we report results of genome-wide association studies (GWAS) for concentrations of 25OHD, 1,25 (OH) 2 D 3 , VDBP, BAVD and iPTH and the traits of volumetric BMD (vBMD) determined by quantitative computed tomography (QCT) in the thoracic and lumbar spine in AA-DHS participants.

Study participants
AA-DHS includes African Americans with T2D recruited from two Wake Forest School of Medicine (WFSM) studies: the family-based Diabetes Heart Study (DHS) and unrelated individuals in the AA-DHS. DHS is a cross-sectional study of European American and African American families with siblings concordant for T2D. AA-DHS started after DHS and enrolled unrelated African Americans using identical inclusion criteria. AA-DHS objectives are to improve understanding of ethnic differences in coronary artery calcification and calcified plaque in populations of African and European ancestry. T2D was diagnosed in all participants after the age of 30 years in the absence of diabetic ketoacidosis and defined as fasting blood glucose �126 mg/dL or a random glucose �200 mg/dL, history of physician diagnosis of diabetes, or use of insulin or an oral hypoglycemic agent. The study was approved by the WFSM Institutional Review Board (BG05-033) and all participants provided written informed consent.

Phenotypes
25OHD was measured by Quest Diagnostics Nichols Institute (San Juan Capistrano, CA) and LabCorp using liquid chromatography and mass spectrometry and an immunochemiluminometric assay performed on the DiaSorin LIASON1 instrument, respectively, with a high paired sample concordance (n = 14, r 2 = 0.92). 1,25(OH) 2 D 3 was measured at Quest or Lab-Corp1 using liquid chromatography mass spectrometry and column chromatography, radioimmunoassay, respectively. iPTH was measured at LabCorp (Burlington, NC) using an electrochemiluminescence immunoassay. VDBP was measured in EDTA plasma samples that had been continuously stored at -80˚C without thawing since collection at baseline visits. Frozen plasma was thawed in a 37˚C water bath for 15 minutes, placed on ice, and then centrifuged at 1700 x g (2800 rpm) for 30 minutes at 4˚C. VDBP was determined using a Quantikine1 Human Vitamin D BP ELISA (cat. no. DDBP0; R&D Systems; Minneapolis, MN) according to manufacturer's instructions. Intraassay and inter-assay coefficients of variation were <10%. All assays were performed using a single lot of reagents and calibrators at WFSM. BAVD was computed as described by Powe et al. [6] using dissociation constants for the VDBP genotype variants as described by Arnaud and Constans [15].
Trabecular vertebral volumetric BMD (vBMD in mg/cm3) in the thoracic spine (T8-T11) and lumbar spine (T12-L3) were measured during computed tomography (CT) examination of the chest and abdomen using a standardized scanning protocol as previously described [16].

Genotyping
DNA from AA-DHS participants was extracted from peripheral blood using the PureGene system (Gentra Systems, Minneapolis, MN). The AA-DHS GWAS utilized the Illumina 5M chip. Quality control (QC) checks in AA-DHS led to the exclusion of 12 individuals from the analyses: six had call rates <95%, two had discordant self-reported and genetically determined sex, one had a heterozygosity score outside of the mean ±4 times the standard error interval, two had the same sample identifiers and one had 100% European ancestry. GWAS analyses were performed on up to 697 individuals. Imputation was performed using IMPUTE2 with phased haplotype data obtained from SHAPEIT2 [16]. Directly genotyped variants with a Hardy-Weinberg Equilibrium (HWE) pvalue �1x10 -04 and a call rate �95% were included and imputation was based on 3,436,913 autosomal variants. The multi-ethnic 1,000 Genomes Phase I integrated variant set release (v3) was used as the reference panel. Statistical analysis was performed for imputed variants that had an info score above 50%, a minor allele frequency (MAF) >1% and a HWE p-value �1x10 -04 .

Statistical analysis
Analyses were run using 25OHD, 1,25(OH) 2 D 3 , VDBP, BAVD, iPTH, thoracic vBMD, and lumbar vBMD as continuous traits. Phenotype correlations were explored using Spearman correlation coefficients. Phenotypes were transformed to best approximate a normal distribution, i.e. square root of 1,25(OH) 2 D 3 plus 1, thoracic vBMD, and lumbar vBMD; natural log of 25OHD plus 1 and iPTH; and log BAVD and VDBP. Age, sex, body mass index (BMI), kidney function based on the Chronic Kidney Disease-Epidemiology Collaboration (CKD-EPI) estimated glomerular filtration rate (eGFR), multivitamin use, and global African ancestry proportions were included as covariates. Global African ancestry proportion estimates were obtained by averaging the local ancestry estimates across the genome obtained from LAM-P-ANC and HAPMIX [17,18]. This approach utilized a linkage disequilibrium (LD) pruning algorithm that was applied with an r 2 threshold of 0.8 to select a subset of variants. Observed data at these variants were then combined with HapMap phase 3 genotypes obtained from Yoruba and CEPH samples; the HapMap samples were used as anchoring populations and were not included in the analysis. Results obtained from LAMP-ANC and HAPMIX were comparable, i.e. Spearman correlation estimates ranged between 0.88 and 0.97 for the 22 autosomes. For these continuous outcomes, linear mixed models were fitted using the Genomewide Efficient Mixed Model Analysis (GEMMA) software [10] and each variant was tested for association using an additive model. To account for inter-assay variability for 25OHD and 1,25(OH) 2 D 3, results from the individual platforms were combined via meta-analysis. Genome-wide significance was defined as a p-value <5.0x10 -8 . In addition, a principal components analysis (PCA) was used to determine the number of phenotypes needed to explain 95% of the observed variation. This number was then used to establish a more conservative threshold for significance. For statistically significant variants, a secondary analysis was performed with the inclusion of sex as a three level variable to control for menopause status. Conditional analysis for significantly associated variants was performed with inclusion of variant genotypes as covariates in the analysis, i.e. VDBP and BAVD. For variants rs7041 and rs4588, interaction was tested in a model which included both variants as main effects and a centered interaction term with adjustment for the same set of covariates as the main effect only models. In addition to the discovery analysis, previously reported variants associated with 25OHD [19], VDBP [20], iPTH [21] and lumbar vBMD [22] were examined in their respective AA-DHS GWAS.

Replication
Replication of associations observed for 25OHD and 1,25(OH) 2 D 3 was performed among Insulin Resistance Atherosclerosis Family Study (IRASFS) African Americans [23]. IRASFS is a population-based cohort in which to evaluate generalizability of findings from a diabetesenriched population. IRASFS African Americans were recruited from Los Angeles, CA with race/ethnicity determined by self-report. Informed consent was obtained from all participants.
Levels of 25OHD and 1,25(OH) 2 D 3 were measured as previously described [24]. Samples (n = 554) were genotyped on the Multi-Ethnic Genotyping Array (MEGA; Illumina, San Diego, CA) [25]. A total of 1,705,970 variants were successfully called for QC and analysis. Sample QC was performed to remove individuals with low call rates (<0.95), sex discordance, DNA contamination, or non-African ancestry. Duplicate samples were compared, and one of each duplicate pair was removed. Variants with missing position or alleles, allele mismatch, call rates <95%, departure from HWE (P�10 −04 ), frequency difference >0.2 compared with the 1000 Genomes Project phase 3 reference panel, and monomorphic variants were removed. Variants and samples that passed QC were used to perform pre-phasing with SHAPEIT2 [26] and imputation with IMPUTE2 [27] using a combined haplotype reference panel from the 1000 Genomes Project phase 3 (1000 Genomes Project Consortium 2010) and the African Genome Variation Project (AGVP) [28]. Variants with imputation info scores >0.4 were included in analysis. Tests of association between the seven genetic variants associated with 25OHD and 1,25(OH) 2 D 3 were computed using the Wald test from the variance component model implemented in Sequential Oligogenic Linkage Analysis Routines (SOLAR) [29]. Phenotypes were transformed to best approximate the distributional assumptions of conditional normality and homogeneity of variance, i.e. both were square root transformed. Admixture estimates were calculated using maximum likelihood estimation of individual ancestries as implemented in ADMIXTURE [30]. Genetic associations were calculated adjusting for age, sex, and admixture estimates. The primary inference was the additive model.
To model the joint effects of variants rs7041 and rs4588 on VDBP, a model that included both variants as main effects and a centered interaction term was examined. Mean VDBP levels by genotype combination at rs7041 and rs4588 are presented in S3 Table. The parameter estimate associated with the interaction term was nominally significant (-0.13 ± 0.06, P = 0.028), suggesting that the joint effect of risk variants at the two variants was not linear.
Similarly, five variants in or near the GC locus were associated with BAVD ( Table 2). These variants were moderately correlated (r 2 = 0.15-0.76) and ranged in association from P = 3.3x10 -19 to 2.1x10 -08 . Analysis conditioned on genotypes at variant rs7041 abolished association of additional variants (P = 0.013-0.85). Further examination of analyses conditional on the genotypes at variant rs7041 revealed only nominal evidence of association genome-wide (P>3.3x10 -07 ). Six variants were significantly associated with 25OHD (P = 2.2x10 -10 -1.4x10 -08 , Table 2). Adjustment for seasonal changes in 25OHD and dietary use of vitamin supplements did not significantly impact these observations. The most significant signal, rs116788687, was a nonsynonymous variant which resulted in the substitution of a methionine for isoleucine in the heparan sulfate proteoglycan 2 gene (HSPG2). This is a rare coding mutation (MAF = 0.0079%) that was observed in ten heterozygotes and one homozygote. Presence of the G allele resulted in a decrease of 25OHD (16.8 vs 20.5 ng/mL). Additionally, variant rs143555701 in the TRAF2 and NCK interacting kinase gene (TNIK) was associated with 25OHD (P = 9.0x10 -09 ). The minor allele of this rare variant was associated with a decrease in 25OHD levels. Moreover, variant rs117075918 was associated with 25OHD (P = 1.4x10 -08 ). This intronic variant had a MAF of 0.99% and the minor allele (G) was associated with a decrease in 25OHD levels. The remaining three variants associated with 25OHD were intergenic (rs116950775, P = 9.3x10 -09 and rs111955953, P = 1.4x10 -08 ) or in a hypothetical gene region (rs114001906, P = 1.3x10 -08 ).
Results from previously reported genome-wide significant loci associated with 25OHD [19], VDBP [20], iPTH [21] and lumbar vBMD [22] are presented in S4 Table. Although few studies have targeted recruitment of African Americans with T2D and measures of vitamin D, replication was attempted in the IRASFS African American cohort. In total, six variants were evaluated for replication with 25OHD and one variant was evaluated for replication with 1,25(OH) 2 D 3 . Consistent with the allele frequencies observed in AA-DHS, all variants were low-frequency, i.e. MAF<5%. Although no variant was identified with nominal evidence of association (P<0.05), a consistent direction of effect was observed for rs143555701, rs114001906, and rs117075918 with 25OHD and rs80068476 with 1,25(OH) 2 D 3 (S5 Table).

Discussion
A GWAS was performed for key clinical phenotypes targeting vitamin D metabolism and bone mineral density in an African American cohort with T2D. Ten independent genomewide significant loci were identified. Of these phenotypes, seven rare variants were significantly associated with 25OHD or 1,25(OH) 2 D 3 consistently resulting in a decrease of their circulating levels. In addition, a detailed analysis of VDBP not only confirmed the striking association of a non-synonymous mutation in GC, but also provided evidence for a second association signal at this locus.
The most significant association signal observed from this series of GWAS was association of a non-synonymous variant, rs7041, in exon 11 of GC with VDBP (P = 9.35x10 -86 ). GC encodes the VDBP which is the primary vitamin D carrier protein which regulates the bioavailability of 25OHD. In contrast to the seminal GWAS publications [20,31] performed in European-derived populations where the T allele frequency was 35%, in AA-DHS the frequency was 16%. Despite this difference, rs7041 was associated with a 4.5-fold reduction in VDBP, i.e. 271.9, 149.4, and 49.9 ug/mL for 0, 1, and 2 copies of rs7041-A, respectively, consistent with prior results in a European-derived population [20], but with approximately two fold lower levels of VDBP by genotype in this African American population with T2D. Despite the fact that vitamin D level is an important determinant of BMD [32], this variant was not significantly associated with the measures of vBMD assessed herein, i.e. lumbar vBMD, P = 0.12 or thoracic vBMD, P = 0.04. Conditional on association at rs7041, we observed a second independent signal of association at the GC locus. A second non-synonymous variant, rs4588, also located in GC was associated with VDBP (P = 1.4x10 -14 ). Notably, this variant was only nominally associated with VDBP (P = 0.031) prior to conditional analysis. This mutation was initially reported by immunophelometric methods for the quantification of the three most common VDBP isotypes [31]. In AA-DHS the minor allele rs4588-T (MAF = 11%) was associated with an increase in VDBP, i.e. 80.7, 81.9, and 100.2 ug/mL for 0, 1, and 2 copies, respectively. An interaction model including the effects of rs7041 and rs4588 was evaluated with VDBP as the outcome. Although not previously evaluated, we posited that variant lowering alleles would interact and result in reductions in VDBP. Upon formal analysis and consistent with the single variant results, the interaction model was associated with decreased VDBP levels (P = 0.028).
Consistent with its role as the primary vitamin D carrier protein, variants at the GC locus were also associated with BAVD in the AA-DHS cohort. The majority of 25OHD in circulation is bound to carrier protein, i.e. VDBP (85-90%) and albumin (10-15%) [33], leaving less than 1% of 25OHD free and available to cells (BAVD). The most significant variant associated with BAVD was the non-synonymous variant, rs7041 (P = 3.3x10 -19 ), associated with a decrease in BAVD, i.e. 6.8, 4.4, and 3.0 ng/mL for 0, 1, and 2 copies of rs7041-C, respectively. In total, five variants were significantly associated with BAVD at the GC locus; however, association was abolished after analyses conditioned on rs7041. Interestingly, variant rs4588 was not associated with BAVD (P = 0.68) or after analysis conditioned on rs7041 (P = 0.19). Taken together, these findings represent the first GWAS for BAVD in the African American population and compliment previously published findings in a European American population [20].
Six loci were significantly associated with 25OHD. The most significant observation was a non-synonymous variant, rs116788687 (P = 2.2x10 -10 ) in HSPG2. HSPG2 encodes a ubiquitous heparan sulfate proteoglycan, perlecan. Perlecan is a large multidomain proteoglycan involved in diverse biological functions. In knockout mouse models, perlecan has been identified as a major component of the pericellular material surrounding osteocyte processes where it functions to maintain structural integrity which is crucial for uninhibited interstitial fluid movement [34]. This variant was not associated in IRASFS African Americans (P = 0.92), although few had T2D (n = 76).
Intronic variant rs143555701 in TNIK, located on chromosome 3, was associated with decreased 25OHD levels (P = 9.0x10 -10 ). TNIK encodes a serine/threonine kinase that functions as an activator of the Wnt signaling pathway [35]. Both vitamin D, specifically 1,25 (OH) 2 D 3 , and the Wnt signaling pathway have been implicated in osteoblastic proliferation and differentiation [36]. Despite this relationship, this variant was not associated with variation in 1,25(OH) 2 D 3 levels (P = 0.53). Based on these observations, the TNIK locus could represent a cross talk pathway between vitamin D and the Wnt signaling pathway. A consistent direction of effect was observed in IRASFS, whereby increasing copies of the minor allele, rs143555701-T, were associated with decreases in 25OHD although this finding was nonsignificant (P = 0.45). Additional loci significantly associated with 25OHD were located intergenically (rs116950775, P = 9.3x10 -09 and rs111955953, P = 1.4x10 -08 ) or in genes of relatively unknown function (rs114001906, P = 1.3x10 -08 and rs117075918, P = 1.4x10 -08 ). While all associated variants were directly genotyped in AA-DHS, suggesting high quality genotype data, none were replicated in IRASFS.
The largest GWAS meta-analysis of circulating 25OHD conducted to date was from the Study of Underlying Genetic Determinants of Vitamin D and Highly Related Traits (SUNLIGHT) Consortium which included 33,996 European ancestry participants from 15 cohorts [37]. In total, four loci were significantly associated with 25OHD. The most significant signal observed was rs2282679 in the GC locus which was not associated with 25OHD in AA-DHS (P = 0.13). Notably, rs2282679 and rs4588 are in complete LD (D'/r 2 = 1.0) in European ancestry populations but not in African-derived populations (YRI, D' = 0.40, r 2 = 0.11), i.e. rs4588 but not rs2282679 was associated with 25OHD in AA-DHS. Index variants for the remaining three loci were not directly assessed in AA-DHS; however, regional examination of DHCR7 (P>0.17), CYP2R1 (P>0.098) and CYP24A1 (P>0.028) failed to identify signals of association. These results are not surprising given the large difference in sample size and differential recruitment strategies. Specifically, AA-DHS focused on individuals with recent African ancestry who, on average, have lower levels of 25OHD levels than European populations, i.e. 20.4 ng/mL in AA-DHS vs 27.92 ng/mL in SUN-LIGHT European-derived populations. Moreover, AA-DHS is enriched for T2D which is associated with altered vitamin D homeostasis, i.e. low 25OHD [38]. These observations support the potential for a differential genetic basis between studies.
Circulating vitamin D is hydroxylated in the kidney to form biologically active 1,25 (OH) 2 D 3 [39]. In contrast to the multiple genetic associations observed with 25OHD, only a single variant was significantly associated with 1,25(OH) 2 D 3 . Variant rs80068476 (P = 6.2x10 -09 ) was associated with decreased 1,25(OH) 2 D 3 levels. This variant is an intergenic variant located between microRNA 4532 (MIR4532), a non-coding RNA involved in post-transcriptional regulation of gene expression, and protein phosphatase 4, regulatory subunit 1-like (PPP4R1L). This variant was not associated in IRASFS African Americans (P = 0.35) but had a consistent direction of effect whereby increasing numbers of the minor allele rs80068476-T were associated with decreasing levels of 1,25(OH) 2 D 3 . In comparison with prior studies, only one GWAS for 1,25(OH) 2 D 3 has been published. In a Hispanic American sample from the IRASFS, a single variant, rs6680429, in the Dab, reelin signal transducer, homolog 1 gene (DAB1) was associated in a pilot study of 229 individuals with replication in the total cohort (n = 1190) [40]. This variant was not associated in the AA-DHS cohort (P = 0.64), despite being a relatively common genetic variant in AA-DHS African Americans (MAF = 0.31).
GWAS remains a powerful tool to understand the genetic architecture of complex biological phenotypes, but it is not without limitations. While the sample size associated with contemporary GWAS, e.g. type 2 diabetes, blood lipids, blood pressure, etc., are comparatively large, the vitamin D and BMD phenotypes examined here represent traits that are not routinely assessed. Moreover, the study design for the AA-DHS focuses exclusively on African Americans with T2D; thus, limiting our ability to identify an identically ascertained replication cohort. In addition, this is among the first genetics studies focused exclusively on the African American population where differences in bone mineral density and the vitamin D metabolism often oppose those in European-derived populations. Despite identification of coding variants related to vitamin D metabolism, additional loci located in intronic and intergenic regions challenges our ability to make causal inferences. This suggests potential contributions from variants not currently assessed, despite use of a dense genotyping array and imputation to 1000 Genomes.
In conclusion, a GWAS of phenotypes focused on vitamin D metabolism and bone mineral density was conducted in the AA-DHS cohort. Recruited on the basis of T2D, the AA-DHS cohort may be enriched for novel genomic loci as insulin resistance is associated with decreased vitamin D concentrations [41,42]. Indeed, eight novel loci were identified for their contribution to 25OHD, 1,25(OH) 2 D 3 , and BAVD while two independent loci were replicated for association with VDBP. This study represents the first GWAS to explore the broad range of vitamin D metabolites. As such, additional studies are required to replicate associated loci. With identification of putative functional variants in coding regions, functional studies are needed to validate their role and provide insight into physiology.

Conclusions
Ethnic-specific differences in bone mineral density and vitamin D metabolism phenotypes could provide insight into the developmental mechanisms of subclinical atherosclerosis and osteoporosis.