Genome-wide association study identifies 143 loci associated with 25 hydroxyvitamin D concentration

Abstract

Vitamin D deficiency is a candidate risk factor for a range of adverse health outcomes. In a genome-wide association study of 25 hydroxyvitamin D (25OHD) concentration in 417,580 Europeans we identify 143 independent loci in 112 1-Mb regions, providing insights into the physiology of vitamin D and implicating genes involved in lipid and lipoprotein metabolism, dermal tissue properties, and the sulphonation and glucuronidation of 25OHD. Mendelian randomization models find no robust evidence that 25OHD concentration has causal effects on candidate phenotypes (e.g. BMI, psychiatric disorders), but many phenotypes have (direct or indirect) causal effects on 25OHD concentration, clarifying the epidemiological relationship between 25OHD status and the health outcomes examined in this study.

Fig. 1. Outline of key analytic steps described in this study.

BMI body mass index, eQTL expression quantitative trait locus, FUMA functional mapping and annotation, GREML genomic relationship restricted maximum likelihood, GSMR generalised summary-based MR, GWAS genome-wide association study, GxE genotype-by-environment interaction, LDSC linkage disequilibrium score regression, MR Mendelian Randomisation, SMR summary-based MR, QC quality control, UKB UK Biobank, vQTL variance quantitative trait locus.

Fig. 2. Heritability, SNP-based heritability and variance explained in out-of-sample prediction.

Heritability (left panel) and SNP-based heritability (middle panel) estimates and the variance explained in out-of-sample prediction (right panel). Heritability and SNP-based heritability estimates are presented with 95% confidence interval. GCTA-GREML was used to estimate heritability from a UKB subset that included all pairs of individuals related with coefficient of relationship > 0.2 (N = 58,738 relatives). GCTA-GREML was used to estimate SNP-based heritabilities labelled GREML summer or winter using samples of ~50 K participants randomly drawn from the UKB. The SBayesR SNP-based heritability is estimated from the GWAS summary statistics (N = 417,580). In out-of-sample prediction into the QIMR and the UKB replication (UKBR) samples, polygenic risk scores (PRS) calculated by the standard P-value threshold method (P + T) were outperformed by using SNP effect estimates calculated from GWAS summary statistics using the SBayesR or SBayesS methods. Bars of the same colour used the same methodology (noting that SBayesR generates an estimate of SNP-based heritability as well as SNP effect sizes in prediction analysis). The numbers on top of the bars are −log10 P-value of the regression of 25OHD on 25OHD PRS. COJO conditional and joint, rg genetic correlation, s.e. standard error.

Fig. 3. Manhattan plot of the 25OHD GWAS in the UK Biobank.

Manhattan plot showing the −log10 P-values from the fastGWA association test of 25-hydroxyvitamin D (25OHD) with genome-wide SNPs. Red dots represent independent variants identified as genome-wide significant with conditional and joint analysis (COJO) applied to the GWAS summary statistics. The horizontal axis shows each chromosome, with 23 representing the X chromosome. The vertical axis is restricted to −log10 P-values < 150. Five COJO SNPs were associated at P < 1 × 10⁻¹⁵⁰ (four on chromosome 11 and one on chromosome 4; Supplementary Data 3) and have approximate locations represented by three red triangles at the top edge of the plot.

Fig. 4. Summary of selected variants associated with 25OHD in the UK Biobank.

Top panel shows loci associated with skin integrity, lipid and lipoprotein pathways and CYP450 and steroid-associated enzymes. Lower panel shows selected variants and putative mechanisms related to 25-hydroxyvitamin D (25OHD) concentration. For selected inter-genic loci, the nearest (upstream or downstream) gene is shown in brackets. The distance between the loci and the nearest gene is shown in base pairs. *CYP24A1 was also the closest gene for an additional three inter-genic loci with distances between 32,865–55,282 base pairs. ns non-synonymous variant, x2 or x3 two or three loci found within the gene. We note that nearest gene annotation should be interpreted recognising that these are not proof of a causal relationship between the associated SNP and expression of the gene (none-the-less supporting evidence for causal relationships is given by the SMR analyses for some loci).

Fig. 5. Mendelian randomisation analysis estimates from GWAS results from 25OHD and selected phenotypes.

Bidirectional Generalised Summary-data level Mendelian Randomisation (GSMR) between 25-hydroxyvitamin D (25OHD) concentrations and selected phenotypes, by three types of adjustments for body mass index (BMI) and with/without HEIDI filtering of pleiotropic loci. Panel a shows the estimate of the causal effect estimate (${\widehat{b}}_{xy}$; dots), and 95% confidence intervals (bars), of 25OHD concentration on selected phenotypes. ${\widehat{b}}_{xy}$ (and respective P-values) was obtained with the generalised summary-data-based Mendelian randomisation (GSMR) method. Negative ${\widehat{b}}_{xy}$ indicates that variants associated with increased 25OHD concentration were associated with a smaller value/reduced risk for the phenotypes of interest. Panel b shows the estimate of the causal effect (and 95% confidence intervals) of the same selected phenotypes on 25OHD concentration. Results are presented with (upper half) and without (lower half) filtering of pleiotropic associations with the heterogeneity in dependent instruments (HEIDI) test, respectively. The numbers above each effect size indicate the number of SNP instruments used in each analysis. For each set of analyses, we show GWAS results: (i) without adjustment for body mass index (BMI), (ii) with BMI included as a covariate and (iii) conditioned on BMI using mtCOJO. The GMSR estimates and 95% confidence intervals are shown in three colours according to the P-value thresholds: grey (non-significant; PGSMR > 0.05), yellow (nominally significant; PGSMR < 0.05) and red (significant after Bonferonni correction for multiple testing; PGSMR < 0.003). All ${\widehat{b}}_{xy}$, respective standard errors and P-values can be found in tabular form in Supplementary Data 11.