Search tips
Search criteria 


Logo of narLink to Publisher's site
Nucleic Acids Res. 2011 November; 39(21): 9181–9193.
Published online 2011 August 16. doi:  10.1093/nar/gkr654
PMCID: PMC3241659

Nuclear hormone 1α,25-dihydroxyvitamin D3 elicits a genome-wide shift in the locations of VDR chromatin occupancy


A global understanding of the actions of the nuclear hormone 1α,25-dihydroxyvitamin D3 (1α,25(OH)2D3) and its vitamin D receptor (VDR) requires a genome-wide analysis of VDR binding sites. In THP-1 human monocytic leukemia cells we identified by ChIP-seq 2340 VDR binding locations, of which 1171 and 520 occurred uniquely with and without 1α,25(OH)2D3 treatment, respectively, while 649 were common. De novo identified direct repeat spaced by 3 nucleotides (DR3)-type response elements (REs) were strongly associated with the ligand-responsiveness of VDR occupation. Only 20% of the VDR peaks diminishing most after ligand treatment have a DR3-type RE, in contrast to 90% for the most growing peaks. Ligand treatment revealed 638 1α,25(OH)2D3 target genes enriched in gene ontology categories associated with immunity and signaling. From the 408 upregulated genes, 72% showed VDR binding within 400 kb of their transcription start sites (TSSs), while this applied only for 43% of the 230 downregulated genes. The VDR loci showed considerable variation in gene regulatory scenarios ranging from a single VDR location near the target gene TSS to very complex clusters of multiple VDR locations and target genes. In conclusion, ligand binding shifts the locations of VDR occupation to DR3-type REs that surround its target genes and occur in a large variety of regulatory constellations.


The vitamin D receptor (VDR) is an endocrine member of the nuclear receptor superfamily, since it is activated already by sub-nanomolar concentrations of its natural ligand 1α,25-dihydroxyvitamin D3 (1α,25(OH)2D3) (1). The classical, physiological role of 1α,25(OH)2D3 is the regulation of calcium and phosphate homeostasis and bone mineralization (2), but there is both epidemiological and pre-clinical evidence that VDR ligands have also anti-proliferative and immuno-modulatory actions (3,4). VDR is nearly ubiquitously expressed (5) and many cell types [for example, bone, skin and monocytes (6)] are capable of metabolizing the main circulating form of vitamin D, 25-hydroxyvitamin D3 (25(OH)2D3), to 1α,25(OH)2D3 by the enzyme 25(OH)D3-1α-hydroxylase, encoded by the gene CYP27B1, i.e. 1α,25(OH)2D3 can have in these tissues both auto- and paracrine roles.

In recent years, an increasing amount of data has shown a physiological function of 1α,25(OH)2D3 in the immune system, both in defense and in maintenance of self-tolerance (7). VDR ligands affect the immune system in many different ways, but the most important seems to be the enhanced capacity of innate anti-bacterial defense and a more tolerogenic profile toward autoimmune phenomena (4). One important target cell of 1α,25(OH)2D3 is the monocyte/macrophage, where the nuclear hormone promotes the anti-bacterial defense system by inducing the secretion of anti-bacterial proteins, such as cathelicidin anti-microbial peptide (CAMP) (8). Importantly, decreased anti-bacterial and anti-viral defenses are obvious in vitamin D deficiency, and restoration of the vitamin D-sufficient state restores immune responses (4).

As a member of the nuclear receptor superfamily VDR acts as a transcription factor that binds to specific vitamin D response elements (VDREs) within the regulatory regions of its primary target genes (1). VDR forms a heterodimer with retinoid X receptor (RXR), another member of the nuclear receptor superfamily, and binds preferentially to VDREs that are formed as direct repeat of two hexameric core binding motifs with 3 nucleotides spacing, so-called DR3-type VDREs (9). However, also other types of VDREs, such as direct repeats with 4 intervening nucleotides (DR4) (10) or everted repeats (ERs) with a spacer of 6–9 nucleotides (ER6-9) (11,12) have been reported. Single gene analyses suggest that VDR target genes often contain multiple functional VDREs (13–16).

VDR belongs to those nuclear receptors that bind to their genomic targets already in the absence of ligand. Under this constellation VDR has been shown to associate via co-repressor proteins with histone deacetylases and to repress its target genes (17). Ligand binding induces a conformational change that leads to dissociation of co-repressors and to the association of different types of co-activator proteins that on one hand lead to chromatin opening and on the other hand are part of the mediator complex, which forms a bridge to the basal transcription machinery associated with the transcription start sites (TSSs) of primary VDR target genes (18).

Thus far, both single gene and microarray analyses in various tissues and cells from different species have suggested a long list of VDR target genes. However, the actions of VDR seem to be highly cell-specific, since there is only very little overlap of VDR targets between the different tissues (19–22). Moreover, such studies have often been performed at such late time points that they did not allow distinguishing primary VDR targets from secondary ones. Recently, Ramagopalan et al. (23) reported the first genome-wide view on VDR locations in human lymphoblastoid cell lines suggesting 2776 genomic VDR locations after ligand treatment and 229 target genes. However, also these data were obtained from a very late time point (36 h after ligand stimulation) and therefore may not fully represent the primary actions of VDR.

In this study, we used undifferentiated THP-1 human monocytic leukemia cells as a model system, since they reflect reasonably well the 1α,25(OH)2D3 response of primary human monocytes (24,25). To gain a genome-wide view of the primary actions of VDR in these cells, we determined its genomic binding sites after 40 min stimulation with 1α,25(OH)2D3 and in the unstimulated state and compared it with the VDR target genes after 4 h ligand treatment. We chose human monocytes as cellular model out of interest in the biological role of 1α,25(OH)2D3 signaling in the initial events of immune response. However, the main focus of this study is of mechanistic nature showing that VDR is associated with chromatin even in the absence of ligand, but that the balance of VDR binding shifts from locations that are gene-centered and lacking classical VDREs to locations that are more distal to the genes and contain one or more DR3-type response elements (REs).


Cell culture

THP-1 human monocytic leukemia cells were grown in RPMI 1640 medium supplemented with 10% fetal calf serum, 2 mM l-glutamine, 0.1 mg/ml streptomycin and 100 U/ml penicillin and the cells were kept at 37°C in a humidified 95% air/5% CO2 incubator. Prior to mRNA or chromatin extraction, cells were grown overnight in phenol red-free medium supplemented with charcoal-stripped fetal bovine serum. Then, for RNA extractions cells were treated for 4 h with 1α,25(OH)2D3 (Sigma-Aldrich) or with solvent (ethanol, finally 0.1%) or for chromatin immunoprecipitation (ChIP) assays for 40 min with 1α,25(OH)2D3 or left unstimulated.

RNA extraction and cDNA synthesis

Total RNA was extracted using the Mini RNA Isolation II Kit (Zymo Research, HiSS Diagnostics, Freiburg, Germany). cDNA synthesis was performed for 60 min at 37°C using 1 µg of total RNA as a template, 100 pmol oligo(dT)18 primers (Eurogentec, Seraing, Belgium), 500 µM dNTPs, 40 U Ribolock Ribonuclease Inhibitor and 40 U MMuLV reverse transcriptase (Fermentas, St Leon-Rot, Germany). Prior to real-time quantitative PCR (qPCR) reaction, the cDNA was diluted 10-fold.


