Prediagnostic 25-Hydroxyvitamin D Concentrations in Relation to Tumor Molecular Alterations and Risk of Breast Cancer Recurrence

Vitamin D is a multifunctional steroid hormone that can be synthesized cutaneously from 7-dehydrocholesterol in response to sun exposure or obtained from natural and fortified foods and supplements (1, 2). Vitamin D undergoes two hydroxylation steps in its activation path to the active hormone, first converted to 25-hydroxyvitamin D [25(OH)D] in the liver and ultimately to calcitrol [1α,25(OH)₂D]—the most biologically active form of vitamin D (1, 2). Although 1α-hydroxylation primarily takes place in the kidney, it also occurs in many, if not most, extrarenal sites such as the breast and breast cancer cells, where vitamin D signal transduction is also modulated in a paracrine or an autocrine manner in nonskeletal sites mediating noncalcium-related actions (1, 2). In all target cells, 1α,25(OH)₂D forms a ligand–bound complex with the vitamin D receptor–retinoid X receptor heterodimer, and together with other comodulators binds to vitamin D response elements to initiate transcriptional cascades that stimulate or reduce selective genes that mediate the actions of the hormone (1–4).

Despite long-standing knowledge of its role in bone mineralization and calcium homoeostasis, compelling experimental studies showed that vitamin D and its analogs exhibit anticarcinogenic properties and can regulate cell proliferation, differentiation, apoptosis, inflammation, invasion, and metastasis in breast tumors (5–7). In vitro and in vivo studies also showed that vitamin D can potentiate the anticarcinogenic effects of chemotherapeutic agents in tumors (8–10). However, in prospective population-based studies and randomized clinical trials, the association between prediagnostic vitamin D levels with breast cancer survival remains inconclusive. Although several meta-analyses showed statistically significant associations between 25(OH)D and improved survival among patients with breast cancer (11–19), a recent large-scale randomized controlled trial (RCT) showed no clear benefit of high-dose vitamin D (2,000 IU/day) on breast cancer–specific mortality over a 5-year follow-up (20). A better understanding of the molecular mechanisms linking prediagnostic vitamin D concentration with breast cancer recurrence would provide important human evidence supporting the role of vitamin D in preventing tumor progression and thereby improving survival.

To assess the potential prognostic effects of vitamin D, we prospectively evaluated the association between prediagnostic circulating 25(OH)D with breast cancer recurrence from invasive breast cancer cases drawn from the Nurses' Health Study (NHS) and NHSII. We next sought to explore molecular underpinnings of 25(OH)D-associated gene expression regulation for breast tumor progression using transcriptome-wide profiling. Since our previous work showed that postdiagnosis vitamin D supplement use was only associated with decreased risk of recurrence among ER-positive, but not ER-negative breast cancer (21), we will examine the associations also by ER status.

This study includes participants from two large-scale prospective longitudinal cohorts of registered female nurses in the United States—the NHS and the NHSII. Established in 1976, the NHS recruited 121,701 women, ages 30 to 55 years, who completed and returned an initial questionnaire. In 1989, the NHSII was initiated, which enrolled 116,429 women, ages 25 to 42 years, who completed and returned an initial questionnaire. Both cohorts followed participants via questionnaires mailed biennially to update information on lifestyles, medications, and certain outcomes. The cumulative follow-up rates were greater than 90% for the NHS and NHSII (22). Invasive breast cancer cases were initially identified by participants' responses to the biennial questionnaires from the start of follow-up to 2012, or through searching the National Death Index for participants who did not respond. With participants' permission, we were able to link 96% of breast cancer cases to relevant medical records. Using established protocols, centralized medical record review confirmed over 99% of reported breast cancer cases. All breast cancer cases included in this analysis had no reported previous history of cancer. In 1993, we started to collect archived formalin-fixed paraffin-embedded (FFPE) breast cancer blocks for participants with primary incident breast cancer and were able to create tissue microarrays (TMA) from 5561 NHS and NHSII participants (Dana-Farber/Harvard Cancer Center Tissue Microarray Core Facility, Boston, MA). Given limited funds, we prioritized women with existing genetic and circulating biomarker measurements, and were able to perform tumor tissue gene expression microarrays for 954 breast cancer cases. The study protocol was approved by the institutional review boards of the Brigham and Women's Hospital and Harvard T.H. Chan School of Public Health, and those of participating registries as required.

