Epigenome-wide effects of vitamin D and their impact on the transcriptome of human monocytes involve CTCF

Abstract

The physiological functions of vitamin D are mediated by its metabolite 1α,25-dihydroxyvitamin D₃ (1,25(OH)₂D₃) activating the transcription factor vitamin D receptor (VDR). In THP-1 human monocytes we demonstrated epigenome-wide effects of 1,25(OH)₂D₃ at 8979 loci with significantly modulated chromatin accessibility. Maximal chromatin opening was observed after 24 h, while after 48 h most sites closed again. The chromatin-organizing protein CTCF bound to 14% of the 1,25(OH)₂D₃-sensitive chromatin regions. Interestingly, 1,25(OH)₂D₃ affected the chromatin association of CTCF providing an additional mechanism for the epigenome-wide effects of the VDR ligand. The 1,25(OH)₂D₃-modulated transcriptome of THP-1 cells comprised 1284 genes, 77.5% of which responded only 24 h after stimulation. During the 1,25(OH)₂D₃ stimulation time course the proportion of down-regulated genes increased from 0% to 44.9% and the top-ranking physiological function of the respective genes shifted from anti-microbial response to connective tissue disorders. The integration of epigenomic and transcriptomic data identified 165 physiologically important 1,25(OH)₂D₃ target genes, including HTT and NOD2, whose expression can be predicted primarily from epigenomic data of their genomic loci. Taken together, a large number of 1,25(OH)₂D₃-triggered epigenome-wide events precede and accompany the transcriptional activation of target genes of the nuclear hormone.

INTRODUCTION

