Vitamin D may protect against poor health outcomes, including heart disease, diabetes, certain cancers, and overall mortality [1–5]. Its biological properties include regulation of cell proliferation and immune function, as well as increased cell differentiation and apoptosis [6–10]. These mechanisms are controlled by the active metabolite 1,25-dihydroxyvitamin D (1,25(OH)₂D) and the vitamin D receptor (VDR), often in conjunction with retinoid X receptor alpha (RXRA) [11]. This 1,25(OH)₂D-VDR-RXRA complex binds to vitamin D response elements that can activate or repress gene transcription [12].

Circulating vitamin D levels could affect DNA methylation via transcriptional regulation or other mechanisms [13]. In mammals, DNA methylation is an epigenetic process by which a methyl group is transferred onto the C5 position of a cytosine, forming 5-methylcytosine. Increased methylation at CpG sites in promoter regions is associated with gene inactivation and transcriptional repression, while increased methylation at CpGs in gene bodies is associated with actively transcribed genes [14, 15]. Examples of other environmental exposures associated with methylation changes include smoking (for both smokers [16] and their offspring [17–19]), as well as body mass index (BMI) [20, 21], alcohol consumption [22], and nutrients such as folate, vitamin B12, and retinoic acid [23–27].

Some empirical evidence supports a link between vitamin D exposure and DNA methylation. Candidate gene approaches have observed that vitamin D is associated with methylation of CYP24A1 [28, 29], BMP2 [30], PTEN [31], and DKK1 [32]. Additionally, one epigenome-wide association study (EWAS) conducted among adolescent African-American males identified two sites (cg16317961 (MAPRE2) and cg04623955 (DIO3)) that were significantly associated with serum levels of the stable precursor to 1,25(OH)₂D, 25-hydroxyvitamin D (25(OH)D) [33]. However, those findings did not replicate in a subsequent EWAS conducted among Caucasian men, nor did that subsequent EWAS identify any novel associations [34]. Another EWAS observed no noteworthy associations between maternal 25(OH)D levels and methylation in cord blood [35], and an epigenome-wide in-vitro study identified no detectable methylation changes in blood mononuclear cells treated with vitamin D [36]. Several studies of the association between vitamin D and LINE-1 global methylation levels have also been negative [37–39].

To further investigate a possible link between vitamin D and DNA methylation, we studied the relationship between serum 25(OH)D and CpGs in or near seven vitamin D-related genes (VDR, RXRA, CYP2R1, CYP24A1, GC, CYP27B1, and DHCR7/NADSYN1) using a random sample of women from a large prospective cohort (“subcohort”). Based on our previous finding that serum 25(OH)D was associated with a 21% reduction in the hazard of breast cancer over 5 years of follow-up [3], and other research observations that methylation status can modify the responses of individuals to vitamin D treatment [28, 29], we also examined 25(OH)D-methylation interactions in relation to breast cancer risk. We additionally conducted an EWAS of serum 25(OH)D.

Women who developed breast cancer during the 5-year follow-up period were slightly older than those in the subcohort (58.7 years versus 55.7 years) and had lower prediagnosis 25(OH)D levels (32.3 versus 32.7 ng/mL). Cases were more likely to have more than one first-degree relative with breast cancer, to be postmenopausal, to be obese, or to be currently taking hormone therapy (Additional file 1: Table S1).

25(OH)D and methylation of vitamin D-related genes in the subcohort

Of the 198 CpGs from vitamin D-related genes, cg21201924 (RXRA) had the lowest p value for association with 25(OH)D in the subcohort (p = 0.0004; Table 1 and Additional file 1: Table S2). Twenty-two other candidate CpGs were significantly associated with 25(OH)D, all but one of which were located in RXRA, NADSYN1/DHCR7, or GC. The large overall contrast between our results and those expected by chance is illustrated by a quantile-quantile plot (Fig. 1a).

Rank	CpG	Gene / location type	Chromosome: position	Mean methylation level (SD)	Association with 25(OH)Da
					β	p value