Prediagnostic circulating 25(OH)D levels were measured in nested case–control studies of breast cancer in six batches, as published previously (23–25). In brief, 25(OH)D levels were measured using either a high-affinity protein-binding assay (Dr. Michael F. Holick laboratory, Boston University School of Medicine, Boston, MA, or Dr. Bruce W. Hollis laboratory, Medical University of South Carolina, Charleston, SC) and a radioimmunoassay with radioiodinate tracers (Heartland Assays, Inc.). Neither assay could distinguish between 25(OH)D₂ and 25(OH)D₃ and therefore we use the term 25(OH)D without subscript to connote total circulating vitamin D₂ and D₃ metabolites. Coefficients of variances (CV) ranged from 6.0% to 17.6% for blinded quality control samples in each batch (23–25). We controlled for potential batch effects by recalibrating 25(OH)D distribution for each batch to an “average” batch, independent of age, body mass index (BMI), case–control status, and season of blood collection (25, 26).

RNA was extracted from multiple cores of 1 or 1.5 mm taken from FFPE breast cancer tumor tissue (n = 1–3 cores) and normal-adjacent (n = 3–5 cores) tissues using the Qiagen AllPrep RNA Isolation Kit. Normal-adjacent tissues were obtained >1 cm away from the edge of the tumor. A detailed protocol has been published previously (microarray accession number: GSE115577; refs. 27–29). In brief, we profiled transcriptome-wide gene expression using Affymetrix Glue Grant Human Transcriptome Array 3.0 (hGlue 3.0) and Human Transcriptome Array 2.0 (HTA 2.0) microarray chips (Affymetrix). We used robust multiarray averages to perform normalization [RMA; Affymetrix Power Tools (ATP)], log-2 transformed the data, and conducted sample quality control with Affymetrix Power Tools probeset summarization based metrics. A total of 1,577 samples (882 tumor tissues and 695 normal-adjacent tissues) from 954 invasive breast cancer cases passed quality control. Participant characteristics were similar among breast cancer cases with and without gene expression measurement (29). For genes that were mapped by multiple probes, we selected the most variable probe to represent the gene. Our current analyses included 17,791 (70%) genes that were profiled in both platforms. Technical variabilities (batch) were controlled using ComBat—an empirical Bayes method used to control for known batch effects (30). Genes with low expression (<25th percentile) were removed from the analyses.

We obtained information on age, BMI, fasting status, and seasonality via questionnaires collected at the time of blood collection. Race, menopausal hormone therapy (MHT), year of breast cancer diagnosis, characteristics (tumor stage, tumor grade), and treatment information (chemotherapy, radiotherapy, and hormone therapy) were obtained via the biennial NHS and NHSII questionnaires. ER status was determined via central review of breast cancer TMAs, and filled with pathological reports if missing. We computed B cell, CD4⁺ T cell, CD8⁺ T cell, neutrophil, macrophage, and dendritic cell proportions with Tumor IMmune Estimation Resource (TIMER), which uses constrained least squares fitting on informative genes to predict the abundance of immune cell infiltration in bulk tumor tissues (31, 32). Our gene expression analysis included women with complete information on 25(OH)D measurements, transcriptional profiling, and covariates, resulting in 659 breast cancer cases (604 tumor tissues, 473 normal-adjacent tissues; 1,077 total samples). Analysis for recurrence was restricted to subjects with tumor TMAs.

Women who provided a blood sample were quite similar to the overall NHS cohort. For example, we observed (blood vs. overall cohort): mean age (56 years vs. 56 years), BMI (25.5 kg/m² vs. 25.7 kg/m²), parity (3.0 vs. 3.0), oral contraceptive use duration (51 months vs. 50 months), and current smokers (13% vs. 17%). Because NHS cohorts consist of primarily white women, race distributions were also similar comparing the nested case–control participants (96% white) versus the full cohort (97% white). For physical activity (metabolic equivalent tasks), mean is 19.02 (hour/week) for the blood nested case–control and 20.16 (hour/week) for full cohorts.

