Plasma 25 hydroxyvitamin D level and blood gene expression profiles: a cross-sectional study of the Norwegian Women and Cancer Post-genome Cohort.

BACKGROUND/
OBJECTIVES: Vitamin D deficiency has been associated with increased risk of developing several diseases, but much is unknown about the molecular effects involved. Gene expression technology is increasingly being used to elucidate molecular mechanisms related to nutritional factors, and in this study of free-living, middle-aged Norwegian women, we aimed at identifying gene expression pathways in the blood associated with vitamin D status. SUBJECTS/
METHODS: Blood samples and questionnaires were collected as a part of the Norwegian Women and Cancer Post-genome Cohort (500 invited subjects, 218 included). Plasma 25 hydroxyvitamin D (25(OH)D) concentrations were measured using high-performance liquid chromatography, and we compared groups with sufficient versus deficient vitamin D status (25(OH)D >50 nmol/l (n=66) versus <37.5 nmol/l (n=83)), to identify differences in gene expression profiles obtained using full-genome microarrays.
RESULTS: In a targeted pathway-level analysis, several immunological processes, immune cell functions and major signaling pathways were differentially regulated according to vitamin D status (P<0.01). To a certain degree, results from in vitro studies reported in the literature were reflected in this population setting.
CONCLUSIONS: We conclude that vitamin D status measured as 25(OH)D was associated with molecular pathways that may ultimately affect the potential onset of diseases. The use of gene expression analysis in a population setting may give valuable input to the study of effects of nutritional factors.

Introduction

Vitamin D is synthesized in the skin following sun exposure, or provided through the diet mainly from fatty fish, fortified foods or supplements. The first hydroxylation reaction required for biological activation takes place in the liver to form 25 hydroxyvitamin D (25(OH)D), and a second reaction occurs in the kidneys to yield the active 1,25 dihydroxyvitamin D (1,25(OH)2D). 1,25(OH)2D binds the vitamin D receptor (VDR) and exerts its biological function through induction of target genes. The first hydroxyl intermediate (25(OH)D) is routinely used as a biomarker for vitamin D exposure.[1, 2] Classical mechanisms influenced by vitamin D include calcium and phosphate metabolism. However, in recent years, an increasing number of collectively termed non-classical mechanisms have been discovered. These include regulation of immunological functions, cell growth, cell differentiation and apoptosis.[1, 3] Non-classical mechanisms may affect the pathogenesis of a wide range of diseases: high levels of 25(OH)D have been associated with reduced risk for diabetes type 1, as well as beneficial effects on autoimmunity, carcinogenesis and vascular health.[4, 5, 6] However, causal relationships have not been established. Owing to uncertainties regarding optimal vitamin D level and the potential for disease prevention, the questions of dosage and supplement use by the general population are heavily debated.[1, 7]

It has been estimated that vitamin D, via the VDR, may alter the expression levels of up to 5% of genes in the human genome.[4] The VDR not only activates gene transcription in a ligand-dependent manner with its dimerization partner retinoid X receptor, but may also repress gene transcription.[8] Multiple genetic and epigenetic elements act to fine-tune the transcriptional activity of VDR,[9] and the specificity and dynamics of these mechanisms are under continuous investigation.

Technologies that allow high-throughput, genome-wide approaches are increasingly being used to study the molecular impact of nutritional components. Using these methods, it has become evident that VDR action is tightly regulated according to cell type and physiological conditions,[10] but the multiple variables influencing potential health effects are not fully understood. Moreover, in vitro and animal studies, as well as investigations of specific sub-populations or groups using defined supplement regimens, may yield results regarding vitamin D mechanisms that are not directly applicable to the general population. More information is needed, especially regarding non-classical effects and potential for disease prevention in the general population. To explore molecular-level mechanisms at the population level, we investigated associations between plasma 25(OH)D and blood gene expression profiles, in a cross-section of middle-aged Norwegian women. We have previously shown that lifestyle factors such as body mass index (BMI) and smoking,[11] hormone therapy use[12] and environmental pollutants[13] are mirrored in the blood gene expression profiles of women in this cohort.

Materials and methods

Study participants, blood sample collection and inclusion criteria