1	cg21201924	RXRA / body	9: 137251825	0.76 (0.042)	−0.020	0.0004
2	cg02127980	RXRA / body	9: 137252116	0.40 (0.067)	−0.015	0.0004
3	cg17559402b	NADSYN1 / body	11: 71187890	0.97 (0.006)	−0.017	0.003
4	cg02059519	RXRA / body	9: 137250935	0.83 (0.026)	−0.012	0.003
5	cg09997530	GC / body	4: 72636217	0.91 (0.017)	0.015	0.005
6	cg04329455b	RXRA	9: 137215364	0.96 (0.008)	0.014	0.007
7	cg00268518b	NADSYN1 / TSS200	11: 71164106	0.01 (0.001)	0.011	0.008
8	cg03146219	NADSYN1 / body	11: 71189514	0.47 (0.097)	−0.012	0.009
9	cg13510651b	RXRA / body	9: 137227772	0.94 (0.008)	−0.010	0.01
10	cg03490288b	DHCR7 / body	11: 71146658	0.97 (0.006)	0.015	0.01
11	cg05785753	NADSYN1 / body	11: 71189490	0.59 (0.074)	−0.010	0.01
12	cg13687497	RXRA / body	9: 137249839	0.80 (0.023)	−0.010	0.01
13	cg07793224b	NADSYN1 / body	11: 71183180	0.97 (0.008)	0.017	0.02
14	cg26044621b	DHCR7 / 3´ UTR	11: 71145665	0.94 (0.011)	0.012	002
15	cg04837494	GC / 3´ UTR	4: 72608149	0.85 (0.033)	0.013	0.03
16	cg14236758	RXRA / body	9: 137252129	0.48 (0.066)	−0.009	0.03
17	cg04774822	NADSYN1 / body	11: 71165839	0.78 (0.033)	−0.009	0.03
18	cg16151558b	DHCR7 / TSS1500	11: 71159853	0.02 (0.066)	0.013	0.04
19	cg20372759b	CYP27B1 / TSS1500	12: 58162287	0.97 (0.006)	−0.010	0.04
20	cg24806812	GC / body	4: 72635202	0.93 (0.018)	0.012	0.04
21	cg14154547b	RXRA / body	9: 137293309	0.92 (0.010)	−0.007	0.04
22	cg07099121b	DHCR7 / 3´ UTR	11: 71146096	0.98 (0.004)	−0.010	004
23	cg16910670b	NADSYN1	11: 71215361	0.97 (0.005)	−0.011	0.05

TSS200 within 200 basepairs upstream of the transcription start site, TSS1500 within 1500 basepairs upstream of the transcription start site, UTR untranslated region

^aEstimated change in methylation (logit(β)) per 10 ng/mL change in serum 25-hydroxyvitamin D (25(OH)D)

^bIntraclass correlation coefficient < 0.5

Epigenome-wide association study of 25(OH)D in subcohort or cases

Within the subcohort, 25(OH)D was associated with methylation levels at three CpGs at q < 0.05 (Fig. 2a and Table 3). The CpG with the smallest p value was cg24350360 (EPHX1; p = 3.4 × 10⁻⁸), followed by cg06177555 (SPN; p = 9.8 × 10⁻⁸) and cg1324316 (SMARCD2; p = 2.9 × 10⁻⁷). Two other CpGs had q < 0.10: cg23761815 (SLC29A3; p = 5.1 × 10⁻⁷) and cg10401362 (DNAJB6; p = 1.0 × 10⁻⁶). The quantile-quantile plot (Fig. 2b) demonstrates that the observed p values systematically deviated from what was expected under the null hypothesis. For all five CpGs with q < 0.10, increases in serum 25(OH)D were associated with decreased methylation (Table 3; Additional file 1: Figure S2). All except cg24350360 (EPHX1) had ICCs > 0.5.