Breast cancer recurrences are ascertained through supplemental questionnaires sent to living cohort members with a prior confirmed diagnosis of breast cancer. An internal validation study found 92% sensitivity and 92% specificity for the self-report of breast cancer recurrence compared with medical record review. In addition, reported new cancers of the liver, bone, brain, and lung on biennial questionnaires, subsequent to a breast cancer diagnosis, are also assumed to be breast cancer recurrences, as these are common sites of metastasis. Finally, for cohort members who died from breast cancer without having previously reported a recurrence, a recurrence is assumed to have occurred. If the death from breast cancer occurred more than 4 years after the initial diagnosis, the date of recurrence is assigned to 2 years prior to the date of death (the median survival time for metastatic breast cancer). If the death from breast cancer occurred within 4 years of the initial diagnosis, the date of recurrence is assigned to be halfway between the dates of initial diagnosis and death. Approximately 40% of breast cancer recurrences in the cohorts are identified through self-report in questionnaires or new diagnosis at common sites of metastasis, while the remaining 60% are assumed to occur in women who died from breast cancer without reporting a recurrence. Recurrence-free survival is then defined as the time from diagnosis to breast cancer recurrence or end of follow-up (December 2015), censoring individuals who died from causes other than breast cancer at the time of death.

Association of prediagnostic circulating 25(OH)D with breast cancer recurrence was examined by using Cox proportional hazards regression to calculate hazard ratios (HR) and 95% confidence intervals (CI). Recurrence-free survival (RFS) is defined as the time (in months) that elapsed between breast cancer diagnosis and breast cancer recurrence, diagnosis of cancer in common sites of recurrence (i.e., bone, brain, liver, and lung), or death from breast cancer without reported recurrence (29). Prediagnostic circulating 25(OH)D was modeled continuously and as a dichotomized variable at 30 ng/mL, as previous observational studies and Endocrine Society Guidelines suggest that a relatively high plasma level (≥30 ng/mL) is needed to achieve physiologic benefits other than bone health (33–35). We defined plasma 25(OH)D ≥30 ng/mL as “sufficient,” and <30 ng/mL as “insufficient.” In the multivariable Cox regression models, we controlled for age at blood collection (continuous), BMI at blood collection (continuous), year of diagnosis (continuous), race (white, non-white), tumor stage (continuous; I, II, III), tumor grade (continuous; low, intermediate, high), chemotherapy (ever use, never use, missing/unknown), radiotherapy (ever use, never use, missing/unknown), and hormone therapy (ever use, never use, missing/unknown). We assessed statistical interactions by ER status, and additionally stratified the analyses by ER status. The proportional hazards assumption was satisfied by evaluating a time-dependent variable, which was a product term between prediagnostic circulating 25(OH)D and log survival time (P > 0.05). Kaplan–Meier survival curve was constructed stratified by prediagnostic circulating 25(OH)D (≥30 ng/mL or <30 ng/mL). As sensitivity analyses, we restricted study participants to women with early-stage (I and II) breast cancer and excluded women with short lag time between diagnosis and recurrence (within 3 years). We additionally adjusted for season with categorical variables (spring, summer, fall, winter) or with Fourier series terms cos(2π * doy/365.25) and sin(2π * doy/365.25), where doy represents day of year (36). Since our previous study from the same population showed no statistically significant interaction between menopausal status and 25(OH)D levels, we did not stratify the analysis by menopausal status.

For each individual gene, we evaluated the association between prediagnostic circulating 25(OH)D with transcriptome-wide gene expression using covariate-adjusted linear regression (limma; ref. 37). We adjusted in each regression model the following covariates selected a priori: age at blood collection (continuous), year of diagnosis (continuous), menopausal status and MHT use (postmenopausal not using, postmenopausal using, premenopausal, or unknown), BMI at blood collection (continuous), fasting status at blood drawn (fasting, not fasting, not known), season (spring, summer, fall, winter), and surrogate variables generated from the transcriptome data (the leek method from Bioconductor sva package in R; ref. 38). We also explored an alternative method of adjusting for season with Fourier series terms as described previously (36). All analyses were performed separately for tumor and normal-adjacent tissues, and we stratified the analysis by ER status of the tumor. We considered a gene to be significant transcriptome-wide if it met FDR P value of P_BH < 0.05 (39).