The Norwegian Women and Cancer Study (NOWAC[14]) consists of a representative study population of 1 72 000 women. From the original cohort, more than 50 000 women born in 1943–57 were randomly recruited to the NOWAC Post-genome Cohort,[15] from which 500 women were randomly selected for the current study.[15] The study was approved by the Norwegian Data Inspectorate and the Regional Committee for Medical Research Ethics, and all participants provided written informed consent. Information about anthropometric and lifestyle factors was extracted from a two-page questionnaire answered at the time of phlebotomy (spring 2005). Sex hormones and blood gene expression were investigated in the same samples;[12] therefore, the study group included post-menopausal women only. Inclusion criteria comprised successful donation of two blood samples: one PAXgene Blood RNA tube (Preanalytix, Qiagen, Hilden, Germany), which stabilizes the gene expression profile of all circulating cells,[16] and one tube of citrate-buffered blood plasma. The blood samples were required to be received at the study center and frozen within 3 days after blood draw. Total RNA was extracted using the PAXgene Blood RNA Isolation Kit (Preanalytix), and sufficient RNA quantity, integrity or purity was required for inclusion. Further, at least 40% of the microarray probes had to have signal to noise ratio (S/N) ⩾3, and the analyses of plasma 25(OH)D had to be successful. Inclusion criteria were fulfilled by 249 women. Marine fatty acids (eicosapentaenoic acid and docosahexaenoic acid) were measured in plasma and treated as a potential confounders. Those who had smoked during the week before blood draw were defined as smokers. Thirty-one women were removed from the data set due to fasting before blood draw, or missing information on either fasting status or BMI. The number of included participants was 218.

Microarray analysis and preprocessing of data

Methods for RNA extraction, microarray experiments and preprocessing of data were described in detail elsewhere.[11] Full-genome mRNA expression levels were analyzed using the Applied Biosystems expression array system (Foster City, LA, USA). Briefly, total RNA was amplified, labeled and hybridized to AB Human Genome Survey Microarray V2.0. The AB Expression System software was used to export signal intensities, signal to noise ratios and flagging values. Gene-wise intensities were adjusted for technical variability.[11]

Plasma 25(OH)D measurements

Analysis for 25(OH)D was performed according to a modified version of the method described by Aksnes.[17] Briefly, 0.25 ml plasma samples were spiked with 3H-25(OH)D3 for calculation of recovery, and 25(OH)D was extracted with methanol and n-hexane. The n-hexane phase was collected, evaporated to dryness and ejected into a reverse-phase high-performance liquid chromatography system. 25(OH)D was eluded with methanol/water (85:15, v/v) and the elute was monitored at 265 nm by a diode-array detector (UV6000; ThermoFinnigan, San Jose, CA, USA) with a 5-cm detector cuvette. External quality controls from the vitamin D External Quality Assessment Scheme (DEQAS, London, UK) and Chromsystems (Munich, Germany) were used. Mean recovery of 25(OH)D was 77.2% (s.d. 3.9%) and the inter-assay variation was 6%, with a detection limit of 6.0 nmol/l. Study participants were divided into three categories according to 25(OH)D concentration: deficient vitamin D status (<37.5 nmol/l, n=83), sub-optimal (38–50 nmol/l, n=69) and sufficient (>50 nmol/l, n=66). The cut-offs reflect vitamin D biology and public health effects, and are in line with the dietary reference intakes.[1]

Statistical analysis

Gene expression profiles of the sufficient and deficient groups were compared using R version 2.13.1 (www.cran.r-project.org) with the Bioconductor package. The analyses were adjusted for technical variability according to Dumeaux et al.[11] (time for phlebotomy to freezing of the sample, RNA extraction date and microarray lot number). Gene-wise linear models (limma) were used to evaluate differences in single-gene expression levels, and the Global Test R package[18] was used to test for overall significance. As recommended in the Global Test documentation (www.bioconductor.org) when overall significance is low, we proceeded with a targeted, pathway-level analysis of 78 selected gene sets, also using the Global Test. Gene sets were either curated from the literature or related to general hypotheses on vitamin D function extracted from review articles (overview in Table 1, described in Supplementary Methods, gene set details in Supplementary Table S1, gene symbols in Supplementary Table S2). Gene sets from experimental publications (12 publications provided 1–9 gene sets each) were derived from cell lines, primary cultures or xenograft mice treated with vitamin D metabolites/analogs (Table 1). The general hypotheses were used to generate a list of 41 gene sets (Table 1), and corresponding gene symbols were extracted from Molecular Signatures Database (Broad Institite).[19] The Global Test was adjusted for multiple testing using Benjamini–Hochberg's false discovery rates.[20] Comparative P-values were calculated for each gene set and indicate the proportion of random gene sets of the same size as the tested gene set being significant by chance.[21] Gene dendrograms from the Global Test were used to evaluate which single genes contribute the most to a significant gene set result.[18] Importantly, a single gene that significantly influences the gene set result is not necessarily differentially expressed above the significance threshold, when examined on its own using linear models. However, the gene dendrograms may be used to aid functional interpretation (see www.bioconductor.org for statistical details). Differences in age, BMI, marine fatty acids, use of medication and smoking between the sufficient and deficient groups were evaluated using linear models and χ2 tests. The Global Test was adjusted for the variables that were significantly different between the two groups (BMI and marine fatty acids). The data set has been submitted to Gene Expression Omnibus, accession number GSE15289.