CpG site	Location (Chr:position)	Location type	Gene	Effect estimatea (95% CI)	p value	q value
cg24350360	1:225997662b	TSS200	EPHX1	−0.04 (−0.06 to −0.03)	3.4 × 10−8	0.01
cg06177555	16:29678624	3′ UTR	SPN	−0.02 (−0.03 to −0.01)	9.8 × 10−8	0.02
cg13243168	17:61915833	Body	SMARCD2	−0.02 (−0.03 to −0.01)	2.9 × 10−7	0.04
cg23761815	10:73083123	Body	SLC29A3	−0.02 (−0.03 to −0.01)	5.1 × 10−7	0.05
cg10401362	7:157185402	Body	DNAJB6	−0.02 (−0.03 to −0.01)	1.0 × 10−6	0.09

CI confidence interval, TSS200 within 200 basepairs upstream of the transcription start site, UTR untranslated region

^aEstimated change in methylation (logit(β)) per 10 ng/mL change in serum 25-hydroxyvitamin D (25(OH)D)

^bIntraclass correlation coefficient < 0.5

No CpGs were associated with 25(OH)D in case-only analyses (Additional file 1: Figure S3). When we compared the results of 25(OH)D-methylation association tests for the subcohort versus breast cancer cases (Fig. 3), no CpGs had a combined p < 1.2 × 10⁻⁷, the Bonferroni-corrected cut-point for significance. Sixteen CpGs with combined p values < 1.0 × 10⁻⁵ were deemed worthy of further investigation; all but two of which had ICCs > 0.5 (Table 4). Of the sixteen, nine had RHR p values < 0.05. Three of the latter nine were associated with 25(OH)D at q < 0.10 in the initial EWAS: cg13243168 (SMARCD2), cg23761815 (SLC29A3), and cg24350360 (EPHX1). Most of the CpGs with small Fisher combined p values were inversely associated with 25(OH)D in both cases and the subcohort.

Rank	CpG	Gene/ Location	25(OH)D-methylation association, subcohort		25(OH)D-methylation association, cases		Fisher combined p value	HR (95% CI) for methylation-breast cancer association	RHRs (95% CI)c	Interaction p value
			Cofficientb	p value	Cofficientb	p value
1	cg08092930	PPFIA1	−0.03	1.3 × 10−5	0.02	0.04	8.5 × 10−6*	1.01 (0.98–1.05)	1.15 (1.06–1.24)	6.3 × 10−4
2	cg23761815	SLC29A3	−0.02	5.1 × 10−7	− 0.01	0.03	2.7 × 10−7	1.13 (1.07–1.18)	1.22 (1.08–1.38)	0.002
3	cg13243168	SMARCD2	−0.02	2.9 × 10−7	− 7 × 10−4	0.87	4.0 × 10−6	1.05 (0.98–0.91)	1.29 (1.09–1.51)	0.002
4	cg15544721	PPP1R9A	−0.01	0.09	−0.03	6.4 × 10−6	8.8 × 10−6	0.98 (0.94–1.03)	0.87 (0.79–0.96)	0.008
5	cg11568290	5p15.1	−0.02	1.4 × 10−4	−0.01	0.004	7.7 × 10−6	1.15 (1.08–1.22)	1.23 (1.05–1.44)	0.009
6	cg19420720	P4HB	−0.01	0.002	0.02	2.3 × 10−4	8.1 × 10−6*	1.01 (0.94–1.08)	1.27 (1.06–1.52)	0.01
7	cg23839180	FAM49A	−0.02	0.05	−0.05	3.0 × 10−6	2.3 × 10−6	0.97 (0.94–1.03)	0.92 (0.86–0.98)	0.01
8	cg15320474d	UBD	0.02	2.8 × 10−6	− 0.01	0.16	7.0 × 10−6*	0.98 (0.83–1.04)	0.87 (0.77–0.99)	0.03
9	cg24350360d	EPHX1	−0.04	3.4 × 10−8	0.01	0.48	3.1 × 10−7*	1.02 (0.99–1.05)	1.08 (1.01–1.16)	0.04
10	cg22488164	PLBD1	−0.03	1.5 × 10−4	−0.03	5.0 × 10−4	1.3 × 10−6	1.06 (1.02–1.10)	0.93 (0.85–1.01)	0.07
11	cg10401362	DNAJB6	−0.02	1.0 × 10−6	−0.01	0.27	4.4 × 10−6	1.02 (0.97–1.07)	1.10 (0.98–1.24)	0.10
12	cg06177555	SPN	−0.02	9.8 × 10−8	−0.003	0.54	9.3 × 10−7	0.97 (0.91–1.03)	1.12 (0.98–1.28)	0.10
13	cg11277126	TRPC4AP	−0.02	7.7 × 10−5	−0.02	4.8 × 10−4	6.7 × 10−7	1.04 (0.98–1.11)	1.07 (0.92–1.24)	0.39
14	cg21527411	GLYAT	0.02	2.1 × 10−6	−0.01	0.30	9.6 × 10−6*	1.02 (0.97–1.07)	0.95 (0.83–1.08)	0.40
15	cg09914444	DMBX1	0.02	5.4 × 10−6	−0.01	0.08	6.9 × 10−6*	1.05 (1.01–1.10)	0.96 (0.87–1.06)	0.44
16	cg23999318	HIPK2	−0.02	2.8 × 10−4	−0.02	3.5 × 10−4	1.7 × 10−6	1.00 (0.94–1.05)	1.00 (0.88–1.14)	1.00