Using a competitive gene set testing procedure, correlation adjusted mean rank (CAMERA), we explored functional enrichment of biological pathways associated with prediagnostic circulating 25(OH)D (40). This method estimates the variance inflation factor associated with intergene correlation, and incorporated into rank-based test procedures. CAMERA correctly controls for type I error rate regardless of intergene correlations, and provides valid testing results based on gene permutation even when genes in the test sets are correlated (40). We chose the 50 “hallmark” gene sets from the Molecular Signature Database (MSigDB; http://www.broadinstitute.org/gsea/msigdb/). These “hallmark” gene sets were thought to reduce redundancy while provide effective summary of most of the relevant information of the original founder sets. In these analyses, we controlled for the same set of covariates as for the single-gene analysis, and chose an intergene correlation of 0.01. Figure 1 shows the flow of the analyses. All analyses were conducted in SAS 9.4 and in R 3.1.4.

Among gene expression pathways that were significant after FDR correction in ER-positive tumor samples, we identified the individual genes that were nominally statistically significantly associated with 25(OH)D (P < 0.05) and created the 25(OH)D gene expression score based on unique genes from the significantly enriched pathways. The 25(OH)D gene expression score was calculated as ∑(z-transformed genes from positively regulated pathways) − ∑(z-transformed genes from negatively regulated pathways). We present three sets of Cox models: model 1 adjusted for age at diagnosis (continuous) and year of diagnosis (continuous); model 2 additionally adjusted for tumor stage (continuous; I, II, III) and tumor grade (continuous; low, intermediate, high); and model 3 additionally adjusted for treatment modality, that is, chemotherapy (ever use, never use, missing/unknown), radiotherapy (ever use, never use, missing/unknown), and hormone therapy (ever use, never use, missing/unknown). We tested the violation of the proportional hazards assumption by evaluating an interaction term between 25(OH)D score and log survival time (P > 0.05). We also conducted a sensitivity analysis by further restricting to women with early-stage (I and II) breast cancer.

To assess if the 25(OH)D gene expression score derived in the NHS was generalizable to other populations, we leveraged an independent dataset, Molecular Taxonomy of Breast Cancer International Consortium (METABRIC). Using the same gene signature derived in NHS (described above) for ER-positive breast cancer, we calculated the 25(OH)D score using the available gene expression data in ER-positive tumor in METABRIC. It is important to note, METBRIC does not have circulating vitamin D and the calculation of the score is based on that from NHS. We then examined if this score was associated with recurrence-free survival. We included 737 stage I to III ER-positive breast cancer cases with complete information on transcriptional profiling, breast cancer–specific survival data, and covariate information. Transcriptional profiling was performed on the Illumina HT-12 v3 platform (41, 42). All clinical and gene expression data are deposited at the European Genome-Phenome Archive (EGA, http://www.ebi.ac.uk/ega/) under accession number EGAS00000000083 (41, 42). We present three sets of Cox regression models, adjusting for the same covariates as described in the previous section. We restricted 25(OH)D score to the same range as for the NHS studies, and further restricted the analyses to women with early-stage (I and II) breast cancer.

Among 604 women who contributed to tumor data, 215 participants had “sufficient” levels of prediagnostic circulating concentration of 25(OH)D (≥30 ng/mL) whereas 389 had “insufficient” levels (<30 ng/mL). Mean age at blood draw was slightly younger in women with sufficient levels of 25(OH)D compared with insufficient levels [mean (SD) = 53.4 (8.6) vs. 54.8 (8.2) years; Table 1]. In the sufficient 25(OH)D group, 17% of participants were postmenopausal and had never used MHT, 33% of women were postmenopausal and used MHT, and 50% were either premenopausal or did not report MHT use. A higher proportion of women in the insufficient 25(OH)D group were postmenopausal and had not used MHT. Average BMI was 24.5 kg/m² (SD = 4.1) in the sufficient 25(OH)D group and was 25.9 kg/m² (SD = 4.9) in the insufficient group. Participants in this study were predominantly white (96%). Dietary vitamin D was higher in the sufficient 25(OH)D group compared with the insufficient group [mean (SD) = 518.6 IU/day (285.5) vs. 445.6 (267.2)]. We observed Gaussian distribution of 25(OH)D levels for the study population (Supplementary Fig. S1). Median time between blood draw and cancer diagnosis was 6.4 years (IQR = 7.2 years), whereas median time between blood draw and tumor microarray assessment was 23 years (IQR = 6.1 years).

Prediagnostic circulating 25(OH)D was not associated with overall (ER-positive and ER-negative) breast cancer recurrence (median follow-up = 13.0 years); however, a statistically significant interaction (P = 0.005) was observed between 25(OH)D and ER status (Table 2). We next stratified the analysis by ER status, and observed that higher prediagnostic circulating 25(OH)D was associated with lower risk of recurrence among ER-positive cases, but not ER-negative breast cancer cases. Every 5 ng/mL increase in 25(OH)D level was associated with a 13% decrease in risk of recurrence in women with ER-positive cancers (HR = 0.87; 95% CI, 0.76–0.99). We also observed a similar reduction in risk in ER-positive diseases when we examined 25(OH)D as dichotomized variable, although it was not significant (≥30 ng/mL vs. <30 ng/mL; HR = 0.72; 95% CI = 0.42–1.22; Table 2; Supplementary Fig. S2). An analysis restricting women to early-stage (I and II) breast cancer yielded similar results (Table 2). After excluding women with short lag time between diagnosis and recurrence (≤3 years), we also observed comparable association between 25(OH)D (dichotomized) and recurrence (HR = 0.72; 95% CI, 0.52–0.99; P-value = 0.05). Adjusting for season using categorical variables or using Fourier series terms also gave comparable results; every 5 ng/mL increase in 25(OH)D level was associated with a 14% decrease in risk of recurrence for ER-positive disease (HR = 0.86; 95% CI, 0.76–0.99). Among ER-negative cases, null associations were observed when we modeled 25(OH)D continuously (HR = 1.07; 95% CI, 0.91–1.27; Table 2) or as a categorical variable (HR = 1.34; 95% CI, 0.54–3.29; Table 2).

Overall, we did not observe a statistically significant association between prediagnostic circulating 25(OH)D and any single gene in tumor or normal-adjacent tissues stratified by ER status (FDR < 5%). We present in Supplementary Fig. S3 differentially expressed genes with relatively large effect size (i.e., genes with log fold change larger than the average of the 2.5th and 97.5th percentile) and had a nominal P value less than 0.001 in tumor and normal-adjacent ER-positive samples.

Competitive gene set enrichment analysis identified 21 pathways significantly enriched in ER-positive tumor tissues, and three significantly enriched in ER-positive normal-adjacent tissues (FDR < 5%; Table 3). Those pathways included genes implicated in proliferation and mitogenesis, cellular stress, immune response, and cancer cell metabolism. Although all those pathways were downregulated in ER-positive tumor tissues except for KRAS signaling, gene sets involved in cellular stress and proliferation were upregulated in normal-adjacent ER-positive tissues. In ER-negative tumor and normal-adjacent tissues, we found nine pathways that were differentially enriched (FDR < 5%) in relation to high 25(OH)D (Table 4). The biological processes included proliferation, cellular stress, tumor microenvironment, and immune response. However, we observed heterogeneous molecular response between ER-positive and ER-negative tumors—overlapping pathways (n = 5) identified in ER-positive and ER-negative tumors tended to have opposite directions of association.

We considered alternative cut-points for defining ER positivity. When we defined ER positivity as ≥10% of cells staining positive, the association tended to become stronger between 25(OH)D and molecular pathways (measured by the magnitude of FDR-corrected P values), suggesting a potential cross-talk between 25(OH)D and ER (Supplementary Table S1). In addition, we observed more CD8⁺ T-cell infiltration in ER-negative tumor as compared with ER-positive tumor (Supplementary Fig. S2). As a sensitivity analysis, we introduced an interaction term [i.e., time between blood draw and cancer diagnosis and 25(OH)D], and found that pathway analysis results were unlikely to be affected by duration between blood draw and diagnosis (P_interaction ranges from 0.15 to 0.98). Adjusting for season using Fourier series term and as categorical variables yielded similar results (Supplementary Table S2).

Among women with stage I to III ER-positive breast cancer, a 10-unit increase in 25(OH)D gene expression score was marginally associated with a 18% reduction in risk of recurrence (HR = 0.77; 95% CI, 0.59–1.01; P-value = 0.058) in the NHS cohorts (Table 5, model 3; Fig. 2). We observed similar magnitude of effect estimates when we restricted the analysis to participants with early-stage (I and II) breast cancer. However, when we examined the generalizability of these findings in another dataset, we did not find that this score was associated with breast cancer recurrence among ER-positive cancer cases in METABRIC (Fig. 2; Supplementary Tables S3 and S4). In the fully adjusted model in METABRIC, HR was 1.02 for every 10-unit increase in 25(OH)D gene expression score (95% CI, 0.92–1.15; P-value = 0.617) for stage I to III breast cancer cases (Supplementary Table S4). Findings remained similar when we restricted 25(OH)D gene expression scores to the same range as for NHS (Fig. 2).

In this prospective population-based study, higher prediagnostic circulating 25(OH)D (measured 6.4 years before diagnosis at one time point) was inversely associated with breast cancer recurrence among women with ER-positive in the NHS and NHSII, but not ER-negative breast cancer. After multiple testing corrections, no single gene was significantly associated with 25(OH)D in either breast tumor or adjacent tissues for ER-positive and ER-negative diseases. However, pathway analysis identified multiple gene sets that were significantly enriched and downregulated at a transcriptome-wide threshold in ER-positive breast tumor, which included genes associated with proliferation, apoptosis, migration and invasion, inflammation, and cancer cell metabolism. Further, among ER-positive breast cancers, we observed a marginally significant inverse association between the composite score summarizing the gene expression of 25(OH)D signatures in tumor tissue and breast cancer recurrence in the NHS cohorts. However, the 25(OH)D score derived on the basis of NHS data did not translate to a publicly available dataset of METABRIC, suggesting that further refinement is needed to improve generalizability.

An increasing body of evidence supports an inverse association between vitamin D and cancer mortality (11–19). Large-scale meta-analyses of 17,332 patients with cancer showed that higher circulating 25(OH)D—measured at or near the time of cancer diagnosis—was associated with improved survival among patients with lymphoma, colorectal, and breast cancer (15). Among 4,413 breast cancer cases in a separate meta-analysis, the pooled HRs comparing highest with lowest quintiles were 0.65 (95% CI, 0.49–0.86) for total mortality and 0.58 (95% CI, 0.38–0.84) for breast cancer–specific mortality (14). A recent case–control study also revealed that serum circulating 25(OH)D, measured at the time of diagnosis, was associated with lower risk of breast cancer–specific mortality (17). Compared with the lowest tertile, women in the highest tertile of 25(OH)D levels had a 28% decrease in total mortality (HR = 0.72; 95% CI, 0.54–0.98; ref. 17). A meta-analysis of current RCTs using high doses of vitamin D supplementation showed significant improvement in breast cancer survival with significantly reduced total cancer mortality, but did not demonstrate a reduction in total cancer incidence (25). The large vitamin D and omega-3 trial (VITAL)—which assessed the effect of high-dose vitamin D supplementation (2,000 IU/day) and omega 3 (1 g/day) on cancer incidence and mortality—showed no association between prediagnostic vitamin D supplementation and total death or breast cancer–specific death (20). However, a protective effect of vitamin D was observed for overall cancer mortality after excluding cases that occurred in the first 2 years of follow-up (20). Because our median follow-up time was 12.8 years, a longer follow-up for VITAL may be needed to observe potential beneficial effects of vitamin D on breast cancer survival.

Consistent with previous findings in the NHS cohorts showing that postdiagnostic vitamin D supplement use was associated with improved breast cancer survival in ER-positive diseases (21), we observed that higher prediagnostic circulating 25(OH)D was inversely associated with breast cancer recurrence only among ER-positive breast cancer cases. A heterogeneous molecular response was also observed between ER-positive and ER-negative breast tumors. In addition to the relatively small sample size of the ER-negative subgroup, we postulate several biological reasons for the observed discrepancy with the ER-positive subgroup. Previous experimental studies have shown that 1α,25(OH)₂D suppress aromatase and COX-2 expression in both cancerous epithelial cells and stromal adipose fibroblasts in the breast (43). Aromatase and COX-2 are key enzymes in the estrogen signaling pathway, and help convert androgenic precursors and arachidonic acid to potent proliferative compounds of estrogen and prostaglandin E2 (PGE2). These findings support a potential cross-talk between 1α,25(OH)₂D and inhibition of estrogen signaling in promoting cellular proliferation. Supported by our data, the strength of the association between 25(OH)D and proliferation-related gene expression pathways became more profound (measured by FDR) when we changed the definition of ER-positivity from ≥1% of tumor nuclei stained positive to ≥10% of tumor nuclei stained positive. It is well documented that ER-negative tumor is more heterogeneous than ER-positive tumor (44). In our study, RNA was extracted from 1 or 1.5 mm cores from bulk tumor tissues. Therefore, it is difficult to rule out the possibility that a more heterogeneous tumor subtype gives rise to less robust biological signals. Using TIMER, a mathematical algorithm to deconvolute tumor cellular heterogeneity, we found more infiltrated CD8⁺ T cells in ER-negative tumors as compared with ER-positive tumors.

Given that higher prediagnostic circulating 25(OH)D was associated with better survival among ER-positive breast cancers in our data, we derived a gene expression score for 25(OH)D using differential enriched pathways and found that downregulation of the 25(OH)D gene expression signature was associated with lower risk of breast cancer recurrence in the NHS cohorts. Several in vitro experiments also studied the impact of 1α,25(OH)₂D on differential gene expression at a transcriptome-wide level (45, 46). Consistent with our findings, in cell lines derived from normal mammary tissue and breast cancer cells, gene set enrichment analysis highlighted several pathways (including cell cycle and proliferation, immunomodulation, cell metabolism, apoptosis, and cell death, as well as cellular adhesion, migration, and invasion) being altered by 1α,25(OH)₂D (45, 46). Nevertheless, our 25(OH)D score may not be generalizable to other datasets, as we observed null association between 25(OH)D score with breast cancer recurrence in METABRIC. This lack of consistent findings between datasets may reflect both real biologic differences in tumors (e.g., stage, grade) and study populations and/or methodologic differences. The METABRIC dataset consists of samples that were assayed in multiple labs and batches, whereas the NHS/NHSII were conducted in two batches using similar array technology. Additional research in other study populations is warranted to determine whether our findings can be confirmed in other studies.

We hypothesize that prediagnostic factors not only influence risk, but also recurrence. This is especially evident for vitamin D, because experimental studies suggest that vitamin D inhibits tumor growth and proliferation (43, 47–50), but data are limited on the effect of vitamin D on tumor initiation. A previous randomized trial of high-dose vitamin D (VITAL) also showed that vitamin D randomization was associated with reduced overall cancer mortality after excluding cases that occurred in the first 2 years of follow-up (but not risk; ref. 20). Prediagnostic factors may be a reflection of a person's general health status and well-being (immunity, damage, repair mechanisms, etc.), which are all important determinants for one's ability to recover or survive postdiagnosis.

Our study has several strengths. To the best of our knowledge, it is the first study linking prediagnostic circulating 25(OH)D, differential gene expression in breast tumor, and breast cancer recurrence in human population–based studies. Our study included a relatively large sample of ER-positive cases, which enabled us to robustly identify differential survival data and gene expression patterns in relation to circulating 25(OH)D levels. Vitamin D supplementation represents a personal level secondary cancer prevention nutritional factor, which could be easily implemented and integrated into daily life to reduce breast cancer disease burden and added to treatment protocols to improve treatment results.

There are also limitations to our study. Our FFPE tissue blocks were archived for 9 to 30 years, and differential RNA degradation may be of concern. To minimize variations arising from different sample storage conditions and the age of FFPE tissue, we adjusted for year of diagnosis in all regression models. Our pilot work also showed high correlations between ESR1, PGR, and ERBB2 expression with ER, PR, and HER2 IHC staining (27, 28). Our microarray analysis was based on RNA extracted from bulk tumor and adjacent tissues. Since ER-negative tumor is often more heterogeneous than ER-positive tumor, 25(OH)D-associated biological signals may be masked in a more heterogeneous tissue subtype. Future research with larger sample size is warranted to study the effect of vitamin D on disease progression in ER-negative breast tumors. Finally, our study had only a single measurement of each subject's prediagnostic levels of 25(OH)D. It would have been better to have follow-up assays but previous work in this population of subjects noted the feasibility of using a single measurement of multiple biomarkers including 25(OH)D that can reliably estimate average levels over a 1- to 3-year period in epidemiologic studies (51). Further, in a previous study that had prediagnostic plasma 25(OH)D measured approximately 10 years apart (1989–1990 to 2000–2002) from the same subset of controls (n = 238) with blood collected in concordant seasons, intraclass correlation (ICC) were 0.51, 95% CI (0.42–0.60; ref. 52). Potential survivorship bias may occur for analyses with recurrence as an outcome since not all breast cancer cases reported their recurrence status. For those cancer cases who did not return supplementary questionnaires, did not die, and did not report subsequent new cancers at the liver, bone, brain, and lung, we assumed that they did not have recurrence. Nonetheless, such misclassification is likely to be nondifferential because missingness to report a recurrence would not be related to the exposure variable [prediagnostic 25(OH)D] and would bias the results towards the null. Our study did not have information on postdiagnostic circulating 25(OH)D. Vitamin D supplementations may be prescribed to women after cancer treatment due to concerns about bone density; therefore, among breast cancer cases, prediagnostic 25(OH)D level may not truly reflect long-term circulating 25(OH)D concentrations, especially during and immediately after cancer treatment. However, given its temporality, postdiagnostic 25(OH)D would not confound our associations. Future studies are warranted to include postdiagnostic 25(OH)D measurements to further verify our conclusions.

Our findings suggest that higher prediagnostic circulating 25(OH)D was associated with decreased risk of recurrence among patients with ER-positive breast cancer in the NHS, potentially through downregulating signaling pathways in tumor that are primarily involved in cell proliferation, apoptosis, migration and invasion, inflammation, and cancer cell metabolism. Our results lend further support to the use of vitamin D as a potential secondary prevention strategy for reducing cancer progression in ER-positive breast cancer. However, the 25(OH)D score generated on the basis of NHS did not translate to publicly available data, suggesting that further refinement is need to improve the generalizability of the score. Future studies which include postdiagnostic 25(OH)D measures would also be important to verify our findings.

No potential conflicts of interest were disclosed.

Conception and design: C. Peng, R.M. Tamimi

Development of methodology: C. Peng, R.M. Tamimi

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): Y.J. Heng, S.E. Hankinson, A.H. Eliassen, R.M. Tamimi

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): C. Peng, N.C. DuPre, K.H. Kensler, K. Glass, O.A. Zeleznik, P. Kraft, D. Feldman, A.H. Eliassen, R.M. Tamimi