Table 1

Overview of the gene sets used as input for the Global Test

Source of gene set	n publications	n gene lists	References
Primary cultures	4	19	[24, 29, 31, 32]
Cell lines	5	10	[8, 25, 33, 34, 35]
Combination	1	2	[36]
Xenograft mice	1	4	[28]
Review publications	2	2	[37, 38]
General hypotheses	NAa	41	NAa

Abbreviation: NA, not applicable.

Several reviews, scientific reports and text books served as sources for general hypotheses about vitamin D function. Relevant references are included in the Introduction.

Results

Characteristics of the study population are presented in Table 2. Mean concentration of 25(OH)D was 43.4 nmol/l (median 43.2 nmol/l), and the concentrations were normally distributed (Supplementary Figure S1, Supplementary Table S3). Compared with women in the deficient vitamin D group (25(OH)D<37.5 nmol/l), women in the sufficient group (25(OH)D>50 nmol/l) had significantly lower BMI and higher concentrations of marine fatty acids (Table 3).

Table 2

Characteristics of study population (n=218)

	Mean (min−max), or frequencies
Age (years)	55.6 (48–62)
25(OH)D (nmol/l)	43.4 (14.4–91.2)
BMI (kg/m2)	25.4 (16.7–40.5)
Marine fatty acids (mg/ml)a	0.23 (0.07–0.73)
Smoking	55 (25%)
Use of medicationb	122 (56%)

Abbreviations: BMI, body mass index; max, maximum; min, minimum; 25(OH)D, 25 hydroxyvitamin D.

Marine fatty acids: sum of eicosapentaenoic acid and docosahexaenoic acid.

Missing information on medication: n=2.

Table 3

Characteristics of comparison groups

	Deficient	Sub-optimal	Sufficient	P-valuea
25(OH)D cut-off (nmol/l)	<37.5	37.5–50	>50
Group, n	83	69	66
Age (years)	55.9 (48–62)	55.4 (49–62)	55.5 (49–62)	0.79
25(OH)D (nmol/l)	29.3 (14.4–37.4)	44.0 (37.6–49.7)	60.6 (50.3–91.2)	<0.01
BMI (kg/m2)	26.8 (16.7–40.5)	25.0 (18.5–35.1)	24.0 (18.7–31.4)	<0.01
Marine fatty acids (mg/ml)b	0.19 (0.07–0.4)	0.25 (0.09–0.73)	0.26 (0.08–0.55)	<0.01
Smoking	22 (27%)	16 (23%)	17 (26%)	0.89
Use of medication	48 (59%)c	40 (58%)	34 (52%)	0.61

Abbreviations: BMI, body mass index; max, maximum; min, minimum; 25(OH)D, 25 hydroxyvitamin D.

Values are given as mean (min−max), or frequencies.

P-values for the difference between sufficient and deficient group.

Marine fatty acids: sum of eicosapentaenoic acid and docosahexaenoic acid.

Missing information on medication in the deficient group: n=2.

When comparing sufficient versus deficient vitamin D status using gene-wise linear models, no single genes were differentially expressed. However, in a targeted pathway-level approach, the Global Test identified 26 significantly differentially expressed gene sets (Table 4). Ten of the significant gene sets were derived from experimental data from five different publications reporting treatment of cell lines, primary cultures, or xenograft mice with either the biologically active 1,25(OH)2D or vitamin D analogs. Out of the significant gene sets corresponding to general hypotheses about vitamin D function (Table 4), a majority were related to the immune system and immune cell function such as cytokine production (gene set 52), innate immunity signaling (gene set 30), T-cell receptor signaling (gene set 63) and Toll-like receptor (TLR) signaling (gene set 64 and 65). Several gene sets described ubiquitous signaling pathways including the p38/mitogen-activated protein kinase (MAPK) pathway (gene set 70), death signaling (gene set 66) and vascular endothelial growth factor (VEGF) signaling (gene set 44).