qPCR was performed using an ABI 7500 Fast System (Applied Biosystems, Lennik, Belgium). The reactions were performed using 250 nM of reverse and forward primers (for sequence see Supplementary Table S1), 4 µl 1/10 diluted cDNA template and the ABsolute Q-PCR SYBR Green Low ROX Mix (Abgene, Westburg, Leusden, The Netherlands) in a total volume of 20 µl. In the PCR reaction Hotstart Taq polymerase was activated for 15 min at 95°C, followed by 40 amplification cycles of 30 s denaturation at 95°C, 30 s annealing at primer-specific temperatures (Supplementary Table S1) and 40 s elongation at 72°C and a final elongation for 5 min at 72°C. PCR product specificity was monitored using post-PCR melt curve analysis. Relative expression levels were determined with the comparative delta threshold cycle (ΔCt) method. Amplification efficiencies of the primer pairs were taken into account and have been determined on the basis of a standard curve of a cDNA dilution series (with two exceptions, all PCR efficiencies were between 80% and 120%). Relative expression levels of the target genes were normalized to the three most stable of 10 tested internal control genes (b2M, GAPDH, HPRT1). The expression stability of the control genes was determined by the geNorm algorithm (26). Briefly, the arithmetic mean of replicated Ct values for each gene is transformed to a relative quantity (setting the sample with the highest expression as calibrator to 1), using the ΔCt formula Q = EΔCt = E(calibratorCt  sampleCt) (Q = quantity sample relative to calibrator sample; E = amplification efficiency). For normalization, the relative quantities were divided by the normalization factor being the geometric mean of the three reference genes.

Microarray analysis

Total RNA was checked for RNA integrity using a Biorad Experion electrophoresis station (Biorad, Nazareth, Belgium) and analyzed on Sentrix Human-6 v2 Expression BeadChips from Illumina (San Diego, CA, USA) at the Finnish Microarray Centre (Turku, Finland) using protocols recommended by the manufacturer. Raw microarray data is available at Gene Expression Omnibus (GEO) under accession number GSE27270 (SubSeries of GSE27438). Microarray data analysis was performed using R statistical software version 2.11 (R Dev Core) with associated libraries from Bioconductor project version 2.6 (27). Data was normalized using VST transformation and RSN normalization used as standard approach for Illumina arrays with the lumi package (28). Probe sets that were not linked to any known or predicted human gene were filtered out. Linear Models for Microarray Data (limma) package (29) using linear model fitting for statistical testing with empirical Bayes variance smoothing procedure was applied to detection of differentially expressed genes. Obtained P-values were corrected for multiple testing using Benjamini–Hochberg false discovery rate (FDR) procedure (30). To identify the biological processes targeted by the 1α,25(OH)2D3-treatment, we used GeneTrail (31) to analyze the enrichment of up and downregulated VDR target genes into GO categories using default settings.


After ligand or no treatment of THP-1 cells, nuclear proteins were cross-linked to DNA by adding formaldehyde directly to the medium to a final concentration of 1% and incubating at room temperature for 15 min on a rocking platform. Cross-linking was stopped by adding glycine to a final concentration of 0.125 M and incubating at room temperature for 5 min on a rocking platform. The cells were collected by centrifugation and washed twice with ice cold PBS (140 mM NaCl, 2.7 mM KCl, 1.5 mM KH2PO4, 8.1 mM Na2HPO4·2H2O). The cell pellets were resuspended in 1 ml of lysis buffer (1% SDS, 10 mM EDTA, protease inhibitors, 50 mM Tris–HCl, pH 8.1) and the lysates were sonicated to result in DNA fragments of 200–400 bp. Cellular debris was removed by centrifugation. Aliquots of 100 µl of the lysate were diluted 1:10 in ChIP dilution buffer (0.01% SDS, 1.1% Triton X-100, 1.2 mM EDTA, 167 mM NaCl, protease inhibitors, 16.7 mM Tris–HCl, pH 8.1) and 1 µg of anti-VDR antibody (ab3508-100, Abcam, Cambridge, UK) or non-specific IgG (12-370, Millipore, Temecula, CA, USA) were added and the samples were incubated for overnight at 4°C on a rotating platform. The immunocomplexes were collected using 20 µl of Magna ChIP protein A magnetic beads (Millipore) at 4°C for 1 h with rotation. The beads were washed sequentially for 5 min in rotating platform with 1 ml of the following buffers: low salt wash buffer (0.1% SDS, 1% Triton X-100, 2 mM EDTA, 150 mM NaCl, 20 mM Tris–HCl, pH 8.1), high salt wash buffer (0.1% SDS, 1% Triton X-100, 2 mM EDTA, 500 mM NaCl, 20 mM Tris–HCl, pH 8.1) and LiCl wash buffer (0.25 M LiCl, 1% Nonidet P-40, 1% sodium deoxycholate, 1 mM EDTA, 10 mM Tris–HCl, pH 8.1). Finally, the beads were washed twice with 1 ml TE buffer (1 mM EDTA, 10 mM Tris–HCl, pH 8.0) and the immune complexes were eluted twice using 250 µl elution buffer (1% SDS, 100 mM NaHCO3) at room temperature for 15 min with rotation. The supernatants were combined and the immune complexes were reverse cross-linked at 64°C overnight in the presence of proteinase K (Fermentas) in a final concentration of 40 µg/ml. DNA was extracted with phenol/chloroform/isoamyl alcohol (25/24/1) and precipitated with 0.1 volumes of 3 M sodium acetate, pH 5.2 and 2 volumes of ethanol using glycogen as carrier. For ChIP-seq five individual VDR ChIP samples were pooled and concentrated using a SpeedVac. The ChIP templates were sequenced using a Solexa Gene Analyzer II platform (Illumina) at 36 bp read length using standard manufacturer protocols at the Genomics Core Facility in Heidelberg, Germany. For the validation of the ChIP-seq data, selected genomic regions containing VDR peaks were analyzed by quantitative ChIP analysis of independent samples using Maxima™ Sybr Green qPCR Master Mix (Fermentas) and the genomic primers listed in Supplementary Table S2. In addition to VDR and IgG antibodies, immunoprecipitation was performed using 1 µg of anti-RXR antibody (sc-553, Santa Cruz Biotechnologies, Santa Cruz, CA, USA). The qPCR reactions were performed with a LightCycler 480 apparatus (Roche) using the following PCR profile: 10 min at 95°C, followed by 50 cycles of 30 s at 95°C, 30 s at 58°C and 30 s at 72°C. The results were normalized with respect to input. Non-specific IgG results were subtracted by using the formula 2−(ΔCp)*100(specific antibody)  2−(ΔCp)*100(non-specific IgG), where ΔCp is Cp(immunoprecipitated DNA)  Cp(input) and Cp is the fractional cycle number. Statistical significance was calculated in reference to non-treated sample using two-tailed Student’s t-test.

ChIP-seq data analysis

Alignment of sequence reads against the human reference genome version hg19 was done using Bowtie software version 0.12.2. (32) with the following command line arguments: bowtie -n 1 -m 1 -e 70 -l 28 -k 1 -t -p 8 -q -S --best hg19 input_file_name output_file_name. Of note, one Solexa run with 1, 2 and 1 lanes of IgG, unstimulated and ligand-treated VDR ChIP sample sequencing, respectively, suffered from a minor post-sequencing processing error at the sequencing platform that had resulted in an omission of the 16th nucleotide from the large majority of reads. This was corrected for by inserting an N at the missing position in each read and using two modified settings in the Bowtie alignment (-n 2 -m 3). The short read data publically available under accession number GSE27437 (SubSeries of GSE27438) contains the corrected reads. Aligned reads from multiple lanes were combined per sample. In the case of ChIP-seq analysis, the common IgG background read set was scaled so that an individual background was generated for each VDR ChIP-seq sample by randomly removing IgG reads such that the read ratio IgG/VDR became the same for each VDR ChIP-seq sample. The scaled IgG and un-scaled VDR reads were then converted to sorted BAM format using SAM tools (33) and subsequently also to TDF format using igvtools to allow visualization in the IGV genome browser (34).