CI confidence interval, HR hazard ratio, RHR ratio of hazard ratio

^aAfter excluding those with missing covariate information

^bEstimated change in methylation (logit(β)) per 10 ng/mL change in serum 25-hydroxyvitamin D (25(OH)D)

^cChange in the methylation-breast cancer association for being in the 4th quartile of 25(OH)D (> 38.0 ng/mL) versus the first three (≤ 38.0 ng/mL)—a value > 1.00 indicates that the estimated HR for the methylation-breast cancer association is higher among those with higher 25(OH)D levels; similarly, an RHR < 1.00 indicates that the estimated HR for the methylation-breast cancer association is higher among those with lower 25(OH)D levels

^dIntraclass correlation coefficient < 0.5

*Effect estimate going in opposing direction for subcohort versus cases

Among our a priori candidate loci, we found that methylation levels at CpGs in or near RXRA, NADSYN1/DHCR7, and GC were associated with serum 25(OH)D levels. In our larger EWAS analysis, CpGs in EPHX1, SPN, and SMARCD2 had epigenome-wide significant associations with serum 25(OH)D. To our knowledge, we are the first to report a link between these three genes and 25(OH)D and the first to study methylation-25(OH)D interactions in relation to breast cancer risk.

For candidate CpG analyses, the top hit for both the methylation-25(OH)D association analysis and the breast cancer interaction analyses was cg21201924, located in the gene body of RXRA. As previously noted, RXRA acts as a transcription factor with 1,25(OH)₂D and VDR and changes to the gene’s expression and methylation levels could have widespread biological impacts. Changes in expression or methylation of GC, NADSYN1/DHCR7, and the other candidate genes may have less pervasive biological effects, but our findings support the hypothesis that these vitamin D-related genes and proteins may interact with circulating vitamin D levels to influence breast cancer risk.

Of the CpGs in or near vitamin D-related genes, most of those that were either associated with 25(OH)D or that showed evidence of interacting with 25(OH)D to affect breast cancer risk were located within gene bodies. Higher 25(OH)D levels tended to be associated with higher methylation, but there was no clear pattern linking CpG locations to the direction of the RHR in the interaction analysis.

None of the candidate CpGs from the vitamin D-related genes met the stringent criterion for statistical significance in the EWAS analysis, and thus there was no overlap between the genes identified in our EWAS and those reported previously to be associated with serum vitamin D levels. One of the two hits reported by Zhu et al. [33] (cg04623955 near DIO3) was also assessed in our sample, but we found no evidence of an association (p = 0.78). Other CpGs identified in their sample also failed to replicate, including cg23492043 (p = 0.64), cg00864867 (p = 0.15), and cg16826718 (p = 0.62). Eight CpGs reported by Florath et al. [34] were also assessed in our sample, but none were significantly associated with 25(OH)D (p values 0.09–0.85). These discrepancies could be related to differences in race, sex, or study design, but could also be the result of sampling variation.