Table 4

Differentially expressed gene sets when comparing expression profiles from groups with sufficient versus deficient vitamin D statusa

Set no.	Set size	Genes tested	Gene set description	Treatment	Experimental model	Reference	Adj. Pb	FDR (%)	Comp. P
11	12	4	Validated by RT-PCR	1,25(OH)2D, 6, 24, 48 h	RWPE1 prostate epithelial cells	[25]	<0.01	3	0.01
19	50	38	Top 50 upregulated	1,25(OH)2D3, 4 h	THP-1 monocytic leukemia cells	[8]	<0.01	6	0.00
10	215	135	Differentially expressed at all time points	1,25(OH)2D, 6, 24, 48 h	RWPE1 prostate epithelial cells	[25]	0.01	11	0.00
21	100	75	Top 50 up- and top 50 downregulated	1,25(OH)2D3, 4 h	THP-1 monocytic leukemia cells	[8]	0.01	11	0.00
55	33	27	BIOCARTA_IL1R_PATHWAY				0.01	11	0.02
64	102	75	KEGG_TOLL_LIKE_RECEPTOR (hsa04620)	1,25(OH)2D3, 24 h	Primary osteoblast cultures	[24]	0.01	11	0.00
65	37	35	BIOCARTA_TOLL_PATHWAY	1,25(OH)2D3, 24 h	Primary osteoblast cultures	[24]	0.01	11	0.01
70	40	39	BIOCARTA_P38MAPK_PATHWAY	1,25(OH)2D3, 24 h	Primary osteoblast cultures	24	0.01	11	0.01
78	50	20	REACTOME_NUCLEAR_RECEPTOR_TRANSCRIPTION				0.01	11	0.01
66	33	33	BIOCARTA_DEAHT_PATHWAY	1,25(OH)2D3, 24 h	Primary osteoblast cultures	[24]	0.02	12	0.02
71	58	46	BIOCARTA_PPARA_PATHWAY				0.02	12	0.01
52	72	48	GO_CYTOKINE_PRODUCTION (GO:0001816)				0.02	12	0.02
1	15	8	Top 15 upregulated	1,25(OH)2D3, 24 h	Primary osteoblast cultures	[24]	0.02	13	0.06
22	43	34	Validated by RT-PCR	1,25(OH)2D3, 4 h	THP-1 monocytic leukemia cells	[8]	0.03	14	0.05
42	12	5	SMOOTH_MUSCLE_CONTRACTION (GO:0006939)				0.03	14	0.08
49	190	138	KEGG_CHEMOKINE_SIGNALING (hsa04062)	TX527	CD3+ T cells	[29]	0.05	14	0.00
63	108	83	KEGG_T_CELL_RECEPTOR_SIGNALING (hsa04660)				0.03	14	0.02
3	30	15	Top 15 up- and top 15 downregulated	1,25(OH)2D3, 24 h	Primary osteoblast cultures	[24]	0.04	14	0.10
15	93	40	Top 100 downregulated	MT19c	SKOV-3 ovarian xenograft mice	[28]	0.04	14	0.06
27	366	292	REACTOME_SIGNALING_IN_IMMUNE_SYSTEM				0.04	14	0.00
30	136	95	REACTOME_INNATE_IMMUNITY_SIGNALING				0.03	14	0.03
31	137	102	KEGG_INSULIN_SIGNALING (hsa04910)				0.04	14	0.02
44	76	58	KEGG_VEGF_SIGNALING (hsa04370)				0.04	14	0.05
50	20	16	BIOCARTA_GCR_PATHWAY				0.04	14	0.09
13	285	191	All downregulated	MT19c	SKOV-3 ovarian xenograft mice	[28]	0.05	14	0.02
34	62	28	REACTOME_STEROID_METABOLISM				0.05	14	0.10

Abbreviations: adj. P, adjusted P-value; comp. P, comparative P-value; FDR, false discovery rate; 1,25(OH)2D3, 1,25 dihydroxyvitamin D3; RT–PCR, PCR with reverse transcription.

The table is sorted by FDR. Gene sets without additional information (treatment, experimental model and reference) correspond to general hypotheses about vitamin D function.