For ChIP-seq analysis statistically significant peaks were identified using MACS program version (35) with the following command line arguments: macs --pvalue=1e-2 --nomodel --wig -t input_sample_file_name -c input_control_file_name --tsize=36 --format=BAM --name=analysis_topic --mfold=13 --shiftsize=250 --bw=250 --verbose=3. A loose initial P-value cutoff was purposeful used to enable catching all true signals with less regard to possible false positive signals at this early phase of analysis. Subsequent refinement of MACS peaks was done using PeakSplitter program version 1.0 (36) with the following command line arguments: PeakSplitter -p peak_folder_name -w aligned_read_wig_folder_name -o output_folder -c 5 -v 0.6 -f. An in-house R script was used to calculate fold enrichments (FEs) over IgG, P-values and FDR estimates for the subpeaks in an identical manner used internally in MACS Since splitting the original peaks into several subpeaks may produce large sets of weaker flanking residual peaks, thus biasing for example the genomic and FE distributions of the peaks, we kept for further analysis only those peaks that either had FDR <1% or were the best subpeak (by FDR) within their parental MACS peak.

Further data analysis was conducted using an in-house R pipeline containing tools for identifying peaks overlapping in two samples, the analysis of genomic distribution of peaks and the integration of peak and gene expression data sets. Of note, since there are genomic positions that cannot be uniquely assigned to the selected 10 types of genomic elements used in the analysis of peak distributions therein, a prioritization scheme was employed, where the peaks were uniquely overlapped to the elements in a step-wise exclusive scheme, starting from coding elements (UTRs and coding exons), and then moving to introns and outward from the gene.

Detection of transcription factor binding sites

For the de novo identification of the VDR binding consensus sequence within ±100 bp of the summits of all 2340 FDR <1% VDR peaks we employed the rGADEM R-package with the following settings: GADEM(Sequences, genome=Hsapiens, weightType=1, populationSize=500, numGeneration=20). Of note, the default settings did not identify the DR3-like sequence in the full peak set. The resulting candidate position weight matrices (PWMs) were then compared to the set of known PWMs in the JASPAR database using the MotIV R-package, and their locations mapped back to the original peaks using in-house R-scripts. To scan the peak sequences for PWMs representing different types of DRs and ERs with varying spacer lengths, we used ‘matrix-scan (quick)’ that is one of the regulatory sequence analysis tools [RSAT, available at (37)]. For the RSAT scans, background was modeled at Markov order 2 from the ‘upstream-noorf’ subset of the human genome and a P-value cutoff of 1 × 10−4 was used.


Identification of VDR target genes

To identify VDR target genes, THP-1 cells were treated for 4 h with either vehicle (ethanol) or 1α,25(OH)2D3 and collected for the extraction of total RNA that was subjected to gene expression microarray analysis. We chose the 4 h time point as a compromise to obtain reasonable mRNA up or downregulation of VDR target genes, the majority of which still can be considered as primary targets. Statistically significant effect on gene expression levels (adjusted P < 0.05) was detected for transcripts representing 638 separate Ensembl genes (604 separate RefSeq gene symbols), out of which 408 (384) were upregulated and 230 (221) downregulated (Supplementary Table S3). The broad expression range of these 638 genes is indicated in Supplementary Figure S1A. In addition to the larger number of affected genes, also the magnitude of change was higher for the upregulated set of VDR target genes (maximal upregulation 6.65-fold), compared to the downregulated set (maximal downregulation 1.48-fold). Validation of 43 selected genes on independent samples by qPCR showed a high correlation (r = 0.889, P < 1.8 × 10−14) to the microarray results across a wide range of fold changes (Supplementary FigureS1B). The validation results of five upregulated immunity-related genes CAMP, CD14, SP100, LRRCA8 and CD93 and a downregulated immunity-related gene HSH2D are shown in Supplementary Figure S1C.

Gene ontology (GO) enrichment analysis of the 408 upregulated VDR target genes using GeneTrail (31) highlighted the role of 1α,25(OH)2D3 in the regulation of immune response and signaling cascades and in the positive regulation of higher-level processes (Table 1). In particular, a positive effect of 1α,25(OH)2D3 on the production and secretion of the cytokine interleukin (IL) 1 was implicated. The 230 downregulated genes showed association with the negative regulation of many similar higher-level processes as the upregulated genes. In contrast to upregulated genes, GO terms enriched for downregulated genes lacked categories related to immunity and signaling, but interestingly included a category for transcription regulator activity (Table 1). Combining the up and downregulated genes for the enrichment analysis did not point to additional functions (data not shown).

Table 1.
GO analysis of VDR target genes

In summary, we identified in THP-1 cells a set of 638 genes as primary targets of regulation by 1α,25(OH)2 D3, and thus by VDR, that were enriched for genes that are associated with the GO terms related to immunity and signaling.

Genome-wide mapping of VDR binding sites

To identify genomic VDR binding sites, THP-1 cells were treated for 40 min with 1α,25(OH)2D3 or remained unstimulated and ChIP assays were performed with antibodies against VDR and non-specific IgGs. Purified chromatin samples were subjected to deep sequencing using a Solexa GAII platform. Analysis of the resulting sequences, using IgG-precipitated sample as a negative control, revealed the genome-wide landscape of the chromatin occupation of VDR in THP-1 cells. The reads at various stages of the analysis are summarized in Supplementary Table S4. In the 1α,25(OH)2D3-treated sample the number of high-confidence (FDR <1%) peaks indicating genomic VDR binding sites was 1820, in contrast to 1169 sites in the unstimulated sample (Figure 1A). The total number of VDR binding sites was 2340, with 649 (27.7%) sites common in both unstimulated and 1α,25(OH)2D3-treated samples and 520 (22.2%) and 1171 (50.0%) unique in unstimulated and 1α,25(OH)2D3-treated samples, respectively (Figure 1A). The genomic locations of the 2340 binding sites are listed in Supplementary Table S5.

Figure 1.
VDR binding locations and genomic distribution in unstimulated and 1α,25(OH)2D3-treated THP-1 cells. ChIP-Seq was performed from chromatin templates extracted from THP-1 cells that were either unstimulated or treated for 40 min with 10 nM ...

VDR binding locations, grouped to unique and common as in the Venn diagram in Figure 1A, showed differences in their distribution to various genomic elements (Figure 1B). Within both the unstimulated and 1α,25(OH)2D3-treated peak sets, the high-confidence peaks, compared to peaks with lower confidence, were enriched to genomic regions representing proximal promoter and exons, and were depleted from distal regions and introns other than the first (Supplementary Figure S2). However, comparing the distributions of unstimulated and 1α,25(OH)2D3-treated peak sets revealed that upon ligand treatment VDR favors significantly more such binding locations that are distal to genes (>1 kb upstream or >10 kb downstream) and disfavors those representing coding sequences (Figure 1B). The shift of VDR binding from proximal to distal locations upon ligand treatment was also highlighted by the VDR peak proximity to TSS (Figure 1C). The density of VDR binding proximal to the TSS diminished with decreasing peak confidence, i.e. with increasing FDR (peaks with FDR >10% may be considered random-like), for peaks that were either unique in the unstimulated state (left panel) or common in both treatments (middle panel). In contrast, such enrichment was nearly completely lacking in the 1α,25(OH)2D3-treated cells (right panel). Together with Figure 1B, this suggests a non-random enrichment of VDR binding locations to gene-dense regions which, at first glance, appears weaker for the 1α,25(OH)2D3-treated peak set but may in fact indicate enrichment to a subset of genes, such as regulated genes.