Writing, review, and/or revision of the manuscript: C. Peng, Y.J. Heng, D. Lu, N.C. DuPre, K.H. Kensler, K. Glass, O.A. Zeleznik, P. Kraft, D. Feldman, S.E. Hankinson, K. Rexrode, A.H. Eliassen, R.M. Tamimi

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): Y.J. Heng, A.H. Eliassen, R.M. Tamimi

Study supervision: A.H. Eliassen, R.M. Tamimi

This work was supported by the NIH/NCI [UM1 CA186107 (to A.H. Eliassen), UM1 CA176726 (to A.H. Eliassen), P01 CA87969 (to A.H. Eliassen and R.M. Tamimi), U19 CA148065 (to R.M. Tamimi and P. Kraft), R01 CA178263 (to R.M. Tamimi and K.M. Rexrode), and R01 CA166666 (to S.E. Hankinson)], National Institute of Health Epidemiology Education Training Great to NCD [NIH T32 CA09001 (to A.H. Eliassen)], the Susan G. Komen for the Cure [IIR13264020, SAC110014 (to R.M. Tamimi and S.E. Hankinson)], and the Prevent Cancer Foundation (to C.Peng). The authors thank all participants and coordinators of the Nurses' Health Study and the Nurses' Health Study II for their valuable contribution, and the cancer registries in the following states for their help: AL, AZ, AR, CA, CO, CT, DE, FL, GA, ID, IL, IN, IA, KY, LA, ME, MD, MA, MI, NE, NH, NJ, NY, NC, ND, OH, OK, OR, PA, RI, SC, TN, TX, VA, WA, and WY. The authors assume full responsibility for analyses and interpretation of these data.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.