Vitamin D has a number of physiological functions, such as maintaining the balance of calcium and phosphorus homeostasis, controlling innate and adaptive immunity and modulating cellular growth, via its biologically most active metabolite 1,25(OH)₂D₃ (1–3). The nuclear hormone is a high affinity ligand of vitamin D receptor (VDR) (4), which is a transcription factor being expressed in most human tissues (http://biogps.org/#goto=genereport&id=7421). Taking all vitamin D responsive organs and cell types together, a few hundred genes are primary targets of VDR and 1,25(OH)₂D₃ (5). However, for a most comprehensive insight on the physiological actions of vitamin D not only short-term effects of 1,25(OH)₂D₃ on primary target genes need to be understood, but also long-term, secondary effects in the time frame of 1–2 days have to be taken into account.

VDR is a member of the nuclear receptor superfamily (6) that in humans comprises 48 members, many of which are modulated in their activity by small lipophilic molecules in the size of cholesterol (7). The principles of nuclear receptor signaling are well understood and seem to be very comparable for most members of the superfamily (8). For example, for VDR the canonical nuclear receptor signaling model is based on DNA binding sites formed by a direct repeat of two hexameric binding motifs spaced by three nucleotides (DR3), on which the receptor binds as a heterodimer with retinoid X receptor (RXR) (9,10). However, the genome-wide binding profile of VDR, as determined by the method chromatin immunoprecipitation sequencing (ChIP-seq) in six human cell culture models (11), indicated that only a minority of 11.5% of the 23 409 individual VDR binding sites carry a DR3-type sequence below their summit. This suggests that VDR has to use additional mechanisms to recognize its genomic targets (12).

After specific binding of their cognate ligand(s) nuclear receptors undergo a conformational change in their ligand-binding domain, which affects the detailed structure of the domain's outer surface and its potential for protein-protein interaction with other nuclear proteins, such as co-activators and co-repressors (13). These co-factors then form larger complexes with chromatin modifying enzymes, such as histone deacetylases (HDACs), histone acetyltransferases and others (14). In most cases the net result of ligand activation of nuclear receptors is the transient opening of chromatin at specific enhancer and transcription start site (TSS) regions resulting in the stimulation of gene transcription (15).

The main components of chromatin are nucleosomes formed of eight histone proteins around whom genomic DNA is wrapped (16). The default state of chromatin is densely packed heterochromatin that protects genes from un-controlled activation (17). This gene silencing is achieved by histone modifications, such as methylation of lysine residues, and the methylation of genomic DNA at cytosines. The positions of these epigenomic modifications are cell-specific and can change in response to environmental changes, such as the activation of signal transduction cascades (18). Accessible chromatin loci can be detected by methods as formaldehyde-assisted isolation of regulatory elements sequencing (FAIRE-seq) (19) and represent only a few percent of the whole genome, such as TSS and enhancer regions. In addition, the 3-dimensional organization of chromatin into larger and smaller scale loops has an important impact on the coordination of gene expression (20). This chromatin looping does not only bring enhancer regions binding transcription factors, such as VDR, into close vicinity of TSS regions, but also sub-divides the genome into functional units, referred to as chromosome domains (21). A key protein in the latter process is the transcription factor CCCTC-binding factor (CTCF) that together with the protein cohesin defines insulator regions separating chromosomal domains from each other (22,23). Experimentally this can be monitored by the method chromatin interaction analysis by paired-end tag sequencing (ChIA-PET) (24).

In this study, we applied FAIRE-seq to 1,25(OH)₂D₃-stimulated THP-1 human monocytes and demonstrated that the chromatin accessibility of nearly 9000 loci distributed along the genome was significantly modulated by the VDR ligand. We characterized these sites for their (i) dynamics in a time frame of 48 h, (ii) overlap with TSS regions, (iii) VDR occupancy and (iv) occurrence of DR3-type binding sites. A screening for non-DR3-type binding motifs below the summits of FAIRE peaks indicated CTCF motifs as top ranking. CTCF ChIP-seq in THP-1 cells did not only confirm this prediction but also demonstrated that genomic CTCF binding can be modulated by 1,25(OH)₂D₃. RNA sequencing (RNA-seq) in THP-1 cells detected 1284 vitamin D target genes within 24 h after ligand stimulation. Our integration of the 1,25(OH)₂D₃-dependent epigenome was able to accurately predict the transcriptional response of 165 of these genes. In summary, we demonstrated that a large number of vitamin D-dependent epigenome-wide events precede and accompany the transcriptional activation of 1,25(OH)₂D₃ target genes.

MATERIALS AND METHODS

Cell culture

The human acute monocytic leukemia cell line THP-1 (25) is a well responding and physiologically meaningful model system for the investigation of 1,25(OH)₂D₃-triggered physiological processes, such as innate immunity and cellular growth (26–29). The 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 were kept at 37°C in a humidified 95% air/5% CO₂ incubator. Prior to mRNA or chromatin extraction, cells were first grown overnight in phenol red-free medium supplemented with charcoal-stripped fetal calf serum and then treated with vehicle (0.1% ethanol (EtOH)) or 100 nM 1,25(OH)₂D₃ (Sigma-Aldrich).

Dicer substrate siRNA (DsiRNA) silencing

THP-1 cells were transfected with either non-specific control (DsiRNA) oligomers (NC1) or specific DsiRNAs targeting CTCF mRNA (Supplementary Table S6). In addition, mock-transfected controls (no DsiRNA) were prepared. The cells were grown for 2 days in RPMI 1640 medium supplemented with 10% FCS, 2 mM L-glutamine, 0.1 mg/ml streptomycin and 100 U/ml penicillin. The transfection was performed in an Amaxa nucleofector I (Lonza) using the Cell Line Nucleofector Kit V according to the manufacturer's instructions with the following modifications. Per sample, 2 x 10⁶ cells were harvested by centrifugation and re-suspended in 100 μl Nucleofector Solution V. A total of 200 pmol of a mixture of three different DsiRNA oligonucleotides was added and the cell suspension was transferred to an electroporation cuvette. The Nucleofector T-12 program was applied. Immediately after the pulse, 500 μl pre-warmed RPMI 1640 medium supplemented with 5% charcoal-stripped FCS, 2 mM L-glutamine, 0.1 mg/ml streptomycin and 100 U/ml penicillin was added and the cell suspension was transferred to pre-equilibrated 6-well plates containing 2.5 ml of medium per well. The extent of the knock-down was determined on mRNA level after an incubation for 48 h including ligand or solvent treatment for 24 h. The transfection efficiency has been determined after nucleofection of a fluorescently labeled DsiRNA oligo (IDT, Leuven, Belgium) by confocal microscopy to be ≥ 88%.

RNA isolation, cDNA synthesis and PCR

Total RNA was extracted using the High Pure RNA Isolation Kit (Roche) according to the manufacturer's instructions. RNA quality was assessed by native agarose gel electrophoresis. cDNA synthesis and qPCR were performed as described previously (30). qPCR reactions were performed with the LightCycler^® 480 System (Roche) using 250 nM of reverse and forward primers, 2 μl cDNA and the LightCycler 480 SYBRGreen I Master mix (Roche). Primer-specific temperatures and sequences are listed in Supplementary Table S7. Relative mRNA expression levels were determined using the formula 2^−(ΔCt), where ΔCt is Ct_{(target gene)} – Ct_{(reference gene)}. The genes B2M, GAPDH and HPRT1 were used as references as described previously (28).

CTCF ChIP

ChIP assays were performed as described by Zhang et al. (31) with some modifications. After treatment of 20 x 10⁶ THP-1 cells, nuclear proteins were cross-linked to genomic DNA by adding formaldehyde directly to the medium to a final concentration of 1% and incubating at room temperature for 10 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 10 min on a rocking platform. The cells were collected by centrifugation and washed twice with ice cold phosphate-buffered saline (PBS). The cell pellets were subsequently resuspended twice in 10 ml cell lysis buffer (0.1% SDS, 1 mM EDTA, 150 mM NaCl, 1% Triton X-100, 0.1% sodium deoxycholate, protease inhibitors, 50 mM HEPES-KOH, pH 7.5) and once in 10 ml nuclear lysis buffer (1% SDS, 1 mM EDTA, 150 mM NaCl, 1% Triton X-100, 0.1% sodium deoxycholate, protease inhibitors, 50 mM HEPES-KOH, pH 7.5). After two washes with cell lysis buffer, the chromatin pellet was resuspended in 700 μl of SDS lysis buffer (1% SDS, 10 mM EDTA, protease inhibitors, 50 mM Tris-HCl, pH 8.1) and the lysates were sonicated in a Bioruptor Plus (Diagenode) to result in DNA fragments of 200 to 500 bp. Cellular debris was removed by centrifugation. 340 μl aliquots of the lysate were diluted 1:5 in IP dilution buffer (1% Triton X-100, 2 mM EDTA, 150 mM NaCl, protease inhibitors, 20 mM Tris-HCl, pH 7.5). 4 μl anti-CTCF antibody (Millipore, 07–729) was coated to 60 μl Dynabeads Protein G (Invitrogen) in an overnight incubation at 4°C. The pre-formed bead-antibody complexes were then washed twice with beads wash buffer (0.1% Triton X-100, PBS, protease inhibitors) and added to the chromatin aliquots. The samples were incubated for overnight at 4°C on a rotating wheel to form and collect the immuno-complexes. The beads were washed sequentially for 5 min on a rotating wheel with 1 ml of the following buffers, each: twice cell lysis buffer, once high salt buffer (0.1% SDS, 1% Triton X-100, 1 mM EDTA, 350 mM NaCl, 0.1% sodium deoxycholate, 50 mM HEPES-KOH, pH 7.5), once ChIP wash buffer (250 mM LiCl, 1% Nonidet P-40, 0.5% sodium deoxycholate, 1 mM EDTA, 10 mM Tris-HCl, pH 8.0) and twice TE buffer (1 mM EDTA, 10 mM Tris-HCl, pH 8.0). Then, the immune complexes were eluted using 250 μl ChIP elution buffer (1% SDS, 10 mM EDTA, 50 mM Tris-HCl, pH 7.5) at 37°C for 30 min with rotation. The elution was repeated with a 10 min rotation and the supernatants were combined. The immune complexes were reverse cross-linked at 50°C for 2 h in the presence of proteinase K (Fermentas) in a final concentration of 160 μg/ml. The genomic DNA was isolated with the ChIP DNA Clean&Concentrator Kit (Zymo Research).

FAIRE

FAIRE analysis was performed according to the protocol published by Giresi et al. (32) with some modifications. In short, 18 x 10⁶ THP-1 cells were cross-linked identically as for ChIP. After 10 min cross-linking and stopping with glycine the washed cell pellets were resuspended and incubated sequentially in 2 ml of buffer L1, 2 ml of buffer L2 and 700 μl of buffer L3. The lysates were sonicated in a Bioruptor to result in DNA fragments of 200 to 500 bp and cellular debris was removed by centrifugation. Input samples were reverse cross-linked overnight at 65°C. The FAIRE samples and reverse cross-linked reference samples were subjected to two sequential phenol/chloroform/isoamyl alcohol (25/24/1) extractions, resuspended in 10 mM Tris-HCl (pH 7.4) and treated with 1 μl of RNase A (10 mg/ml) for 1 h at 37°C. The genomic DNA was purified using the ChIP DNA Clean&Concentrator Kit.

ChIP-seq, FAIRE-seq and RNA-seq analyses

ChIP and FAIRE DNA templates and total RNA samples were sequenced at 50 bp read length using standard manufacturer protocols at the Gene Core at the EMBL (Heidelberg, Germany). In parallel, public ENCODE ChIP-seq data sets (33) from K562 cells for CTCF (GSM749690) and the transcription factors ETS1 (GSM803442), GABPA (GSM803524), NR2F2 (GSM1010782) and NR4A1 (GSM777637) were downloaded. All ChIP-seq data were (re)analyzed at harmonized settings: alignment with the human reference genome version hg19 using Bowtie software version 1.1.1 (34) with the following essential command line arguments: bowtie -n 1 -m 1 -k 1 -e 70 –best. The aligned input and CTCF reads were converted to sorted BAM format using samtools (35) and, after merging the read sets per sample, converted to TDF format using igvtools, in order to allow efficient visualization in the Integrative Genomics Viewer (IGV) genome browser (36). Statistically significant ChIP-seq peaks were identified using Model-based Analysis of ChIP-Seq data (MACS) version 2 (37) with the following essential command line arguments: macs2 callpeak –bw 150 –keep-dup 1 –g hs –qvalue = 0.01 –m 5 50 –bdg. Otherwise, default parameters were used.

For FAIRE-seq data analysis statistically significant peaks were identified using the Zinba program package version 2.02 by setting the mean fragment length at 200 bp and using other settings as recommended for FAIRE-seq in the Zinba website (http://code.google.com/p/zinba/wiki/UsingZINBA) including peak refinement (38). To enable comparisons across the full time series for all positions with a significant FAIRE peak in at least one of the time points, common peak borders were identified by detecting the local maxima of original peak border density and filtering out all resulting intervals that were > 2 kb. The peak statistics (except the false discovery rate (FDR)) were recalculated as done in MACS version 1.3.7.1 and all intervals where the maximal fold enrichment over the entire data set was < 5 were filtered out.

RNA-seq reads were mapped to the reference genome (hg19) using Bowtie (http://bowtie-bio.sourceforge.net/index.shtml). Then, raw counts and fragments per virtual kb million were obtained for profiling gene expression. Finally, differential expression was calculated with Cuffdiff (http://cole-trapnell-lab.github.io/cufflinks/cuffdiff) based on a t-test over the ratio of normalized counts for the treated and untreated conditions using a cutoff FDR of 0.5%. CTCF ChIP-seq, FAIRE-seq and RNA-seq raw data are available at Gene Expression Omnibus (GEO; www.ncbi.nlm.nih.gov/geo) at GSE69303.

Functional characterization of FAIRE peak loci

More than 80 000 FAIRE-seq peaks were identified at one or more time points. In order to obtain higher accuracy in data integration (see below), we were interested in persistent peaks and filtered out some 18 000 peaks (22% of all) that were not present at all four time points (0, 2, 24 and 48 h). The summit regions (±100 bp) of these 62 231 peaks were specifically searched for DR3-type sequences and other transcription factor binding motifs by using HOMER (39) with score 7. For a better characterization of the dynamics of the 8979 sites of 1,25(OH)₂D₃-modulated chromatin accessibility along the four time points, we applied self-organizing maps (SOM). SOM (40) is an unsupervised method that allows the identification of representative profiles of the inspected data. In the context of this contribution, SOM received as input vectors the peak intensity along the four time points, or the fold change (FC), of the aforementioned 8979 peaks. Each component of the vector represented a specific time point, and thus the input space is four-dimensional. SOM classifies these peaks in one of nine possible profiles in a 3 × 3 lattice, which allows a visual inspection of the similarities of the chromatin accessibility as a function of time.

Gene ontology analysis

Gene ontology analysis was conducted through the Ingenuity Pathway Analysis (IPA) software (www.qiagen.com/ingenuity). The genes with a significantly differential expression as calculated by Cuffdiff as well as their logFC were provided as input to the software. Three different data sets were constructed, one for each time point (2.5, 4 and 48 h).

Statistics

From 62 231 FAIRE-seq peaks we computed a significant modulation (ligand-effect) comparing the intensity of the peaks at 2, 24 and 48 h to the co-located peak intensity at time point 0 h. Three replicate experiments were conducted for each of the four time points allowing to compute statistical significance. At each time point, a peak was declared to have a ligand-effect, if the P-value of a paired Student's t-test applied to the peak intensity of the three replicates at that time point and the three replicates at 0 h was lower than 0.05. We focused on the group of 8979 peaks that presented a significant ligand-effect at one or more time points. Since the distribution of the FC does not fit to a Gaussian distribution (Supplementary Figure S2), a non-parametric analysis of variance by the Kruskal–Wallis test was calculated (H(3) = 780.61, 0.0063). This test offers the basis to reject the null hypothesis that all three data are coming from the same distribution, and at least one sample dominates at least one of the other two. Next, we applied a pairwise Mann–Whitney–Wilcoxon test and found also that the FC was significant in all cases (U(2 h, 24 h) = 87.3, 0.00038; U(2 h, 48 h) = 94.2, 0.00092; U(24 h, 48 h) = 102.8, 0.00024).

RESULTS

Genome-wide changes in chromatin accessibility in response to 1,25(OH)₂D₃

In three independent experiments THP-1 human monocytes were stimulated for 0, 2, 24 and 48 h with 1,25(OH)₂D₃ and chromatin was isolated for FAIRE-seq analysis. This approach identified 62 231 genomic regions that showed at all four time points significantly (FDR < 0.1%) accessible chromatin. Importantly, 8979 (14.4%) of these FAIRE peaks were at least at one time point significantly (P < 0.05) modulated by VDR ligand treatment (Supplementary Figure S1A and Table S1). In reference to untreated cells (0 h), 7716 (85.9%) regions were ligand-responsive only at one, 1197 (13.3%) at two and 66 (0.7%) at three time points. Already after 2 h stimulation with 1,25(OH)₂D₃ 3323 genomic regions were affected, of which 2754 (82.8%) showed increased chromatin accessibility, and after 24 h even 4262 (92.9%) out of 4586 ligand-responsive chromatin regions opened up (Supplementary Figure S1B). In contrast, after 48 h only 689 (28.6%) out of 2399 modulated FAIRE peaks were induced by VDR ligand, i.e. at this late time point the majority of the 1,25(OH)₂D₃-responsive chromatin loci were less accessible than at their basal state. Four representative example regions illustrate the increase of chromatin accessibility (Supplementary Figure S1C) or chromatin closing (Supplementary Figure S1D) in response to VDR ligand treatment.

The global effect of 1,25(OH)₂D₃ treatment on the number of open chromatin regions and their grade of ligand inducibility at different time points was illustrated by histogram plots (Supplementary Figure S2). The average increase in chromatin accessibility of all significantly ligand-responsive FAIRE peaks after 2 h ligand treatment was already 1.57-fold and increased even to 1.76-fold after 24 h (Figure 1A). In contrast, after 48 h the average induction of the ligand-responsive regions was only 0.89-fold, i.e. they were 1.11-fold repressed compared to the basal state. This suggested that the average chromatin accessibility of the genome in human monocytes (i) increased significantly (P < 0.001) after 2 h stimulation with 1,25(OH)₂D₃, (ii) was maximal at 24 h and (iii) after 48 h the chromatin was even slightly less accessible than before onset of stimulation. A very comparable result was obtained when all 22 autosomal chromosomes and chromosome X were analyzed individually (Figure 1B), i.e. on each chromosome there was always more 1,25(OH)₂D₃-mediated chromatin opening after 24 h than after 2 h, while after 48 h chromatin was closed again. A more detailed analysis using 5000 equally sized (621 kb) windows over the whole human genome (3.1 Gb, Figure 1C), indicated 390 regions that (i) contain one to six ligand-responsive FAIRE peaks and (ii) in average were more than 2-fold up- or down-regulated (Supplementary Table S1). Only two (0.5%) of these regions were observed at all three time points and 19 (4.9%) at two time points, i.e. the vast majority (369 (94.6%)) of them were observed only at one time point. In addition, we identified 95 regions that contain six or more ligand-modulated FAIRE peaks (Supplementary Table S1).

In summary, in THP-1 cells 8979 chromatin regions, i.e. approximately every seventh accessible site, were significantly modulated by 1,25(OH)₂D₃. At earlier time points (2 and 24 h) most regions opened up, but later (48 h) 1,25(OH)₂D₃ reduced the accessibility of the majority of them.

Properties of 1,25(OH)₂D₃-modulated chromatin regions

From 62 231 regions of the human genome that in THP-1 cells were located within accessible chromatin, 12 515 (20.1%) overlapped with TSS regions (Supplementary Figure S3A). Also in the ligand-responsive subsets of 3323 (2 h), 4586 (24 h) and 2399 (48 h) FAIRE peaks the TSS overlap rate was 18.9%, 19.8% and 20.0%, respectively, being very similar in all time points. Next, we performed SOM analysis of the time course response of the absolute chromatin accessibility of all 1,25(OH)₂D₃-responsive genomic regions (Supplementary Figure S3B). This approach subdivided the 8979 time course profiles into nine groups (I to IX) containing 69 to 1779 genomic sites, each. The groups III, VI and VII showed a more than 2.5-times higher overlap with TSS regions (50.7–55.8%) than the average (20%). In total, these three groups contained 723 genomic regions that co-locate with 400 TSS regions displaying the overall highest levels of open chromatin. In addition, these regions were characterized by a rapid increase of chromatin accessibility within the first 2 h after ligand treatment. In contrast, groups V and VIII showed with 7.3% and 7.6% a more than 2.5-times lower TSS overlap rate than average, i.e. the 3062 members of these groups may primarily represent enhancer regions. Furthermore, rather low absolute levels of open chromatin and minor ligand-dependent increases of chromatin accessibility characterized these regions. The remaining four groups (I, II, IV and IX) represented 5220 genomic regions that showed average TSS rates (16.4–29.6%) and maximal chromatin accessibility 24 h after onset of stimulation with 1,25(OH)₂D₃.

Traditional characteristics of 1,25(OH)₂D₃ signaling are the presence of VDR and the involvement of DR3-type binding sites at the concerned genomic loci (12). Based on VDR ChIP-seq data from THP-1 cells (11,29) only 824 (1.3%) of the 62 231 sites of open chromatin were occupied by VDR (Figure 2A). In parallel, HOMER screening identified DR3-type motifs only within 3626 (5.8%) accessible chromatin regions. Of these 3626 DR3 sites only 170 overlapped with VDR binding loci. This means that at only 20.6% of all 824 VDR-occupied sites the receptor had access to a DR3-type sequence for the formation of heterodimers with RXR. In the subsets of 1,25(OH)₂D₃-responsive chromatin regions at 2, 24 and 48 h the VDR occupancy rate increased to 1.7%, 2.6% and 1.9%, respectively, while the DR3 rate raised slightly to 6.7%, 7.0% and 7.2% (Figure 2B). In addition, the DR3 overlap rate of the VDR binding sites at 2, 24 and 48 h ligand treatment increased to 25.5%, 33.9% and 34.8% (in average 29.4%), respectively. However, this means that the majority of all accessible DR3-type sites (87.1–93.7%) were not occupied by VDR and that the receptor used in some 70% of its binding events an alternative mode of complex formation than VDR-RXR heterodimers on DR3-type sequences (Figure 2C).

In order to identify the most interesting subsets of 1,25(OH)₂D₃-responsive chromatin sites, we repeated the SOM analysis with the 8979 FC-derived profiles of the time course. This again resulted in nine groups (A to I) with 60 to 2565 members (Figure 2D). Groups A and B showed with 18.3% and 10.2% not only a clearly increased VDR rate but with 30.0% and 11.1%, respectively, their DR3 rate was far higher than average. The 168 members of both groups were characterized by an in average more than 4-fold increase in chromatin accessibility displaying a maximum at 24 h after ligand treatment. The 2204 members of groups E and H showed a similar behavior, but displayed with VDR rates of 3.3% and 3.4% and DR3 rates of 9.4% and 8.6%, respectively, less pronounced enrichments of these key characteristics of 1,25(OH)₂D₃ signaling. The 119 members of group D showed a prominent induction of chromatin accessibility within the first 24 h, but they did not overlap with any VDR binding site and displayed with 6.7% a rather mean DR3 rate. The 1292 members of group G also showed with 0.3% a very low VDR rate and with 4.6% a slightly reduced DR3 rate, but they represented the only subset of the time course series that displayed continuous chromatin closing at all time points after ligand treatment. Finally, groups C, F and I showed with VDR rates of 1.5%, 1.5% and 1.1% and DR3 rates of 6.6%, 5.3% and 7.0%, respectively, rather average characteristics. The 5196 members of these groups demonstrated less prominent opening of chromatin at time points 2 and 24 h and in contrast to all other time courses their chromatin accessibility returned after 48 h back to the basal state. However, in accordance with the results presented in Figure 1, all time course subsets showed lower chromatin accessibility after 48 h ligand treatment than after 24 h. Chromatin regions representing groups A and B, respectively, are illustrated in Figure 2E (left example: group A, right example: group B).

Taking together, SOM analysis allowed the distinction of different profile subsets within the group of 8979 sites displaying 1,25(OH)₂D₃-modulated chromatin accessibility. Regions overlapping with TSS regions showed high levels of chromatin accessibility and were quickly inducible, while enhancer-type regions were less prominent and responded delayed. Furthermore, on regions with a high VDR occupancy rate chromatin was opening prominently, whereas low VDR rates mostly indicated ligand-dependent chromatin closing.

1,25(OH)₂D₃-inducible CTCF binding

Only 608 (6.8%) of the 8979 ligand-modulated chromatin sites contained a DR3-type motif and even only 170 of all sites (1.9%) were occupied by VDR (see Figure 2B). In fact, only 50 (29.4%) of these 170 VDR binding sites contained a DR3-type sequence. This suggested that (i) the very most of the effects of 1,25(OH)₂D₃ on chromatin opening or closing could not be explained by VDR binding to the same site and (ii) even in cases where VDR was occupying such a ligand-dependent chromatin locus there must be further mechanistic explanations for the effects of 1,25(OH)₂D₃ in addition to VDR-RXR heterodimers binding to DR3-type sequences. Therefore, we used HOMER to screen the sequences below the summits (±100 bp) of all 8979 ligand-modulated FAIRE peaks for transcription factor binding motifs. The 10 most significantly enriched motifs were binding sites for the insulator proteins CTCF and BORIS, the pioneer factor PU.1 and several members of the ETS transcription factor family (Supplementary Figure S4A). This suggested that these transcription factors could occupy 1,25(OH)₂D₃-modulated, accessible chromatin sites and may contribute to the mechanisms of action of the VDR ligand on the epigenome.

We focused on the top ranking transcription factor CTCF and performed ChIP-seq analysis with chromatin from THP-1 cells that were stimulated for 24 h with 1,25(OH)₂D₃ in comparison to vehicle-treated cells. In total, we found 71 993 CTCF peaks, but from the 46 691 and 46 922 CTCF binding events in presence and absence of 1,25(OH)₂D₃, respectively, only 21 690 (≈46%) were common (Figure 3A and Supplementary Table S2). 52.2% of all 71 993 CTCF peaks in THP-1 cells were found at the same position as in the ENCODE tier 1 cell line K562 (33), for which our harmonized re-analysis of the CTCF ChIP-seq data indicated 98 515 peaks (Supplementary Figure S4B). In the absence of ligand 65.9% of the 46 920 CTCF sites overlapped with K562 data and for the 21 690 common CTCF sites the rate was even 91.7%.

From the 8979 ligand-dependent FAIRE peaks 1249 (13.9%) overlapped with these common CTCF sites (Supplementary Table S1). These 1249 CTCF binding sites co-located with 234 (18.7%) TSS regions and 21 (1.7%) VDR sites, i.e. these were average percentages (compare Supplementary Figure S3 and Figure 2) not indicating any enrichment for these types of sites. Interestingly, 5915 common CTCF sites were 2-fold or more induced or reduced after 24 h treatment with 1,25(OH)₂D₃ (Figure 3B and Supplementary Table S2). Moreover, 1752 common CTCF sites overlapped with ligand-dependent FAIRE peaks, 597 of which were also sensitive to 1,25(OH)₂D₃ (Figure 3B and Supplementary Table S2). Four representative examples illustrate the co-location of ligand-modulated FAIRE peaks with sites of either up- or down-regulated CTCF association after treatment of the cells with 1,25(OH)₂D₃ (Figure 3C).

From the 5915 ligand-dependent CTCF sites, 609 (10.3%) overlapped with a TSS region, which was a 2.3-fold increased percentage compared to all common CTCF sites (4.4%, Supplementary Table S2). Similarly, the percentage of overlap with 1,25(OH)₂D₃-modulated FAIRE peaks was increased 2.5-fold from 4.0% for all 71 993 CTCF sites to 8.8% for the 21 690 common CTCF sites. In contrast, comparing both groups, there was no significant increase for the rate of DR3-type sequence (7.4 versus 8.6%). Moreover, a HOMER analysis for enriched motifs below the CTCF peaks did not provide any evidence for an enrichment of VDR binding events. In parallel, there were only 37 (0.6%) ligand-dependent CTCF sites overlapping with a VDR binding site (Supplementary Table S2). Compared to the transcription factors ETS1, GABPA and NR2F2 (also called COUP-TF2), which in the ENCODE tier 1 cell line K562 (33) showed overlap rates with CTCF of 10.1%, 7.7% and 15.1%, respectively (Supplementary Figure S4C), VDR showed a very low rate. Also from the perspective of the rather low total number of 1142 VDR binding sites in THP-1 cells 37 VDR peaks overlapping with CTCF represent only 3.2% of all, which was even 4-fold lower than respective rates for the nuclear receptor NR4A1 (also called NUR77). This indicates that VDR is not directly involved in the 1,25(OH)₂D₃-dependent association or dissociation of CTCF with its binding sites, i.e. most likely VDR and CTCF are not part of the same nuclear protein complex.

The functional impact of CTCF for the regulation of 1,25(OH)₂D₃ target genes by the ligand was tested by silencing CTCF mRNA expression (Figure 3D and Supplementary Figure S5). After the transfection of THP-1 cells with DsiRNA oligonucleotides against CTCF, the cells were treated with 1,25(OH)₂D₃ for 24 h. The knock-down of CTCF mRNA expression was 70% and 67% for solvent- and ligand-treated cells, respectively, compared to cells transfected with a negative control DsiRNA (NC1). The additionally performed mock-transfected controls differed in average by 9.6% from the NC1 control (standard deviation 9.6%, data not shown). From in total 20 tested 1,25(OH)₂D₃ target genes, six genes showed significant modulations of either their relative mRNA expression or their ligand-inducibility due to the CTCF silencing. We observed both increased (CAMP and ALOX5) and decreased (NOD2, ZFP36, THBD and MYC) mRNA levels and the changes ranged from –34% (NOD2) to +37% (CAMP). Likewise, the induction of the 1,25(OH)₂D₃ target genes by the ligand was both positively and negatively modulated by the CTCF silencing. The CAMP gene showed with 33% the strongest attenuation of its ligand response, while the induction of the ALOX5 gene increased by even 69%. In contrast, the ZFP36 mRNA was reduced to a similar extent for both solvent- and ligand-treated cells (22% and 24%, respectively), thus its inducibility by 1,25(OH)₂D₃ did not change. Interestingly, at the ligand-dependent CTCF sites (Figure 3B) closest to the TSS regions of the six genes, which showed sensitivity to CTCF knock-down, CTCF association was more prominently induced or reduced by 1,25(OH)₂D₃ than at those of the remaining 14 genes that were not affected by CTCF silencing. An exception to this observation is the CTCF site 86 kb downstream of the TMEM37 gene TSS, which is nearly 5-fold induced (Supplementary Table S2). However, it is possible that it lies outside of the chromatin domain of the TMEM37 gene and therefore does not influence the gene's expression.

In summary, ligand-modulated chromatin sites were not only highly significantly enriched for CTCF binding motifs, but 13.9% of these 8979 loci also bound the transcription factor. Surprisingly, the association of CTCF with some of these genomic sites could be up- or down-regulated by 1,25(OH)₂D₃, which provided an additional mechanistic explanation for the epigenomic effects of the VDR ligand. VDR is not directly involved in the ligand-modulated binding of CTCF, though. The actions of CTCF affected 30% of all tested 1,25(OH)₂D₃ target genes.

Transcriptome-wide effects of 1,25(OH)₂D₃ stimulation

Next, we performed triplicate RNA-seq experiments with THP-1 cells that were for 2.5, 4 and 24 h stimulated with 1,25(OH)₂D₃ (Supplementary Table S3). From 23 622 human genes 9220 (39.0%) were not expressed, but from the 14 402 expressed genes 1246 (8.7%) were at least at one time point significantly (P < 0.0001) modulated by the VDR ligand (Supplementary Figure S6A). Of note, 38 genes, such as the well-known VDR target gene CYP24A1 (41), were at basal level below the expression threshold (0.3) but passed it after stimulation with 1,25(OH)₂D₃ (gray in Supplementary Table S3). Therefore, in the following we refer to in total 1284 1,25(OH)₂D₃ target genes. Forty-six of these genes were already responding 2.5 h after onset of 1,25(OH)₂D₃ stimulation, 288 after 4 h and 1204 after 24 h (Figure 4A). The expression of only 45 genes was affected by the VDR ligand at all three time points (Supplementary Table S3). Interestingly, after 2.5 h stimulation all 46 differentially expressed genes were up-regulated, after 4 h 170 (59.0%) out of 288 and after 24 h 663 (55.1%) out of 1204 genes (Figure 4B, Supplementary Figure S6B and S6C). This suggests that, like in the epigenome analysis (Figure 1), at later time points the 1,25(OH)₂D₃-dependently down-regulated part of the transcriptome had a far larger impact. A gene ontology analysis using IPA software of the 1,25(OH)₂D₃-regulated genes at each of the three time points indicated that the top-ranking biofunction represented by the respective genes shifts from ‘anti-microbial response’ at 2.5 h to ‘endocrine system disorder’ at 4 h and ‘connective tissue disorders’ at 24 h (Figure 4C and Supplementary Table S4). The second ranking term changed from ‘gene expression’ to ‘renal and urological disease’ and ‘hepatic system disease’.

Taken together, the 1,25(OH)₂D₃-modulated transcriptome of THP-1 cells comprised 1284 genes, of which 681 (53.0%) were up-regulated and 603 (47.0%) down-regulated. The vast majority of the genes (995 (77.5%)) responded with a delay of 24 h to stimulation with 1,25(OH)₂D₃. Interestingly, during the time course of 1,25(OH)₂D₃ treatment the proportion of down-regulated genes increased from 0 to 44.9% and in parallel the main biofunction of the genes shifted from anti-microbial response to connective tissue disorder.

Integration of epigenome- and transcriptome-wide responses to 1,25(OH)₂D₃

After testing several models for epigenome and transcriptome data integration, such as Pearson product-moment between pairs of variables presented in Supplementary Tables S1 and S3, that delivered only low statistical correlations, we followed a machine learning perspective (42). This approach provided the clearest, although not univariate correlation between epigenomic and transcriptomic data for FAIRE peaks overlapping with the TSS regions of 165 genes that were significantly regulated after 1,25(OH)₂D₃ stimulation for 24 h. These 165 genomic loci of accessible chromatin were a subset of 1790 FAIRE peaks that overlapped with at least one TSS region. The logFC of these 165 genes was discretized into one of six states using the k-means algorithm (k = 6, Supplementary Figure S7A). That is, based on their responsiveness to 1,25(OH)₂D₃ stimulation the genes were subdivided into six different categories (Supplementary Table S5) and further characterized. For that, we applied a feature selection algorithm based on random forests, in order to rank the relevance of the epigenomic (Supplementary Table S1) and transcriptomic (Supplementary Table S3) parameters in the description of all 8979 FAIRE peaks (Supplementary Figure S7B) and the prediction of the expression state of the TSS regions of the 165 selected genes (Supplementary Figure S7C). When considering all 8979 peaks, transcriptomic data demonstrated higher impact than the rest of the attributes, while for the description of the 165 selected genes the epigenomic data were more important. Nevertheless, the combined discriminative power of epigenomic and transcriptomic data was higher than the discriminative power of each category alone, with a clearer subsuming effect for the 165 selected genes (Supplementary Figure S7D and S7E). Importantly, for 138 out of 165 genes their inducibility after 24 h 1,25(OH)₂D₃ treatment could be accurately (FDR < 5%) predicted by the six most relevant parameters listed in Supplementary Figure S7C (Supplementary Figure S7F).

Two representative examples of the 165 selected 1,25(OH)₂D₃ target genes were presented in Figure 5. CTCF ChIA-PET data from K562 cells indicated that the HTT gene is co-located in the same chromatin domain with the genes GRK4 and HTT-AS (Figure 5A). However, the HTT gene was more prominently expressed than the two other genes and the only gene that after 4 and 24 h was significantly induced by 1,25(OH)₂D₃ (Figure 5B). The HTT gene is widely expressed and is required for normal development, but in cases when it acquires additional unstable trinucleotide repeats it can cause the neurodegenerative disorder Huntington's disease, which is characterized by loss of striatal neurons (43). At two loci, 450 and 4895 bp downstream of the HTT TSS, VDR binding sites overlapped with ligand-inducible FAIRE sites (Figure 5A). In contrast, the chromatin domain of the NOD2 gene did not contain any VDR binding site, but the chromatin accessibility at its TSS region was induced after treatment with 1,25(OH)₂D₃ (Figure 5C). The NOD2 protein acts as an intracellular pattern recognition receptor binding bacterial peptidoglycans and stimulates the innate immune response (44). The NOD2 gene is co-located with the genes NDK1, SNX20 and CYLD, but only NOD2 responded significantly to 1,25(OH)₂D₃ after 24 h stimulation (Figure 5D).

In summary, the integration of epigenomic and transcriptomic data identified 165 1,25(OH)₂D₃ target genes, the expression of which could be predicted primarily from epigenomic data of their genomic loci. For example, the 1,25(OH)₂D₃ target gene HTT is regulated directly via VDR binding close to its TSS region, while 1,25(OH)₂D₃-depending chromatin opening at the NOD2 TSS contributes to the regulation of the gene.

DISCUSSION

This study describes the epigenome-wide effects of 1,25(OH)₂D₃ on chromatin accessibility of human monocytes. During a time course of 48 h the natural VDR ligand significantly affected the opening and closing of chromatin at nearly 9000 sites, which was one in seven of all accessible genomic regions in this cellular model. Interestingly, both in number of affected loci as well as in average strength of the chromatin opening, the 24 h stimulation with 1,25(OH)₂D₃ was more prominent than the treatment for 2 h or 48 h. Epigenomic effects are assumed to be mediated primarily by chromatin modifying enzymes that act quickly (45). In fact, there was a subset of 1,25(OH)₂D₃-modulated chromatin sites, in particular at TSS regions, that already showed maximal accessibility 2 h after onset of stimulation. However, the majority of the sites reached their maximal levels of chromatin opening delayed after 24 h. This suggests that the process of 1,25(OH)₂D₃-mediated chromatin opening involves multiple steps and that the effects observed 24 h after onset of stimulation are not primary responses but rather secondary effects on chromatin accessibility. The observation that most 1,25(OH)₂D₃-sensitive chromatin sites return after 48 h back to basal levels indicates that the epigenome-wide effects of the stimulation with VDR ligand are transient. This suggests that for a persistent effect a constant vitamin D supplementation is required.

Since VDR is the only nuclear protein that binds 1,25(OH)₂D₃ with high affinity (46), the nuclear receptor must be central to the mechanisms initiating the epigenome-wide effects of its natural ligand. However, in accordance with previous genome-wide studies on vitamin D signaling (11,29) only a very low number of 1,25(OH)₂D₃-sensitive chromatin sites (1.9%) contain a DR3-type binding site that is essential for the well-understood binding of VDR-RXR heterodimers. Even the small subset of 168 rapidly responding genomic regions with high VDR occupancy have a DR3 rate of only 17.9%, which is in the same order as obtained in the VDR ChIP-seq data set meta-analysis (11). This suggests that the canonical model of VDR-RXR heterodimers binding to DR3-type sequences applies only for a subset of the primary responding chromatin sites and already in the majority of the fast responding sites VDR seems to be involved in different types of complexes. In some of these complexes VDR may be partnered with other nuclear proteins than RXR, in order to contact DNA. As suggested already a longer time ago (47,48), these could be other members of the nuclear receptor superfamily. However, motif screenings indicated that primarily rather ubiquitously binding transcription factors of the ETS family co-locate with genomic VDR binding sites (11). It has not yet been shown that VDR forms any heterodimeric complexes with these proteins, but they may act as pioneer factors supporting VDR in opening chromatin regions. In parallel, it is possible that at a number of sites VDR does not associate directly to genomic DNA but may bind ‘piggyback’ on DNA-associated transcription factors (12).

In previous studies we already described that several members of the HDAC family are primary 1,25(OH)₂D₃ target genes (49) and that HDAC proteins and their small molecule inhibitors are involved in the control of VDR target genes (28,50). In this study, we demonstrated that CTCF co-locates with some 14% of all 1,25(OH)₂D₃-sensitive chromatin sites. This provides an additional mechanistic explanation for the epigenome-wide actions of the VDR ligand. However, we could not find any evidence that VDR is directly involved in the association of CTCF with DNA, i.e. VDR and CTCF appear not to be part of the same nuclear protein complex. Interestingly, silencing of CTCF affected some 30% of 1,25(OH)₂D₃ target genes by either further activating or inhibiting them. The TSS regions of these CTCF-sensitive 1,25(OH)₂D₃ target genes are rather close to strongly ligand-sensitive CTCF binding sites. A down-regulation of CTCF expression likely leads to changes in chromatin architecture, so that the contact of TSS regions of 1,25(OH)₂D₃ target genes with VDR binding enhancer regions will be rendered possible or impossible. In MCF-7 cells, where CTCF mediates the association of repressive chromatin regions with the KCNQ5 gene locus, CTCF silencing leads to disruption of this hub and in the following to the down-regulation of KCNQ5 expression (51). Likewise, CTCF-mediated chromatin loops may confine VDR binding sites in context of at least some of its target genes. However, local changes in chromatin organization after knock-down of CTCF could also facilitate interactions of TSS regions with VDR binding loci or even allow contacts to additional VDR sites, which may explain the observed increase in ligand-inducibility of some genes. As CTCF is involved in the maintenance of H3K27me3 levels on the inactive X chromosome (52), CTCF silencing would lead to reduced chromatin compaction, which could provide the mechanism for our observation that some genes are up-regulated after CTCF silencing.

The novel finding that the associations of CTCF with its genomic targets were modulated by 1,25(OH)₂D₃ suggests that the CTCF-dependent organization of chromatin can be modulated by vitamin D. This could affect the size and function of chromosomal domains, i.e. such a mechanism has wider consequences than a 1,25(OH)₂D₃-dependent action at any specific chromatin site. The 1,249 1,25(OH)₂D₃-sensitive chromatin loci that bind CTCF showed only low percentages of VDR occupancy and TSS overlap. CTCF tends to associate preferentially with loci where chromatin accessibility was modulated early after ligand treatment and in particular then the association of CTCF with genomic DNA increased. This suggests that CTCF may be involved in the early epigenomic effects of 1,25(OH)₂D₃. We assume that CTCF affects also the action of other transcription factors, including other nuclear receptors. In fact, CTCF seems to co-locate with a number of these regulatory proteins in the same nuclear complex, i.e. CTCF may use additional mechanisms for the interference with the respective signaling pathways.

Within the first 24 h of stimulation with 1,25(OH)₂D₃ at some 90% of the 7193 genomic loci, which were accessible after 2 h ligand treatment, chromatin accessibility increased while in parallel the RNA level of 681 genes was up-regulated and of 603 genes down-regulated. The more than 5-times higher number of 1,25(OH)₂D₃-sensitive chromatin loci compared to that of differentially expressed genes suggests that a gene may be regulated via several chromatin sites. In parallel, for a number of chromatin loci their modulation by 1,25(OH)₂D₃ may be without any functional consequence. Interestingly, the comparison of these numbers implies that also for the down-regulation of gene transcription, chromatin sites rather needed to be open than closed. In accordance with this, the very early primary transcriptome-wide response to 1,25(OH)₂D₃ stimulation was exclusively up-regulation. This suggests that for molecular mechanisms of down-regulation of gene transcription more time and probably more steps are necessary than in gene activation (53). The latter is a straightforward mechanism that follows for some major sites the canonical model of 1,25(OH)₂D₃ signaling (12), while for down-regulation many different mechanism seem to apply (54).

It can be assumed that the vast majority of the 995 genes whose mRNA level is regulated only after 24 h of stimulation with 1,25(OH)₂D₃, such as NOD2, are secondary targets of VDR. This means that no VDR binding site needs to be found within their regulatory region but instead effects of 1,25(OH)₂D₃ on the accessibility of the respective chromatin regions are detected. This emphasizes the impact of epigenome-wide data collection, such as chromatin accessibility as well as binding of VDR, CTCF and other nuclear proteins, over pure RNA quantifications, in order to understand the mechanisms of the regulation of these 1,25(OH)₂D₃ target genes. Moreover, even in cases, where the regulatory mechanism is rather obvious, such as for the two VDR binding sites close to the TSS of the HTT gene, the 1,25(OH)₂D₃-dependent modulation of chromatin accessibility at these sites provides a confirmation of their impact for the control of the respective genomic region and the genes located at these sites. Nevertheless, the fact that the integration of epigenomic and transcriptomic data can so far only correctly predict the action of 165 of 1284 1,25(OH)₂D₃ target genes, indicates that we are still missing insight into the complexity of vitamin D signaling, i.e. additional data need to be assessed and integrated into the model. The Roadmap Epigenomics Consortium (55) has recently presented the integration of 111 human epigenomes. They based their findings on an impressive set of epigenome-wide data, including different histone modifications and DNA methylation, from a large variety of tissues, but could only provide a model for the basal state of the human epigenome. This means that for the different challenges of the human body, such as a supplementation with vitamin D, additional data need to be collected.

The shift of the main physiological function ‘anti-microbial response’ represented by the primary vitamin D targets to ‘connective tissue disorders’ by secondary targets characterizes and reflects the functional profile of VDR and its ligand 1,25(OH)₂D₃. The primary function of ancestral, probably orphan nuclear receptors was the control of metabolism (56) and nowadays still many members of the nuclear receptor superfamily have this property. During evolution the superfamily diversified and the different members not only acquired high affinity ligand binding but also additional roles in physiology. As a side effect of evolutionary and developmental pressures, VDR became a key controller of innate immunity but also to affect cellular growth in many of its target tissues (57).

In conclusion, in this study we demonstrated the dynamic nature of the epigenetic landscape of human monocytes in response to stimulation with 1,25(OH)₂D₃. The characterization of the epigenomic changes within the chromosomal domain of 1,25(OH)₂D₃ target genes contributes to the understanding of the involved mechanisms. The variety of observed effects of 1,25(OH)₂D₃ is an indication for a larger number of different mechanisms comprising many nuclear proteins, such as the ones demonstrated here for CTCF. The results of this study present mechanisms and factors that significantly extend the canonical model of 1,25(OH)₂D₃ signaling.

This work was supported by the Academy of Finland (#267067) and the Juselius Foundation. Kind thanks to the Gene Core Facility at the EMBL in Heidelberg, Germany, for massively parallel sequencing services.

FUNDING

Funding for open access charge: Academy of Finland and Juselius Foundation.

Conflict of interest statement. None declared.

REFERENCES