Taken together, we have identified in THP-1 cells 2340 genomic VDR binding locations, a surprisingly large proportion of which display at least some binding also in the unstimulated state and not only upon 1α,25(OH)2D3-treatment.

Transcription factor binding sequences at genomic VDR locations

De novo binding site search using rGADEM and MotIV R/Bioconductor packages identified the classical DR3-type RE consensus sequence for VDR–RXR heterodimers as the most highly enriched within ±100 bp of the VDR peak summits (E-value 2.65 × 10−12 for the JASPAR consensus matrix MA0074.1 for RXRA::VDR). Of note, while the JASPAR consensus has near-exclusive T at positions 4, 12 and 13, more variability is allowed at these positions in the de novo DR3-type REs actually occupied by VDR in THP-1 cells (Figure 2A). However, only 31.7% (742) of all VDR peak summits included one or more de novo DR3-type RE. The quality of the identified de novo DR3-type REs appears high, since the more each of the de novo DR3-type REs resembled their common consensus, i.e. had lower E-value, the higher was the FE of the respective VDR peak (Pearson r = −0.217, P < 3 × 10−15, Figure 2B). Hoping to identify any atypical VDR binding consensus sequences for the remaining 68.3% (1598) of VDR peaks, we re-ran the de novo analysis on them only, but even the most frequent consensus nrGGrGswGGrrn, corresponding best to a SP1 binding site in the JASPAR database (MA0079.2, E-value 1.03 × 10−7), was found in 538 peaks (23.0%). Another consensus sequence, nsAGGAAGn, was even less frequent with 274 peaks (11.7%). The latter had high similarity to the consensus for ETS transcription factor family member SPI1/PU.1 (MA0080.2; E-value 8.50 × 10−7). De novo analysis of even more focused subgroups, like the peaks unique to unstimulated cells, did not reveal additional consensus sequences (data not shown). Using a set of DR- and ER-type VDREs where the half sites were based on the identified de novo DR3-type RE and which covered all spacer lengths from 0 to 16 nucleotides, VDR peak sequences lacking the de novo DR3-type RE were scanned using RSAT (37). This scan failed to identify any VDRE with statistically significant enrichment to ±100 bp of the peak summit versus the two 150 bp flanks of that summit region (data not shown). It should be noted that we found no genome-wide evidence that VDR binds to any other type of VDRE than DR3. However, the possibility can not be excluded that a few individual non-DR3 sites exist that bind VDR.

Figure 2.
De novo identified DR3-type REs. De novo analysis was performed on sequences ±100 bp of the FDR <1% VDR peak summits. (A) The consensus DR3-type RE identified by the rGADEM de novo search is displayed in comparison to the public ...

In summary, only about a third of the genomic VDR binding sites in THP-1 cells contained a DR3-type RE and no common element was found to consistently explain the presence of VDR at locations lacking the classical VDRE.

Ligand-induced genome-wide shift in DR3-type RE usage by VDR

While the de novo DR3-type RE was, at first look, apparently relatively infrequent at sites of VDR chromatin occupation in THP-1 cells (Figure 2), on closer inspection its presence strikingly correlated with the 1α,25(OH)2D3-dependent effects on VDR binding at each location (Figure 3A). The de novo DR3-type RE was present in about 20% of genomic VDR locations where the VDR binding in the unstimulated sample exceeded that in the 1α,25(OH)2D3-treated sample by >10 FE. Conversely, when the VDR binding after a short 40 min 1α,25(OH)2D3 treatment superseded that of the unstimulated state, the percentage of peak locations with a DR3-type RE climbed to >90%. Furthermore, the number of de novo DR3-type REs near the peak summit tended to associate with increased peak FE within the common and 1α,25(OH)2D3-exclusive peaks, but not within the unstimulated-exclusive peaks (Figure 3B). The presence or absence of either the de novo SPI1/PU.1 or SP1-like consensus sequences, irrespective of coinciding de novo DR3-type REs, did not affect the FE of VDR peaks (data not shown).