Vitamin D status measured as plasma 25(OH)D concentration, sufficient group: 25(OH)D>50 nmol/l (n=66), deficient group: 25(OH)D<37.5 nmol/l (n=83).

P-values adjusted for body mass index and marine fatty acids (sum of eicosapentaenoic acid and docosahexaenoic acid).

For the two most significant gene sets (gene set 11 and 19, false discovery rate<10%), clustering graphs and dendrograms showing the single genes that contribute significantly to the result of the respective gene sets are given in Figure 1. Among the 26 significant gene sets, a total of seven single genes were identified as significantly contributing to the gene set results: cluster of differentiation 14 (CD14, significant in eight gene sets from five different sources), thrombomodulin (THBD, significant in three gene sets from two different sources) and the following genes found in one gene set each: mitogen-activated protein 14 (MAPK14), prokineticin 2 (PROK2), Toll-like receptor 2 (TLR2), tumor necrosis factor (TNF) receptor superfamily member 8 (TNFRSF8) and guanine nucleotide exchange factor vav3 (VAV3). These are the single genes that were most clearly associated with vitamin D status within the significant gene sets. All of the seven genes were associated with vitamin D deficiency, that is, more highly expressed in the deficient group compared with the sufficient group.

Figure 1

Clustering graphs and dendrograms of the two most significant gene sets, showing significant single genes associated with vitamin D status. Adjusted P<0.05 is indicated by the solid black line in the clustering graph. Vitamin D status was measured as plasma 25(OH)D concentration, and defined as sufficient (25(OH)D>50 nmol/l) or deficient (25(OH)D<37.5 nmol/l). (a) In gene set 11, two genes (CD14 and VAV3) were significantly associated with vitamin D deficiency (dark gray). The gene set was based on validated differentially expressed genes after 1,25(OH)2D treatment of prostate epithelial cells by Kovalenko et al.[25] (b) In gene set 19, two single genes (CD14 and THBD) and a group of four genes were significantly associated with vitamin D deficiency (dark gray). The gene set was based on upregulated genes after 1,25(OH)2D treatment of monocytic leukemia cells by Heikkinen et al.[8] IFNGR2, interferon-γ receptor 2; LRG1, leucine-rich-α-2-glycoprotein 1; PFKFB4, 6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 4; ZFP36, zinc finger protein 36.

Discussion

Considering the large number of VDR-responsive genes and the myriad of downstream pathways, genome-wide approaches are necessary when investigating transcriptional effects of vitamin D. We explored the association of vitamin D status (measured as plasma 25(OH)D concentration) and whole blood genome-wide expression profiles in a cross-section of middle-aged Norwegian women (n=218). Using a targeted, pathway-level analysis, we identified 26 gene sets differentially expressed when comparing sufficient versus deficient vitamin D status. A majority of the gene sets were related to immune system regulation, anti-inflammatory mechanisms and modulation of immune cell function and signaling. To our knowledge, the data presented here provides the first insight into the impact of vitamin D status on genome-wide transcription in the blood, using a representative free-living population.

Overarching gene sets such as signaling in immune system (gene set 27), innate immunity signaling (gene set 30), cytokine production (gene set 52) and chemokine signaling (gene set 49) were differentially expressed according to vitamin D status. The majority of genes in these gene sets were associated with vitamin D deficiency, in line with the hypothesis that vitamin D limits pathological immune responses that may ultimately lead to hypersensitivity or autoimmunity.[6, 22] Three pathways related to innate immunity were more highly expressed in the vitamin D-deficient group: TLR signaling (gene set 64 and 65), a major pathway governing the inflammatory response to infection, and IL-1R pathway (gene set 55), which increases migration of leukocytes to sites of infection. In vitro studies have shown that 1,25(OH)2D or vitamin D analogs have anti-inflammatory properties by downregulating inflammatory gene expression in monocytes/macrophages. In monocytes, the downregulation of TLR2 by 1,25(OH)2D in a time- and dose-dependent manner has been demonstrated by Sadeghi et al.[23] and others,[22] and TLR2 was associated with vitamin D deficiency in gene set 65. Contrasting the findings by Sadeghi et al.[23] and other in vitro studies,[8, 24, 25] we found that CD14 significantly contributed to the results of several gene sets (gene set 1, 3, 10, 11, 19, 22, 64, 65) and was associated with vitamin D deficiency. CD14 is a pattern recognition receptor for lipopolysaccharides and other pathogen compounds, and it acts as a co-receptor with TLR4 to trigger cytokine production and inflammatory response to microbial infection. Its association with vitamin D deficiency in our data set is in line with the general anti-inflammatory actions of 1,25(OH)2D.