As previously noted, vitamin D plays a role in immune response, including regulation of innate and adaptive immunity [9, 10], as well as detoxification [51]. Possible mechanisms for these actions could be through the 1,25(OH)₂D-VDR-RXRA complex and its effects on gene transcription [12]. Although there is no established link between 25(OH)D or vitamin D metabolism and SPN, SMARCD2, SLC29A3, or DNAJB6 specifically, the observed associations between these genes and 25(OH)D could be related to VDR or other components of immune function.

EPHX1 encodes epoxide hydrolase, an enzyme that breaks down epoxides from xenobiotic aromatic compounds (e.g., polycyclic aromatic hydrocarbons, benzene) [52]. Further, EPHX1 regulation of detoxification via CYP450 enzymes has been shown to modulate the immune response in mice [53]. Although the direct mechanisms linking EPHX1 and vitamin D are unclear, an in-vivo study showed that 1,25(OH)₂D₃ increased the expression of EPHX1 and other phase I and phase II biotransformation enzymes in the intestinal tissue of vitamin D-deficient rats [54]. There is no known association between EPHX1 and breast cancer risk [55]. We do note that our results should be interpreted with caution as the EPHX1 CpG site that was strongly associated with 25(OH)D in our sample had a low ICC (< 0.5), meaning that the within-subject variability was larger than the between-subject variability.

SPN encodes a sialoglycoprotein expressed on leukocytes and platelets. Cell culture models have demonstrated that vitamin A and D metabolites upregulate SPN expression [56, 57]. SMARCD2 encodes a critical component of the SWItch/Sucrose Non-Fermentable (SWI/SNF) chromatin remodeling complex, which uses ATP-derived energy to unwrap or restructure chromatin [58]. SMARCD subunits serve as a link between the SWI/SNF core complex and transcription regulators, including nuclear receptors such as VDR and RXR [59–61]. Although we found no prior reports of a direct link between these sites and breast cancer risk, recent studies have demonstrated that genes encoding for the SWI/SNF chromatin-remodeling complex are mutated in approximately 20% of all human tumors [62] and are considered to be critical tumor suppressors [58] and epigenetic regulators of tumorigenesis [63].

The effect measures we estimate for the interaction analysis, RHRs, measure the extent to which the hazard ratio for the 25(OH)D-breast cancer association depends on the epigenetics, as measured by methylation at a particular CpG locus. Because methylation and 25(OH)D were measured in the same blood samples, we cannot address the temporality of the identified associations to determine whether 25(OH)D influences methylation, methylation influences 25(OH)D, or a third factor influences both. Similarly, we can only assess whether the relationship between methylation and 25(OH)D is different for those who later developed breast cancer, suggesting multiplicative interaction, and not whether 25(OH)D or methylation acts as an effect modifier or mediator. Repeated measures of methylation and 25(OH)D would be needed to address temporality. Future studies could also help to determine the most appropriate cut-point for 25(OH)D levels in gene-by-environment interaction studies. We selected 38.0 ng/mL based on our previous results [3] and other findings supporting a threshold effect of similar magnitude [64], but we cannot be sure what levels have the most biological relevance for breast cancer risk.

We limited our sample to non-Hispanic white women to minimize the influence of genetic ancestry. As such, our results may not be fully generalizable. Our sample is also selective in that all participants had a sister diagnosed with breast cancer, and had, on average, approximately twice the risk of developing breast cancer themselves. Our findings are internally valid, but overrepresented risk factors (e.g., germline genetic or early childhood exposures) may inflate the magnitude of effect estimates if they influence the associations evaluated here.

Serum 25(OH)D concentrations were associated with methylation levels at candidate CpGs in vitamin D-related genes and three genes with links to immune function or regulation of VDR. Methylation levels of some of these CpGs may interact with vitamin D to affect breast cancer risk. These results contribute to our understanding of the relationship between vitamin D and DNA methylation and the impact of vitamin D on breast cancer risk.