Figure 3.
Ligand-dependent VDR occupation on the de novo DR3-type REs. The usage of the de novo DR3-type RE is strongly ligand dependent. (A) The running mean of the difference in peak FE at VDR binding locations (x-axis; 1α,25(OH)2D3 minus unstimulated, ...

Taken together, even without the ligand VDR occupies many chromatin sites, but most of these do not contain a DR3-type RE. In contrast, even a short ligand-treatment focuses the bulk of VDR binding to sites that do contain at least one DR3-type RE.

Proximity of VDR binding sites to 1α,25(OH)2D3 target genes

Combining gene expression data with the genomic locations of VDR offered an avenue to explore the mechanisms of target gene regulation by VDR. First we mapped the VDR binding locations into the neighborhoods of each gene regulated by the 4 h 1α,25(OH)2D3 treatment. For the 408 upregulated genes the bulk of VDR regulation seems to occur via sites within only 30 kb from the target gene TSS. However, a considerable proportion of regulation appears to be orchestrated by VDR binding sites up to 400 kb away, especially when the sites are unique to the 1α,25(OH)2D3 treatment (Figure 4A, top panel and Supplementary Figure S3A, top right panel). In contrast, the 230 downregulated genes had only slight enrichment for common and 1α,25(OH)2D3-unique peaks near the TSS (Figures 4A, Supplementary Figures S3A and S3B, bottom panels). The frequency of VDR peaks around the TSS was also higher for upregulated than for downregulated genes (2.04 ± 2.54 versus 1.05 ± 2.16 peaks per gene, Wilcoxon rank sum test P-value 5.01 × 10−13). The strong enrichment of 1α,25(OH)2D3-treated peaks to the 1α,25(OH)2D3-regulated gene neighborhoods confirms regulated genes as the likeliest source to the apparent depletion of these peaks from the proximity of any TSS seen in Figure 1B and and1C.1C. Thus, VDR seems to shift from sites in the vicinity of the TSS of non-regulated genes to the proximity of VDR-regulated genes.

Figure 4.
VDR peak height, location at <400 kb from the TSS, and the presence of de novo DR3-type REs predicts strong upregulation of 1α,25(OH)2D3 target genes. (A) Common and treatment-unique VDR binding sites in the neighborhoods of up- ...

Since the presence of the de novo DR3-type REs strongly associated with ligand-dependency of the VDR binding (Figure 3), we investigated whether this is also reflected by the VDR target gene regulation (Figure 4B). We identified for each regulated gene the strongest VDR binding location in 1α,25(OH)2D3-treated cells within 400 kb of the gene TSS and defined them as the most likely regulating VDR locations (Supplementary Table S6). These sites were further split to those in proximal (0 to ± 30 kb of the TSS) and distal regions (±30 to ±400 kb). If no VDR peak was present within the 400 kb limit, we used from the outlying ‘desert’ the closest peak of any FE to represent a potential very distant regulator. Genes were divided, separately for both up and downregulated sets, into equally sized tertiles according to their stimulated fold change [4 h 1α,25(OH)2D3 versus vehicle treatment]. The most likely regulating VDR locations were also divided at median FE to equally sized ‘low’ and ‘high’ VDR peak sets and further to peaks with or without any DR3-type RE near the peak summit. Gene counts in these 12 sub-categories suggest that the higher the VDR peak is, the more likely it is to strongly upregulate its target gene, especially when it contains a DR3-type RE. As already seen in Figure 4A and Supplementary Figure S3, strong regulation by VDR is possible from both proximal and distal sites, but less likely from sites beyond ± 400 kb. Also among the genes for which the most likely regulating VDRE was in the ‘low’ category, the proportion of DR3-type RE containing peaks versus peaks without such sites increased toward the high fold change genes. In contrast to upregulated genes, the VDR peak FE or the presence of DR3-type REs associated much less with the fold change of downregulated target genes.

Taken together, 72% of the 408 upregulated 1α,25(OH)2D3 target genes showed VDR binding within 400 kb of their TSS, while this applied only for 43% of the 230 downregulated genes.

Examples of VDR target genes

Examples of previously unidentified VDR target genes that were upregulated by a 4 h 1α,25(OH)2D3 treatment, and for which the most likely regulating VDRE could be easily proposed, are shown in Figure 5. SP100, a gene encoding SP100 nuclear antigen that has been identified as a co-repressor for transcription factors, such as ETS1 (38), is characterized by a single, strong and ligand-dependent VDR binding region with a single de novo DR3-type RE located in the first intron at + 0.5 kb of the gene’s TSS (a category 4 peak, Figure 5A).

Figure 5.
VDR ChIP-seq analysis of genomic loci for the 1α,25(OH)2D3 target genes SP100, ELL and THBD. In (A–C), the peak tracks show the unstimulated and 1α,25(OH)2D3-treated VDR peaks. The high-confidence peaks are indicated as yellow-to-orange-to-red ...

Also for ELL, a gene encoding for an elongation factor of RNA polymerase II that has also been identified as a transcriptional co-factor for transcription factors, such as HIF-1 (39), the single dominant VDRE is located in the first intron at +19.3 kb of the gene’s TSS, is slightly ligand-dependent although with considerable basal binding and contains a DR3-type RE (Figure 5B). The locus carries also a distal 1α,25(OH)2D3 target gene, LRRC25 encoding the protein leucine rich repeat containing 25 of unknown function, for which the ELL VDR binding region, located 105 kb upstream of the LRRC25 TSS, appears the most likely regulator. Therefore, for the ELL gene the VDR peak is of category 4, while for the LRRC25 gene it is of category 5.

The third example is a complex locus containing the three 1α,25(OH)2D3 target genes thrombomodulin (THBD), encoding for the integral membrane protein THBD (also called CD141) (40), CD93, encoding for the C-type transmembrane receptor CD93 (41), and GZF1, encoding for the GDNF-inducible zinc finger protein 1 (42) (Figure 5C). However, the microarray result for GZF1 could not be validated by qPCR. The THBD gene responds most effectively to 1α,25(OH)2D3 treatment, while the strongest ligand-dependent VDR binding region of the genomic locus is 285 kb upstream of the THBD TSS and contains eight distinct DR3-type REs and, i.e. it is a category 5 peak. There is also one other good genomic VDR location at 176 kb upstream of the THBD TSS, but this lacks a DR3-type RE. The locus has still two additional, but weaker VDR binding regions, one with a DR3-type RE and one without, at 77 and 107 kb upstream of the THBD TSS, respectively. Interestingly, there is also an unstimulated-unique VDR binding site within the short, single exon THBD gene. It is tempting to speculate that this VDR binding site could have a repressive role that could, together with the strong ligand-dependent VDR binding site in the vicinity, facilitate the relatively strong upregulation of THBD.

For all three examples the VDR binding at the sites indicated with an arrow was validated from independent chromatin samples of unstimulated and 40 min 1α,25(OH)2D3-treated THP-1 cells by qPCR (Figure 5D). In addition to VDR, also the presence of the heterodimerizing partner RXR was verified from the same chromatin samples. In all cases the 1α,25(OH)2D3 treatment significantly increased the amount of both VDR and RXR at the respective VDR binding region. In addition, the upregulation by 1α,25(OH)2D3 of SP100, ELL and THBD mRNA was confirmed on independent samples (Figure 5D).

Additional examples of VDR target gene and VDR binding site scenarios are shown in Supplementary Figure S4. A strong double VDR peak with each a DR3-type VDRE was found very close to the TSS of the uncharacterized VDR target gene in chromosome 20 called open reading frame 197 (C20orf197), which also seems to serve as regulatory site for the VDR target genes C20orf177 and PPP1R3D (Supplementary Figure S4A). The gene ArfGAP with SH3 domain, ankyrin repeat and PH domain (ASAP2) is regulated by three VDR binding sites, of which one is proximal and two are distal (Supplementary Figure S4B). The location of the myosin IXB (MYO9B) gene represents a similar complex locus as that of the THBD gene, as the MYO9B gene seems to be regulated together with the genes OCEL1 and ABHD8 by two proximal VDR binding sites that both contain a DR3-type RE and four proximal sites (Supplementary Figure S4C). As examples of downregulated VDR target genes the genes glucosaminyl (N-acetyl) transferase 1, core 2 (GCNT1) (Supplementary Figure S4D) and SUMO1 pseudogene 1 (SUMO1P1) (Supplementary Figure S4E) are shown having their next VDR binding site nearly 1 Mb and more than 300 kb distant from their TSS regions, respectively. Interestingly, the genomic region of SUMO1P1 also contains the 24-hydroxylase (CYP24A1) gene, which is in most cellular systems a strong VDR target gene with a number of ChIP-confirmed VDR binding locations (13,43).

In summary, the eight exemplary genomic loci showed the variation of gene regulatory scenarios ranging from a single VDR location close to the 1α,25(OH)2D3 target gene TSS, more complex ones with one VDR location distal to the TSSs of one or more 1α,25(OH)2D3 target genes and very complex ones with up to six VDR locations clustering with multiple 1α,25(OH)2D3 target genes.


For a general understanding of the mechanisms of 1α,25(OH)2D3 signaling it is essential to monitor the genome-wide location of VDR in relation to primary 1α,25(OH)2D3 target genes. In this study, we investigated the sites of VDR chromatin occupancy in THP-1 human monocytes after 40 min treatment with 1α,25(OH)2D3 and in the unstimulated state using ChIP-seq, and related these to the VDR target genes identified using microarrays from the same cellular model 4 h after ligand treatment. As expected, GO term analysis of the 638 primary 1α,25(OH)2D3 target genes that we identified in THP-1 cells showed a clear enrichment for immune-related genes. In particular, the positive effect of 1α,25(OH)2D3 on IL1 production and secretion was highlighted, which fits well with previous reports (44,45). However, we did not find IL1A or IL1B within our 638 1α,25(OH)2D3 target genes, but only the gene of the IL1 receptor antagonist (IL1RN) downregulated. The latter suggests that the regulation of IL1 activity by 1α,25(OH)2D3 is indirect via the block of the IL1 receptor. In addition, Kaler et al. (46) suggested that in THP-1 cells 1α,25(OH)2D3 downregulates the transcription factor STAT1, which is central in IL1 regulation.

Since the only nuclear target of 1α,25(OH)2D3 is VDR (47) and the 4 h treatment is a reasonably short time, most of the genes that we found in our microarray analysis are assumed to be primary VDR target genes, i.e. VDR should bind a regulatory chromatin region that loops to the target gene’s TSS region. The majority (64%) of the 1α,25(OH)2D3 target genes were upregulated, which is typical for most VDR target tissues (48,49).

Through the choice of early time points both in microarray and ChIP-seq experiments we aimed to detect mainly primary VDR effects. This is in contrast to the recent study of Ramagopalan et al. (23), who reported VDR ChIP-seq data in HapMap lymphoblastoid cells after a long (36 h) 1α,25(OH)2D3 treatment, focusing more on the association of VDR binding sites with disease SNPs and aspects of human evolution.

In this study, we identified 2340 genomic sites of VDR binding, 520 of which are exclusive to the unstimulated THP-1 cells. By size, our set of 1820 1α,25(OH)2D3-treated genomic VDR binding sites in THP-1 cells is comparable with the 2776 VDR binding sites in 1α,25(OH)2D3-treated lymphoblastoids (23), but the overlap is only 18.2%. On the level of 1α,25(OH)2D3 target genes identified by microarray analysis the overlap between the two cellular models is even lower at 5.6%. Although the poor overlap may be partly explained by differences in cutoffs when considering high-confidence peaks and target genes as well as the rather different duration of 1α,25(OH)2D3 treatment (40 min or 4 h versus 36 h), it also suggests that VDR utilizes considerably variant binding sites in different cells. In silico screening for VDR binding showed an approximately 100-fold higher number of putative REs than actually used (50). The main regulator for the access to these sites is the organization of chromatin, which is very cell specific. Moreover, in cells of primary origin there is a significant inter-individual variation in chromatin organization and gene expression. Therefore, it is not surprising that monocytes and lymphoblasts have a rather small overlap of VDR binding sites. This is fully in line with the diversity of VDR target gene sets between different tissues (19–22).

In the absence of ligand, VDR has been shown to actively repress its target genes, likely via a mechanism involving an interaction with co-repressor proteins (17,51). Interestingly, only 14% of the 520 genomic VDR binding locations that occur uniquely in the absence of ligand contain a DR3-type VDRE. In contrast, up to 90% of the best 1α,25(OH)2D3-dependent VDR binding regions contain at least one DR3-type RE. This means that upon ligand treatment the VDR occupancy is strongly shifted from non-DR3 locations to DR3-type RE locations. It is possible that the non-DR3 locations may serve as a nuclear store of VDR to be utilized rapidly upon the introduction of the ligand, partly substituting for the need to transport VDR into the nucleus from outside. This mechanism may also apply to retinoic acid and thyroid hormone receptors, but not for the classical steroid receptors, since they do not bind to DNA in the absence of ligand.

Interestingly, for 300 (73.5%) of the 408 upregulated 1α,25(OH)2D3 target genes we found VDR binding within 400 kb of their TSS, while this applied only for 104 (45.2%) of the 230 downregulated genes. This suggests that the mechanisms of downregulation of VDR target genes are rather complex and/or varied, and may require gene-specific investigations as demonstrated for the CYP27B1 gene (52). In this case the repressive function of VDR is known to result from indirect interaction with the chromatin, via transcription factor 3, also known as VDR interacting repressor (53). An additional question is whether for 1α,25(OH)2D3 target genes there is a ligand-induced derepression mechanism as described for few other members of the nuclear receptor superfamily (54). From the 408 1α,25(OH)2D3 target genes that we can explain via a VDR peak within 400 kb of their TSS only 11 (2.7%) meet the derepression criteria that they have a VDR peak in the unstimulated sample and no peak in the 1α,25(OH)2D3-treated sample. Additional 32 genes (7.9%) can be called dominantly derepressed, since their main peak shows only in the unstimulated sample. This indicates that for some 10% of all 1α,25(OH)2D3 target genes a derepression mechanism may apply.

The fact that DR3-type REs are detected by our de novo motif search as the most dominating sequence in the genomic VDR peak regions confirms earlier reports (9,55) and the general belief in the field that DR3-type REs are the main type of VDREs. Nevertheless, only 31.7% of the identified VDR binding regions contain at least one DR3-type RE, i.e. the majority of genomic VDR binding sites do not use this type of VDRE. Since the most dominant ligand-induced VDR binding sites are highly enriched for DR3-type REs that, moreover, preferentially locate to the neighborhoods of upregulated 1α,25(OH)2D3 target genes, it is likely that DR3-type REs are primarily used in gene activation. However, this leaves the question about the nature of VDR binding at those 68.3% of locations that lack a DR3-type VDRE. Although we could not identify any non-DR3-type VDRE that would explain any significant proportion of the VDR binding, it does not exclude the possibility that a low number of such sites exist. Indeed, a case of repressive gene regulation via an ER9-type VDRE has been demonstrated by Polly et al. (17). Other possible explanations include sites where VDR binds indirectly to genomic DNA via other transcription factors, such as SP1, although again our de novo motif analysis did not identify specific enrichment of other transcription factor binding sites in VDR binding locations lacking a DR3-type VDRE compared to those that had it. A probably more likely explanation for both the VDR-peak-less target genes as well as the target gene-less VDR peaks comes from the recent demonstration that gene regulation by VDR is a very dynamic process with rapid changes of VDR binding site occupancy and, consequently, on the mRNA accumulation of the 1α,25(OH)2D3 target genes (13,56). Therefore, the time points chosen in this study represent only snap-shots of the early actions of 1α,25(OH)2D3 and its receptor VDR. Due to the fast dynamics of VDR binding site occupancy and resulting fluctuations in target gene mRNA amounts, it is likely that without time-course data a considerable proportion of all transient VDR binding sites remains unknown.

Although most of the upregulated genes may be explained by a more unified mechanism of one or multiple VDREs, as already suggested for various VDR target genes, such as CDKN1A (57), CYP24A1 (13,43), ALOX5 (16) and CCNC (14), we demonstrate here with three exemplary genomic regions that there is quite a variation in the gene regulatory scenarios of upregulated genes. Genome-wide there are only about 20 genes that show, as in the case of the SP100 gene, a single VDR location close to the target gene TSS. More common are situations where either one target gene has multiple VDR binding sites in its regulatory region or, as shown for ELL and LRRC25, a pair of closely located VDR target genes share one or more VDR binding sites. The first constellation we have already demonstrated for the above mentioned VDR target genes, while the second scenario applies for the members of the IGFBP gene family (15). However, the most complex regulatory situations we found in the clusters around the VDR target genes THBD and MYO9B, which each contains five or six VDR binding sites of different characteristics.

In conclusion, our genome-wide characterization of early-responding VDR binding sites close to primary 1α,25(OH)2D3 target genes identified a clear shift of VDR binding from non-DR3 sites to DR3-type REs upon ligand treatment. Furthermore, our data suggest regulatory explanations for a large majority of especially upregulated 1α,25(OH)2D3 target genes in THP-1 cells.


The gene expression microarray and sequence data have been submitted collectively to the NCBI Gene Expression Omnibus ( under accession number GSE27438.


Supplementary Data are available at NAR Online.


Academy of Finland; Juselius Foundation, Télévie (Luxembourg); Luxembourgish National Research Fund (FNR); University of Luxembourg. Funding for open access charge: University of Eastern Finland.

Conflict of interest statement. None declared.

Supplementary Material

Supplementary Data:


Kind thanks to the Finnish Microarray and Sequencing Centre in Turku, Finland for microarray services and members of the Benes team for massive parallel sequencing service. S.H., S.V. and C.C. designed the project; S.V. and S.S. performed the experiments; S.H., P.P. and V.B. analyzed the data and S.H. and C.C. wrote the paper.


1. Carlberg C, Polly P. Gene regulation by vitamin D3. Crit. Rev. Eukaryot. Gene Expr. 1998;8:19–42. [PubMed]
2. DeLuca HF. Overview of general physiologic features and functions of vitamin D. Am. J. Clin. Nutr. 2004;80:1689S–1696S. [PubMed]
3. Ingraham BA, Bragdon B, Nohe A. Molecular basis of the potential of vitamin D to prevent cancer. Curr. Med. Res. Opin. 2008;24:139–149. [PubMed]
4. Verstuyf A, Carmeliet G, Bouillon R, Mathieu C. Vitamin D: a pleiotropic hormone. Kidney Int. 2010;78:140–145. [PubMed]
5. Bouillon R, Carmeliet G, Verlinden L, van Etten E, Verstuyf A, Luderer HF, Lieben L, Mathieu C, Demay M. Vitamin D and human health: lessons from vitamin D receptor null mice. Endocr. Rev. 2008;29:726–776. [PubMed]
6. Stoffels K, Overbergh L, Giulietti A, Verlinden L, Bouillon R, Mathieu C. Immune regulation of 25-hydroxyvitamin-D3-1a-hydroxylase in human monocytes. J. Bone Miner Res. 2006;21:37–47. [PubMed]
7. Baeke F, Etten EV, Overbergh L, Mathieu C. Vitamin D3 and the immune system: maintaining the balance in health and disease. Nutr. Res. Rev. 2007;20:106–118. [PubMed]
8. Gombart AF, Borregaard N, Koeffler HP. Human cathelicidin antimicrobial peptide (CAMP) gene is a direct target of the vitamin D receptor and is strongly up-regulated in myeloid cells by 1,25-dihydroxyvitamin D3. FASEB J. 2005;19:1067–1077. [PubMed]
9. Carlberg C, Bendik I, Wyss A, Meier E, Sturzenbecker LJ, Grippo JF, Hunziker W. Two nuclear signalling pathways for vitamin D. Nature. 1993;361:657–660. [PubMed]
10. Quack M, Frank C, Carlberg C. Differential nuclear receptor signalling from DR4-type response elements. J. Cell. Biochem. 2002;86:601–612. [PubMed]
11. Schräder M, Nayeri S, Kahlen JP, Müller KM, Carlberg C. Natural vitamin D3 response elements formed by inverted palindromes: polarity-directed ligand sensitivity of vitamin D3 receptor-retinoid X receptor heterodimer-mediated transactivation. Mol. Cell. Biol. 1995;15:1154–1161. [PMC free article] [PubMed]
12. Tavera-Mendoza L, Wang TT, Lallemant B, Zhang R, Nagai Y, Bourdeau V, Ramirez-Calderon M, Desbarats J, Mader S, White JH. Convergence of vitamin D and retinoic acid signalling at a common hormone response element. EMBO Rep. 2006;7:180–185. [PubMed]
13. Väisänen S, Dunlop TW, Sinkkonen L, Frank C, Carlberg C. Spatio-temporal activation of chromatin on the human CYP24 gene promoter in the presence of 1a,25-dihydroxyvitamin D3. J. Mol. Biol. 2005;350:65–77. [PubMed]
14. Sinkkonen L, Malinen M, Saavalainen K, Väisänen S, Carlberg C. Regulation of the human cyclin C gene via multiple vitamin D3-responsive regions in its promoter. Nucleic Acids Res. 2005;33:2440–2451. [PMC free article] [PubMed]
15. Matilainen M, Malinen M, Saavalainen K, Carlberg C. Regulation of multiple insulin-like growth factor binding protein genes by 1a,25-dihydroxyvitamin D3. Nucleic Acids Res. 2005;33:5521–5532. [PMC free article] [PubMed]
16. Seuter S, Väisänen S, Radmark O, Carlberg C, Steinhilber D. Functional characterization of vitamin D responding regions in the human 5-lipoxygenase gene. Biochim Biophys Acta. 2007;1771:864–872. [PubMed]
17. Polly P, Herdick M, Moehren U, Baniahmad A, Heinzel T, Carlberg C. VDR-Alien: a novel, DNA-selective vitamin D3 receptor-corepressor partnership. FASEB J. 2000;14:1455–1463. [PubMed]
18. Carlberg C, Molnar F. Detailed molecular understanding of agonistic and antagonistic vitamin D receptor ligands. Curr. Top. Med. Chem. 2006;6:1243–1253. [PubMed]
19. Palmer HG, Sanchez-Carbayo M, Ordonez-Moran P, Larriba MJ, Cordon-Cardo C, Munoz A. Genetic signatures of differentiation induced by 1a,25-dihydroxyvitamin D3 in human colon cancer cells. Cancer Res. 2003;63:7799–7806. [PubMed]
20. Wang TT, Tavera-Mendoza LE, Laperriere D, Libby E, MacLeod NB, Nagai Y, Bourdeau V, Konstorum A, Lallemant B, Zhang R, et al. Large-scale in silico and microarray-based identification of direct 1,25-dihydroxyvitamin D3 target genes. Mol. Endocrinol. 2005;19:2685–2695. [PubMed]
21. Suzuki T, Tazoe H, Taguchi K, Koyama Y, Ichikawa H, Hayakawa S, Munakata H, Isemura M. DNA microarray analysis of changes in gene expression induced by 1,25-dihydroxyvitamin D3 in human promyelocytic leukemia HL-60 cells. Biomed. Res. 2006;27:99–109. [PubMed]
22. White SL, Belov L, Barber N, Hodgkin PD, Christopherson RI. Immunophenotypic changes induced on human HL60 leukaemia cells by 1a,25-dihydroxyvitamin D3 and 12-O-tetradecanoyl phorbol-13-acetate. Leuk. Res. 2005;29:1141–1151. [PubMed]
23. Ramagopalan SV, Heger A, Berlanga AJ, Maugeri NJ, Lincoln MR, Burrell A, Handunnetthi L, Handel AE, Disanto G, Orton SM, et al. A ChIP-seq defined genome-wide map of vitamin D receptor binding: associations with disease and evolution. Genome Res. 2010;20:1352–1360. [PubMed]
24. Matilainen JM, Husso T, Toropainen S, Seuter S, Turunen MP, Gynther P, Yla-Herttuala S, Carlberg C, Vaisanen S. Primary effect of 1a,25(OH)2D3 on IL-10 expression in monocytes is short-term downregulation. Biochim. Biophys. Acta. 2010;1803:1276–1286. [PubMed]
25. Gynther P, Toropainen S, Matilainen JM, Seuter S, Carlberg C, Vaisanen S. Mechanism of 1a,25-dihydroxyvitamin D3-dependent repression of interleukin-12B. Biochim. Biophys. Acta. 2011;1813:810–818. [PubMed]
26. van Antwerp DJ, Martin SJ, Kafri T, Green DR, Verma IM. Suppression of TNF-a-induced apoptosis by NF-kB. Science. 1996;274:787–789. [PubMed]
27. Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J, et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004;5:R80. [PMC free article] [PubMed]
28. Du P, Kibbe WA, Lin SM. lumi: a pipeline for processing Illumina microarray. Bioinformatics. 2008;24:1547–1548. [PubMed]
29. Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat. Appl. Genet Mol. Biol. 2004;3:Article3. [PubMed]
30. Klipper-Aurbach Y, Wasserman M, Braunspiegel-Weintrob N, Borstein D, Peleg S, Assa S, Karp M, Benjamini Y, Hochberg Y, Laron Z. Mathematical formulae for the prediction of the residual beta cell function during the first two years of disease in children and adolescents with insulin-dependent diabetes mellitus. Med. Hypotheses. 1995;45:486–490. [PubMed]
31. Keller A, Backes C, Al-Awadhi M, Gerasch A, Kuntzer J, Kohlbacher O, Kaufmann M, Lenhof HP. GeneTrailExpress: a web-based pipeline for the statistical evaluation of microarray experiments. BMC Bioinformatics. 2008;9:552. [PMC free article] [PubMed]
32. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25. [PMC free article] [PubMed]
33. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25:2078–2079. [PMC free article] [PubMed]
34. Robinson JT, Thorvaldsdottir H, Winckler W, Guttman M, Lander ES, Getz G, Mesirov JP. Integrative genomics viewer. Nat. Biotechnol. 2011;29:24–26. [PMC free article] [PubMed]
35. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W, et al. Model-based analysis of ChIP-Seq (MACS) Genome Biol. 2008;9:R137. [PMC free article] [PubMed]
36. Salmon-Divon M, Dvinge H, Tammoja K, Bertone P. PeakAnalyzer: genome-wide annotation of chromatin binding and modification loci. BMC Bioinformatics. 2010;11:415. [PMC free article] [PubMed]
37. Thomas-Chollier M, Sand O, Turatsinze JV, Janky R, Defrance M, Vervisch E, Brohee S, van Helden J. RSAT: regulatory sequence analysis tools. Nucleic Acids Res. 2008;36:W119–W127. [PMC free article] [PubMed]
38. Yordy JS, Moussa O, Pei H, Chaussabel D, Li R, Watson DK. SP100 inhibits ETS1 activity in primary endothelial cells. Oncogene. 2005;24:916–931. [PubMed]
39. Liu L, Ai J, Xiao W, Liu J, Wang Y, Xin D, He Z, Guo Y, Wang Z. ELL is an HIF-1a partner that regulates and responds to hypoxia response in PC3 cells. Prostate. 2010;70:797–805. [PMC free article] [PubMed]
40. Koutsi A, Papapanagiotou A, Papavassiliou AG. Thrombomodulin: from haemostasis to inflammation and tumourigenesis. Int. J. Biochem. Cell. Biol. 2008;40:1669–1673. [PubMed]
41. Greenlee MC, Sullivan SA, Bohlson SS. CD93 and related family members: their role in innate immunity. Curr. Drug Targets. 2008;9:130–138. [PubMed]
42. Fukuda N, Ichihara M, Morinaga T, Kawai K, Hayashi H, Murakumo Y, Matsuo S, Takahashi M. Identification of a novel glial cell line-derived neurotrophic factor-inducible gene required for renal branching morphogenesis. J. Biol. Chem. 2003;278:50386–50392. [PubMed]
43. Meyer MB, Goetsch PD, Pike JW. A downstream intergenic cluster of regulatory enhancers contributes to the induction of CYP24A1 Expression by 1a,25-dihydroxyvitamin D3. J. Biol. Chem. 2010;285:15599–15610. [PMC free article] [PubMed]
44. Kong J, Grando SA, Li YC. Regulation of IL-1 family cytokines IL-1a, IL-1 receptor antagonist, and IL-18 by 1,25-dihydroxyvitamin D3 in primary keratinocytes. J. Immunol. 2006;176:3780–3787. [PubMed]
45. Liu PT, Schenk M, Walker VP, Dempsey PW, Kanchanapoomi M, Wheelwright M, Vazirnia A, Zhang X, Steinmeyer A, Zügel U, et al. Convergence of IL-1b and VDR activation pathways in human TLR2/1-induced antimicrobial responses. PLoS One. 2009;4:e5810. [PMC free article] [PubMed]
46. Kaler P, Augenlicht L, Klampfer L. Macrophage-derived IL-1b stimulates Wnt signaling and growth of colon cancer cells: a crosstalk interrupted by vitamin D3. Oncogene. 2009;28:3892–3902. [PMC free article] [PubMed]
47. Haussler MR, Haussler CA, Jurutka PW, Thompson PD, Hsieh JC, Remus LS, Selznick SH, Whitfield GK. The vitamin D hormone and its nuclear receptor: molecular actions and disease states. J. Endocrinol. 1997;154(Suppl.):S57–S73. [PubMed]
48. Pike JW, Zella LA, Meyer MB, Fretz JA, Kim S. Molecular actions of 1,25-dihydroxyvitamin D3 on genes involved in calcium homeostasis. J. Bone Miner Res. 2007;22(Suppl. 2):V16–V19. [PubMed]
49. White JH. Profiling 1,25-dihydroxyvitamin D3-regulated gene expression by microarray analysis. J. Steroid Biochem. Mol. Biol. 2004;89–90:239–244. [PubMed]
50. Carlberg C, Dunlop TW, Saramäki A, Sinkkonen L, Matilainen M, Väisänen S. Controlling the chromatin organization of vitamin D target genes by multiple vitamin D receptor binding sites. J. Steroid. Biochem. Mol. Biol. 2007;103:338–343. [PubMed]
51. Sanchez-Martinez R, Zambrano A, Castillo AI, Aranda A. Vitamin D-dependent recruitment of corepressors to vitamin D/retinoid X receptor heterodimers. Mol. Cell. Biol. 2008;28:3817–3829. [PMC free article] [PubMed]
52. Turunen MM, Dunlop TW, Carlberg C, Väisänen S. Selective use of multiple vitamin D response elements underlies the 1a,25-dihydroxyvitamin D3-mediated negative regulation of the human CYP27B1 gene. Nucleic Acids Res. 2007;35:2734–2747. [PMC free article] [PubMed]
53. Murayama A, Kim MS, Yanagisawa J, Takeyama K, Kato S. Transrepression by a liganded nuclear receptor via a bHLH activator through co-regulator switching. EMBO J. 2004;23:1598–1608. [PubMed]
54. Baniahmad A, Ha I, Reinberg D, Tsai S, Tsai M-J, O’Malley BW. Interaction of human thyroid hormone receptor b with transcription factor TFIIB may mediate target gene derepression and activation by thyroid hormone. Proc. Natl Acad. Sci. USA. 1993;90:8832–8836. [PubMed]
55. Umesono K, Murakami KK, Thompson CC, Evans RM. Direct repeats as selective response elements for the thyroid hormone, retinoic acid, and vitamin D3 receptors. Cell. 1991;65:1255–1266. [PubMed]
56. Saramäki A, Diermeier S, Kellner R, Laitinen H, Väisänen S, Carlberg C. Cyclical chromatin looping and transcription factor association on the regulatory regions of the p21 (CDKN1A) gene in response to 1a,25-dihydroxyvitamin D3. J. Biol. Chem. 2009;284:8073–8082. [PMC free article] [PubMed]
57. Saramäki A, Banwell CM, Campbell MJ, Carlberg C. Regulation of the human p21(waf1/cip1) gene promoter via multiple binding sites for p53 and the vitamin D3 receptor. Nucleic Acids Res. 2006;34:543–554. [PMC free article] [PubMed]

Articles from Nucleic Acids Research are provided here courtesy of Oxford University Press