Several cellular signaling pathways that govern immune cell function and viability were associated with vitamin D status. In a non-neoplastic setting, processes such as regulation of cell proliferation and differentiation, induction of programmed cell death (apoptosis), and inhibition of angiogenesis and cell invasiveness are tightly regulated.[4] Anti-neoplastic effects of 1,25(OH)2D have been well documented, but the underlying processes are not fully understood.[4] VEGF signaling (gene set 44, associated with vitamin D deficiency) regulates angiogenesis, and was inhibited by 1,25(OH)2D treatment of cancer cells.[26] Activation of the VEGF receptor leads to signaling through a number of possible downstream pathways with a multitude of cellular effects besides angiogenesis, and p38/MAPK (gene set 70) is one of these pathways. Conversely, p38/MAPK may be activated by a number of signals in addition to VEGF, such as inflammatory and apoptotic factors. In line with our findings, in vitro studies have demonstrated modulation of the p38/MAPK pathway by 1,25(OH)2D in several cell types.[4, 27] Transcriptional induction by the p38/MAPK pathway may result in cytokine production, cell migration or apoptosis to regulate immune cell function. The differential expression of the death pathway (gene set 66) and/or the cytokine production pathway (gene set 52) may be related to changes in the major VEGF and p38/MAPK signaling cascades.

Ten significant gene sets were derived from experimental publications reporting gene expression changes in cells or animal models, after treatment with 1,25(OH)2D or vitamin D analogs. On the basis of analysis of primary osteoblast cultures derived from elderly, female bone fragment donors, Tarroni et al.[24] found several pathway-level changes that were significant in our data set, including TLR signaling (gene set 64 and 65), the p38/MAPK pathway (gene set 70) and the death receptor signaling pathway (gene set 66), all discussed above. In addition, their list of top differentially expressed genes was among our significant gene sets. The downregulated genes reported by Stuckey et al.[28] after treatment of xenograft mice with a vitamin D analog were significant (gene set 13 and 15) and mainly associated with vitamin D deficiency in our data. Three cell line studies[8, 25, 29] provided several significant gene sets, however, the direction of gene expression change was not always in consistency with our findings. This illustrates that in vitro results may give valuable indications concerning vitamin D function at the population level. Still, several gene sets did not come out significant, leading us to conclude that care should be taken when attempting to extrapolate in vitro findings to hypotheses about in vivo function. Discrepancies between reported findings may be related to the differential response of specific cells and tissues, treatment dose and duration, differing methodological and analytical methods, and bias introduced by the authors when reporting only subsets of differentially expressed genes.

The cross-sectional study design allows no conclusions to be made concerning causal relationships, and because only post-menopausal women were included, generalizability to other population strata is limited. Our study failed to identify expression changes at the single-gene level. This was not unexpected, and most likely due to our study design that focused on population-level associations. In the multi-factorial setting of a free-living population, the impact of single nutrients such as vitamin D is not easily detected. Gene expression differences would perhaps be larger if the 25(OH)D concentrations had a wider distribution than in this study population, for example, in an intervention study. Owing to the nonsignificant differences at the single-gene level, and the high correlations of ABI results and quantitative reverse transcriptase PCR (qRT–PCR),[30] validation by qRT–PCR was not performed.

The negative result of the overall test of global significance points to inter-individual variation of gene expression profiles as a main contributor to the lack of clear associations.[18] Variation may stem from dissimilar levels of blood cell sub-populations, of which we had no available information. Residual bias may remain in the data set after adjustment for significant technical and biological factors, for example, related to vitamin A, antioxidant levels or unknown factors. Finally, when using the Global Test, small gene sets may be overrepresented among significant results, as may be the case with gene set 11. Overall, sizes of the significant gene sets vary from 4 to 292 (mean: 65, median: 40), giving no indication of overrepresentation of small gene sets.

We conclude that vitamin D status was associated with pathways related to immune system and immune cell function, as well as major signaling cascades. Furthermore, gene sets reported from experimental in vitro studies of the biologically active 1,25(OH)2D or vitamin D analogs were identified at the population level. This study demonstrates the feasibility of blood gene expression profiling to explore effects of nutritional factors such as vitamin D in the general population. Studies using similar methods are likely to yield increased insight into the molecular effects of nutrients.