Genome-wide analysis of LXRα activation reveals new transcriptional networks in human atherosclerotic foam cells.

Increased physiological levels of oxysterols are major risk factors for developing atherosclerosis and cardiovascular disease. Lipid-loaded macrophages, termed foam cells, are important during the early development of atherosclerotic plaques. To pursue the hypothesis that ligand-based modulation of the nuclear receptor LXRα is crucial for cell homeostasis during atherosclerotic processes, we analysed genome-wide the action of LXRα in foam cells and macrophages. By integrating chromatin immunoprecipitation-sequencing (ChIP-seq) and gene expression profile analyses, we generated a highly stringent set of 186 LXRα target genes. Treatment with the nanomolar-binding ligand T0901317 and subsequent auto-regulatory LXRα activation resulted in sequence-dependent sharpening of the genome-binding patterns of LXRα. LXRα-binding loci that correlated with differential gene expression revealed 32 novel target genes with potential beneficial effects, which in part explained the implications of disease-associated genetic variation data. These observations identified highly integrated LXRα ligand-dependent transcriptional networks, including the APOE/C1/C4/C2-gene cluster, which contribute to the reversal of cholesterol efflux and the dampening of inflammation processes in foam cells to prevent atherogenesis.

INTRODUCTION

Cardiovascular and metabolic disorders have become an enormous burden for human health. Homeostasis of cholesterol at the blood vessel wall is important to cope with elevated cholesterol levels that contribute to increased rates of atherosclerosis and cardiovascular disease (1). A key event in the development of atherosclerosis consists of the uncontrolled uptake of oxidized low-density lipoproteins (oxLDL) by macrophages recruited at the subendothelial space of vessel walls (2). In concert with local inflammatory reactions, this process leads to the formation of lipid-loaded macrophages termed foam cells (3). The liver x (nuclear) receptors, LXRα and LXRβ, are ligand-dependent factors that regulate cholesterol homeostasis and reverse cholesterol transport; hence, they are interesting drug targets for treating cardiovascular diseases (4). Oxysterols are natural ligands and activators of LXRα (5), which is the predominant LXR subtype required for the full anti-atherogenic action of LXR agonists in inhibiting atherosclerosis (6). Evolutionary analysis of the sequences of both LXR subtypes suggested that a single LXR gene duplicated during mammalian evolution (7) and resulted in one general and one more specialized factor. Consistent with this hypothesis, LXRβ was found to be ubiquitously expressed, whereas LXRα is more restricted to cell types that modulate cholesterol and fatty acid metabolism, most importantly in liver and macrophages (8). Hypothesis-driven research revealed important insights, including LXRα regulation at selected genomic loci of target genes, such as the ATP-binding cassette transporters ABCA1 and ABCG1, or apolipoprotein E, to understand mechanistic aspects of foam cell development and atherosclerosis (9).

Here, we applied genome-wide analysis to shed light on the contribution of ligand-based LXRα regulation to these processes. We performed integrative genomic analyses with the aim to decipher LXRα-dependent functional features and transcriptional regulation pathways, including the effects of anti-atherosclerotic drug treatment. Key networks determined by LXRα activation led to biological function such as counteracting lipid-overload in macrophages, which could in principle not be fully derived from single gene events and analysis thereof.

To analyze the role of LXRα modulation in atherosclerosis and related diseases we used a well-validated human macrophage and foam cell model and applied the efficient synthetic LXR agonist T0901317.

MATERIALS AND METHODS

Cell models

Human monocytic leukaemia THP1 cells were maintained in RPMI 1640 medium (Sigma-Aldrich) supplemented with 10% fetal bovine serum (FBS, Biochrom) and differentiated for 48 h using 10e-8 M phorbol 12-myristate 13-acetate (PMA, Sigma-Aldrich). Human primary macrophages were isolated from at least four individual buffy coats donated by healthy volunteers (kindly provided by DRK-Blutspendedienst Ost Gemeinnützige GmbH Institut Berlin). Peripheral blood monocytes (PBM) were isolated from buffy coats with Ficoll Paque and MACS Monocyte Isolation Kit II with MACS LS columns (Miltenyi Biotec). Monocytes were differentiated for 7 days. Primary and THP1 macrophages were treated similarly either with 0.01% dimethylsulfoxide (DMSO) or 1 µM T0901317 (Sigma-Aldrich) for 24 h. Foam cell formation was induced by 100 µg/ml oxLDL (Autogen Bioclear UK Ltd) for 48 h and subsequently treated with 0.01% DMSO or 1 µM T0901317 for 48 h. Cholesterol loading and treatment was controlled by Oil Red O staining; detailed cholesterol composition was assessed with the fluorometric method of Amplex Red Cholesterol Assay Kit (Invitrogen).

LXR knockdown

For siRNA-mediated LXRα/β-knockdown, differentiated THP1 cells were transfected with TransIT-TKO transfection reagent (Mirus) and 15 nM LXRα Silencer Validated siRNA (ID 5458) and 15 nM LXRβ Silencer Select Validated siRNA (ID s14684) or 30 nM Silencer Select Negative Control #1 (all from Ambion). Transfection was carried out for 48 h followed by treatment of macrophages and foam cells with 10 µM T0901317 or 0.1% DMSO for 24 h.

Western blotting

Western blotting was performed with commercial antibodies against LXRα [Abcam, ab 41902, (10,11)], LXRβ [Abcam, ab56237, (12)], RXRα (Santa Cruz Biotechnology, sc-774 X) and β-actin (Santa Cruz Biotechnology, sc-47778, C4). Secondary antibodies were horseradish peroxidase (HRP)-labelled anti-mouse and anti-rabbit (Santa Cruz Biotechnology). Densitometry was performed using ImageQuant TL (GE Healthcare).

Chromatin immunoprecipitation and sequencing

Chromatin immunoprecipitation (ChIP) was performed with Diagenode’s Transcription Factor ChIP Kit (Diagenode). For immunoprecipitation, we used the well-validated antibodies against LXRα (Abcam, ab 41902), LXRβ (Abcam, ab56237), H3K4me3 (pAB-003-050, Diagenode), anti-H4K20me1 (ab9051, Abcam) and negative control IgG (kch-819-015, Diagenode). For one LXRα ChIP-seq reaction, we pooled at least four individual ChIP reactions. Sequencing of LXRα was performed twice (individual biological replicates) for macrophages and T0901317 treated macrophages. Sequence reads of 36 bp were obtained using the second-generation Genome Analyzer and the Solexa Analysis Pipeline (Illumina). Reads were mapped to human genome assembly (February 2009, GRCh37/hg19) using Bowtie (13). Sequencing data were submitted to EBI and can be accessed via www.ebi.ac.uk/ena/data/view/ERP001502. Peak calling and LXRα binding profiles were generated genome-wide using the model based-analysis of ChIP-seq algorithm [MACS, (14)]. Validation of sequencing results was performed with ChIP–quantitative polymerase chain reaction (qPCR) (SYBR Green PCR Master Mix, Applied Biosystems). The relative occupancy of the immunoprecipitated factor at a locus was estimated using the following equation: 2 × e (Ct input − Ct ChIP). LXRα ChIP qPCR was normalized against IgG ChIP qPCR of the same locus. Detailed description of ChIP procedure, peak calling and filtering can be found in the Supplementary Methods and Supplementary Datasets S1.

Peak analysis

To compare different ChIP-seq lanes and to detect differential LXRα binding across investigated cell models, we performed LXRα peak-enrichment normalization (15) for all cell models. Differential binding events were based on a change of 1.5-fold. For visualization, raw tag data was loaded to seqMINER (16). Genomic distributions were determined using the Cis-regulatory Element Annotation System (17).

Annotation of genes controlled by nearby peaks

As we assumed that one LXRα-binding site can regulate multiple surrounding target genes (18,19), we decided to annotate all potential LXRα target genes with a maximal distance of 200 kb from peak centre (Supplementary Datasets S2). Therefore, we used the Peak Center Annotation script (peak2gene) from the Cistrome Analysis Pipeline (cistrome.org). Gene definitions were taken from the UCSC Genome Browser’s RefGene table (20).

Formaldehyde-assisted isolation of regulatory elements followed by sequencing

Formaldehyde-assisted isolation of regulatory elements (FAIRE) was performed as described previously (21). For relative openness analysis, enrichment intensities were extracted, quantile normalized and averaged for each genomic region. Detailed descriptions can be found in Supplementary Methods.

LXRα motif analyses

Motif search was performed de novo with top 100 bound sequences for each defined peak set individually using the MEME-ChIP tool from MEME suite [http://meme.sdsc.edu, (22)] with default settings. Derived motifs with P-value of < 10e-7 were scanned using available position weight matrices from the Transfac database (23). Motif distribution in all LXRα peak set sequences was determined by the find motif tool FIMO (MEME suite tool) with threshold P-value of 10e-4.

Reporter gene assays

Reporter gene assays were performed with natural and mutated LXR response elements (LXREs). Therefore, five copies of the LXREs were cloned into pGL4.31 vector (Promega) with the In-Fusion HD EcoDry cloning system (Clontech Takara Bio Europe). Full-length LXRα and RXRα were cloned from cDNA fragments (Source BioScience clone IRATp970C0271, Gene ID: 7376 and clone IOH39435, Gene ID: 6256) into pBIND vector (Promega). Detailed descriptions can be found in supplementary experimental procedures. For reporter gene analysis, HEK293T cells were co-transfected with the pGL4.31–LXRE–Luc, pBIND–LXRα and pBIND–RXRα in 0.25% Lipofectamine (Invitrogen) for 4 h and subsequently treated with 10 µM T0901317 or 0.1% DMSO for 24 h. Reporter activity was determined with Dual Luciferase Reporter System (Promega). Activity was measured in relative luciferase units and finally visualized as fold-change of T0901317 treatment from vehicle control.

Gene expression

RNA extraction, microarray and quantitative PCR analyses were performed as described elsewhere (24). Differential expression analysis was performed on background subtracted data with cubic spline normalization and Benjamini Hochberg false discovery rate (FDR) correction. Significant data were considered to have a detection P-value of <0.01 and differential P-value of <0.05. Gene expression data were submitted in MIAME-compliant form to the ArrayExpress database under accession number E-MTAB-1106 (www.ebi.ac.uk/arrayexpress). Detailed descriptions and gene lists can be found in Supplementary Methods and Supplementary Datasets S3.

Correlation analyses

For correlation, we compared differentially expressed genes (versus vehicle treated macrophages) with binding site associated genes. For visualization, we chose to plot the normalized LXRα peak enrichment and differential gene expression sorted in 4 up- and 4 down- quantiles according to the expression fold-change. For gene activity prediction, we determined promoter-specific changes of H3K4me3, H4K20me1 and chromatin accessibility at genes with a nearby binding (25). Further, we generated a relative histone and FAIRE-seq signal of our cell models versus untreated macrophage and built the mean signal ±1.5 kb of the target gene transcription start site (TSS) and plotted the mean signal and differentially expressed genes in quintiles for each cell model (26). For prediction, we chose the mean signal fold-change of expression quantile 4. To determine a significant impact of LXRa binding on gene expression, we correlated LXRa peak locations and differential expression of nearby genes and compared this with a randomly shuffled peak set. Initially, we investigated this relationship as described by Boeva et al. (27). The significance of this spatial relationship was tested by three statistical hypothesis tests, hypergeometric test, Kolmogorov-Smirnov test (K-S test) and empirical cumulative distribution function (ECDF)-based test (Anderson-Darling test) (28).

Gene ontology analysis

Gene ontology (GO)-enrichment analysis was performed using the Database for Annotation, Visualization and Integrated Discovery [DAVID, http://david.abcc.ncifcrf.gov/,(29)] for the model-specific and knockdown-validated LXRα target genes with a binding site. We performed functional annotation clustering for GO term Bioprocess_FAT with highest classification stringency.

Pathway analysis

Pathway analysis was also performed with DAVID analysing the Kyoto Encyclopedia of Genes and Genomes (KEGG) and BIOCARTA pathways with the functional annotation chart option. Information on regulation of lipid metabolism by peroxisome proliferators-activated receptor alpha (PPARα) was extracted from Reactome pathway database [www.reactome.org, (30)]. Pathway analysis for T0901317-specific LXRα target genes in foam cells was performed with Ingenuity Pathways Analysis [IPA, www.Ingenuity.com, (31)]. We performed core analysis using the Ingenuity Knowledge Base as reference set and considered direct and indirect relationships with high or experimentally observed confidence.

Association of LXRα-binding sites with genome-wide association studies

Correlation of LXRα-binding data with genome-wide association studies (GWAS) was processed by overlapping the NHGRI GWAS catalogue (32) single-nucleotide polymorphism (SNPs) positions with our defined LXRα-binding sites. We chose a P-value threshold of <10e-5. Additionally, we controlled that the LXRα-binding sites of interest were within the linkage disequilibrium (LD) block [DistiLD Database, (33)] of the SNP. Further, we considered just LXRα peak-associated genes that were also reported in the GWAS. Detailed list can be found in Supplementary Datasets S4.

Functional network analysis

Interaction networks were derived from the FANTOM4-EdgeExpress Database (34) and STRING (35). Networks were visualized using Cytoscape (36). Differentiation between already known and new LXRα target genes and interactors were done with BIOGRID, NEXTBIO, Nuclear Receptor Resource [(37,38); http://nrresource.org] databases and most recent publications on LXRα (19,39).

Statistical analysis

The Student’s t-test was used to calculate the statistical significance. When multiple samples were analyzed the ANOVA test was applied followed by Dunnett’s multiple comparison test with P < 0.05 considered as statistically significant. All results represent the mean ± standard deviation.

RESULTS

LXRα binding is highly ligand dependent in THP1 and human primary macrophage models

Knowledge about the binding of LXRα to the human genome is required to lay the foundation for deciphering specific sets of regulatory regions of this ligand-dependent nuclear receptor. LXRα is by nature a flexible sensor for diverse metabolites in the human body, which can respond quickly to changing ligand concentrations. To generate stringent data sets for further analyses, we considered only LXRα-binding sites, which resulted in differential gene expression and were LXR knockdown sensitive (Supplementary Figure S1A).

First, we examined the amount of LXRα and LXRβ proteins and their heterodimerization partner RXRα. Because of ligand-based activation in human macrophages, LXRα is increasingly expressed in an auto-regulatory fashion, whereas comparably low proteins amounts of the LXRβ subtype and RXRα showed no significant change (Figure 1A, Figure 1B and Supplementary Figure S1B). Notably, mouse macrophages do not show a feed-forward loop of LXRα expression (40), indicating that activation of LXRα in human cells differs strikingly from mouse foam cells. Cholesterol loading and beneficial effects of T0901317 were further confirmed by Oil red O staining (Supplementary Figure S1C) and cholesterol composition analysis (Figure 1C). To determine LXRα binding in THP1-derived macrophages and foam cells in presence and absence of synthetic LXRα ligand T0901317, we applied ChIP using a well-validated and previously applied antibody (10,41) followed by massively parallel deep sequencing. Macrophages and T0901317-treated macrophages were sequenced in biological duplicates and reached correlation values of r = 0.98 and r = 0.92, respectively (Supplementary Figure S1D). In ligand-free macrophages with low amounts of LXRα (Figure 1B), its genome-wide enrichment at potential binding sites was mostly below the defined threshold for detection of significantly enriched LXRα-binding sites (Supplementary Figure S1E). To exclude any bias during ChIP-seq data processing, we selected LXRα-binding sites and validated successfully 21 ChIP-seq peaks by ChIP–qPCR analysis, including confirmation of the extremely low abundance of LXRα in the absence of ligand (Figure 1D and Supplementary Figure S1F and G). Characteristically, LXRα binding was highly induced by its synthetic ligand T0901317 in contrast to ligand-independent, constitutive binding of the LXRβ subtype. As LXRα was 38 times more abundant in T0901317-treated foam cells than its β-subtype, we consistently observed up to 80-fold enrichment of LXRα, whereas we found only 2.5-fold enrichment for LXRβ, as detected at the ABCA1 LXR response element (LXRE) locus (Figure 1D and Supplementary Table S1).

Figure 1.

In human macrophage and foam cell models, LXRα is regulated in an auto-regulatory fashion. (A) Sketch of LXRα auto-regulatory function. LXRα, in complex with RXRα and a specific ligand, is bound to a direct repeat of 4 (DR4) response elements in front of LXRα gene locus. Ligand-dependent activation of LXRα gene expression is regulated via a feed-forward loop mechanism. (B) Western blot analysis of LXRα, LXRβ, RXRα and ACTB (β-actin) in macrophages and foam cells in the presence and absence of T0901317. (C) Cholesterol composition analysis in macrophages, foam cells and T0901317-treated foam cells expressed as micromolar cholesterol per gram protein. Total cholesterol (black), free cholesterol (grey) and esterified cholesterol (white). The average from six biological replicates with standard errors of mean is shown. Cross indicates statistically significant, P < 0.05. (D) ChIP–qPCR validation of LXRα and LXRβ enrichment at ABCA1 LXRE locus in THP1-derived macrophage models and LXRα in human PBM-derived models relative to IgG control.

To validate the THP1 model, we used human PBM-derived macrophages and produced foam cells thereof, and obtained similar LXRα-binding profiles (Figure 1D). Interestingly, we detected spurious ligand-free LXRα-binding background only in PBM-derived macrophages, indicating experimentally interfering lipids observed in the donor blood. THP1 cells were required to gain sufficient and standardized material for further analyses and were, therefore, the experimental resource of choice.

Three sets of LXRα genomic loci define gene regulation in macrophages and foam cells

To understand oxLDL- and/or T0901317-induced cell-specific and shared binding events of LXRα, we applied comparative ChIP-Seq analyses (15) and could categorize three different sets of LXRα genome binding (Figure 2A). Fifty-five per cent of LXRα peaks (1459 peaks) were different among cell models. Thereof, 29% (769 peaks) showed prominent enrichment in T0901317 treated macrophages and foam cells (T0901317 specific set), and we observed an enrichment of 26% (690 peaks) of LXRα sites in foam cells over T0901317-treated cells, which we subsumed as foam cell-specific set. Forty-five per cent of LXRα peaks (1193 peaks) were shared among all cell models, which we termed shared peak set. Open chromatin is associated with active transcriptional regulators and pre-determines transcription factor binding (42,43). Using FAIRE-seq analysis, we observed pronounced chromatin openness of LXRα loci of shared binding sites and T0901317-specific binding sites compared with foam cell-specific–binding sites, similarly in all cell models (Figure 2A and Supplementary Figure S2A). Further, we observed that open LXRα-binding sites were enriched with other co-localized transcription factor-binding sites (TFBS, data extracted from Encode) compared with oxysterol-induced foam cell-specific–binding sites (Supplementary Figure S2A, TFBS). Although T0901317 treatment apparently induced LXRα and co-activator assembly at highly accessible transcriptional hotspots, oxLDL treatment seemed to trigger less pronounced LXRα binding to DNA. As shown in Figure 2B and Supplementary Figure S2B and C, characteristic loci for each set of the LXRα-binding site included the well-known LXRα target genes ABCA1, LXRα itself and the APOE/C2/C4-gene cluster, or unexpected new target genes as poly(ADP-ribose)polymerase-1, PARP1 or pre-B-cell leukemia homeobox 4, PBX4.

Figure 2.

Differential LXRα-binding sites in human macrophage and foam cell models. (A) Comparative LXRα ChIP-seq enrichment heatmap and mean tag density ±1.5 kb from peak center. Peaks were sorted according to their degree of variability, calculated as log2 fold-change between models. The following sets were defined: i) T0901317-specific set (29%); ii) all models that share similar LXRα enrichment are referred to as shared peaks (45%); iii) foam cell-specific peak set (26%). Additionally to the LXRα ChIP-seq, we also plotted the average chromatin openness of the three peak sets derived from FAIRE-seq analyses. Macrophages are represented in grey, T0901317-treated macrophages in brown, foam cells in green, T0901317-treated foam cells in blue and average FAIRE in black. (B) Representative tag alignment tracks of LXRα ChIP-seq, which were normalized against IgG control lane tags of same loci. Set-specific examples of genomic regions with differential or shared LXRα binding in our models. APOC1 locus (T0901317 specific set), LXRα locus (shared set) and PARP1 locus (foam cell specific set). Arrows indicate transcription start sites and orientation of transcription. Black arrows under tracks show peaks.

LXRα peak enrichment at transcriptional start sites of target genes is sharpened by T0901317 treatment

To further characterize shared and differential LXRα-binding sites, we performed de novo motif search for each of the three sets of specific LXRα genomic sites. Only in the shared and T0901317-specific LXRα peak sets we found significant enrichment of the direct repeat spacer 4 (DR4) motif, which is known from targeted gene analyses [Figure 3A, (44)]. In contrast, the foam cell-specific set of LXRα-binding site sequences did not significantly enrich any motif. Similar was true for the targeted scanning of de novo extracted LXR:RXR motif from the shared peak set. With a fairly stringent P-value cut-off of 10e-4, we found that 15% (398 peaks) of all LXRα-binding sites occupied this DR4 motif, but only 2% thereof were derived from foam cell-specific set (Supplementary Figure S3A). To biochemically characterize the derived de novo consensus LXR:RXR motif, we performed transient reporter gene assays (Supplementary Figure S3B). Surprisingly, single-nucleotide mutations of LXR response elements did hardly change reporter gene activity. Major mutations in LXRα halfsites were required to significantly decrease reporter activity. These observations were consistent with other nuclear receptor studies, suggesting that LXRα could bind to fairly degenerated DR elements (19,39).

Figure 3.

Differential binding but shared gene expression properties in T0901317-treated and -untreated foam cells. (A) De novo motif analyses. Motif enrichment was analysed with the top 100 ChIP-seq sequences from the three peak sets followed by TRANSFAC motif search. (B) Genomic positions of LXRα peaks. Fold-change of LXRα peak distribution was compared with genome composition for each differential peak set. The following features are shown from left to right: promoter ≤1 kb, promoter 1–2 kb, promoter 2–3 kb, downstream ≤1 kb, downstream 1–2 kb, downstream 2–3 kb, 5′–untranslated region (UTR), 3′ UTR, coding exon, distal intergenic and intron. (C) Distribution of all LXRα-binding sites over the human genome and association with gene expression. Peaks that are located in an area with no genes ±200 kb (grey), peaks annotated to near genes (distance <200 kb) without indication of differential expression (compared with ligand-free macrophage) (blue). All peaks annotated to near genes (distance <200 kb) with indication of differential expression (white) including peaks near differentially expressed genes (orange), peaks near LXR knockdown-sensitive genes (8%, stringent set) and peaks near genes that showed a relevant chromatin modification change at the transcriptional start site (light blue). (D) Correlation of differential expression fold-change (separated in 4 up- and 4 down- quantiles) with normalised LXRα peak enrichment for T0901317-treated macrophages (brown), T0901317-treated foam cells (blue) and foam cells (green).

We also tested the observed higher degree of conservation in the spacer region at positions 8 and 9. Mutations in the spacer region resulted in a strong decrease of reporter gene expression, indicating potential stereochemical disturbance of LXRα binding.

To further analyse the spatial characteristics of LXRα binding, we determined the genomic positions of LXRα-binding sites. This analysis revealed that T0901317-induced LXRα binding increased enrichments around transcription start sites and promoters, whereas oxysterol induced foam cell specific LXRα binding sites were enriched downstream of genes (Figure 3B). This analysis suggested a more pronounced shaping of transcriptional initiation by the nanomolar binding ligand T0901317 in contrast to only micromolar binding by natural sterols.

To study the effects of genome-wide LXRα binding on transcriptional networks, we performed gene expression and LXR knockdown analyses of the investigated cellular models. By correlating the differentially expressed target genes (P ≤ 0.05 versus macrophage) with the genome-wide LXRα binding map, 19% of the LXRα-binding sites overlapped with differentially expressed genes, thereof 8% were associated with 186 genes that were significantly sensitive to LXR knockdown (Figure 3C, stringent set). The observed enrichment of differentially expressed genes close to LXRα-binding sites was statistically significant (P < 0.005 for K–S test and hypergeometric test and P < 0.01 for ECDF test) and was absent in a randomly shuffled set of binding sites (Supplementary Figure S3C).

In addition to gene expression data, we generated for all cell models genome-wide maps from histone modifications H3K4me3 and H4K20me1, which are important transcriptional initiation or elongation marks, respectively (45). Combined with FAIRE-seq data of open chromatin, we determined promoter-specific changes of potential LXRα target genes that harbour enriched LXRα-binding sites. These changes correlated well with gene expression profiles [Supplementary Figure S3D and E, (25)] and could be used to identify 492 additional LXRα target genes. This approach further increased the total correlation of LXRα-binding sites with target genes to 30% (Figure 3C). Correlation analysis between gene expression and LXRα peak enrichment at differentially regulated LXRα target genes revealed a bimodal regulation of gene expression in all cell models (Figure 3D). Expression levels separated in quantiles showed increased LXRα-binding enrichment in quantiles −2 to −4 and 2 to 3. Remarkably, expression quantile 4 in all cell models was correlated with decreased LXRα peak enrichment, suggesting minor effects of LXRα on this gene set.

Despite similar expression profiles, we observed less enriched LXRα binding at differentially regulated target genes in untreated foam cells compared with T0901317-treated foam cells. This was also consistent with mentioned observations of less efficient positioning of LXRα binding to promoter sites and DR4 sequence motifs. Taken together, we observed a bimodal expression profile in all described cell models with less pronounced binding properties in foam cells, which could be substantially sharpened by T0901317 treatment.

Main functions of LXRα in cholesterol metabolism, apoptosis and interaction with the PPARα signalling pathway

For the functional annotation, we first analysed the stringently validated set of 186 knockdown-sensitive LXRα target genes and deciphered specific pathways enriched in foam cells, T0901317-foam cells and T0901317-treated macrophages (Supplementary Table S2). In vehicle-treated foam cells, we discovered “organic acid biosynthetic process” and “regulation of cholesterol storage”. Further, we found 14 genes associated with the regulation of apoptosis, indicating increasing cell death that could contribute to the formation of atherosclerotic plaques in vivo. Foam cells treated with the synthetic LXRα ligand T0901317 did not significantly enrich this disease-associated pathway. Instead they showed enrichment of the BIOCARTA pathway “nuclear receptors in lipid metabolism and toxicity”, and the bioprocess “negative regulation of cholesterol storage”, which is one of the main functions of this anti-atherogenic compound. Interestingly, the PPARα signalling pathway was enriched in T0901317-treated and -untreated foam cells, indicating tight interactions of LXRα- and PPARα-regulated pathways, as reported recently for macrophages (39). Further investigation of the involvement of LXRα in PPARα signalling pathway in T0901317-treated foam cells revealed a major proportion of shared target genes (Supplementary Figure S4). Strikingly, the LXRα complex also regulated the expression of two major co-activators of PPARα, Cbp/p300 interacting transactivator with Glu/Asp-rich carboxy terminal domain 2 (CITED2) and peroxisomal proliferator-activated receptor A interacting complex 285 (PRIC 285). For set-specific functional annotation, we found similar results as for the stringent set of LXRα target genes (Supplementary Table S3), with the addition of defence response enrichment in shared peak set and the involvement of Ras protein signal transduction in foam cell-specific peaks set.

LXRα binding and genetic variants associated with complex diseases reveal key target genes in transcriptional networks

Numerous genome-wide association studies (GWAS) were successful in associating genetic variants and genomic loci to common diseases (46). But in many cases, genetic variation associated to disease phenotypes could not functionally explain the underlying mechanisms. To discover the potential impact of genetically influenced LXRα-induced transcriptional networks on complex diseases, we searched for significant, disease relevant SNPs. We considered SNPs in linkage disequilibrium of cell model-specific LXRα-dependent–binding sites that showed differentially expressed target genes. As listed in Table 1, we found a number of LXRα-dependent genes with central impact on metabolism and inflammation. The new LXRα target gene homeobox protein PBX4 (Supplementary Figure S2C) showed striking association with LDL cholesterol metabolism, whereas another LXRα target gene, the LDL-associated phospholipase A2 (PLA2G7) was associated with “lipoprotein-associated lipase activity and mass”. The most striking connection between LXRα binding and disease related loci was found with SNP rs4420638, which was associated to “C-reactive protein” with a P-value of 9e-139 and “LDL cholesterol” with a P-value of 1e-60 in meta-analyses (Table 1). The SNP rs4420638 is located close to the transcription termination site of APOC1 and is encircled by two LXRα-binding sites that are in 5.7- and 6.5-kb distance (Figure 4A). LXRα-binding peaks detected at this locus were T0901317 specific (Figures 2B and 4A), whereas the natural ligand oxLDL was not efficient in recruiting LXRα to this locus. Consistently, T0901317 treatment led to an improved expression profile of the whole-gene cluster (Supplementary Figure S5). In this context, it is interesting to note that APOC1 gene is not only located within a cluster of genes that are implicated in plasma lipid metabolism [APOE, APOC4 and APOC2, (47)] but also surrounded by genes (BCL3, RELB, PVRL2) that are involved in biological processes, such as inflammation, immunity, differentiation, cell growth, tumorigenesis and apoptosis (first two via NF-κB signalling). To provide a frame on the global impact of APOC1 and surrounding genes on lipid metabolism and inflammatory response, we performed transcriptional network analysis of this LXRα locus. To embed the gene cluster around APOC1 into a broader metabolic context, we further considered well-known LXRα target genes and generated a network with differential expression data of the foam cell using the STRING database (Figure 4B). This network illustrated the tight relationship between direct and indirect LXRα target genes in foam cells and their involvement in major processes, such as lipid metabolism, inflammation and apoptosis. Thus, we can conclude that LXRα-dependent activation of the APOE/C1/C4/C2 cluster supported lowering cholesterol levels (48) by triggering a defined network of genes.

Table 1.

Correlation of LXRα-binding data with GWAS

Function	Disease/term	SNP	P-value	Genes	Distance to SNP (kb)
Metabolism	LDL cholesterol	rs10401969	2.0E-08	PBX4	311.5
		rs4420638	1.0E-60	APOC1	5.7
		rs2142672	2.0E-08	MYLIP	52.1
	Triglycerides	rs17216525	4.0E-11	PBX4	57.0
		rs4420638	3.0E-13	APOC1	5.7
		rs7679	7.0E-11	PLTP	33.3
	HDL cholesterol	rs7679	4.0E-09		33.3
	Lipoprotein-associated phospholipase A2 activity and mass	rs1805017	6.0E-14	PLA2G7	40.4
Inflammation	Mean platelet volume	rs12485738	6.0E-31	ARHGEF3	16.0
	Rheumatoid arthritis	rs2062583	2.0E-06		84.4
	C-reactive protein	rs4420638	9.0E-139	APOC1	5.7
Heart	Myocardial infarction (early onset)	rs9982601	6.0E-11	MRPS6	17.2
	Coronary heart disease	rs9982601	4.0E-10		17.2
		rs7801190	3.0E-08	SLC12A9	53.6
		rs10512597	8.0E-11	CD300LF	135.3
	D-dimer levels	rs2774920	1.0E-06	ABCA4	164.3
	Resting heart rate	rs314370	6.0E-10	SLC12A9	58.5
	Diastolic blood pressure	rs1530440	1.0E-09	ARID5B	153.4
Brain	Alzheimer’s disease	rs4420638	2.0E-44	APOC1	5.7
	Information processing speed	rs6051520	2.0E-07	TRIB3	6.2
	Multiple sclerosis	rs17174870	1.0E-08	MERTK	93.7
	Amyloid A levels	rs2896526	4.0E-22	LDHA	92.5
Skin	Melanoma	rs3219090	9.0E-08	PARP1	72.9
		rs2230926	1.0E-17	TNFAIP3	36.6
	Psoriasis	rs610604	9.0E-12		33.2

For screening potential overlaps between LXRα-binding sites and validated SNPs that were genetically associated with common diseases, we chose the NHGRI GWAS catalogue. We considered only those LXRα target sites that were within LD to the reported SNPs. Columns indicate the general category, associated diseases or terms and the SNP ID. The P-value indicates the degree of certainty of disease association of the SNP. Next column depicts the reported and potentially disease relevant genes in proximity to the associated SNP. The last column shows the distance between LXRα binding sites and SNPs within the LD block, presented in kilobases (kb).

Figure 4.

Functional connection between LXRα binding to APOE/C1/C4/C2 locus and disease-relevant SNPs. (A) UCSC genome browser picture of the rs4420638 SNP and in linkage disequilibrium (DistiLD Database, Pallejà et al., 2011) surrounded area (chr19: 45247846-45544890). SNP position is indicated with a black arrow on top of normalized tag alignment tracks of LXRα ChIP-seq. LXRα-binding sites (black arrows under track) around apolipoprotein C1 gene (APOC1). LXRα enrichment in macrophages (grey), in T0901317 treated macrophages (brown), in foam cells (green) and in T0901317-treated foam cells (blue) is depicted. (B) STRING interaction network of genes within the extended APOE/C1/C4/C2-gene cluster with addition of well-known LXRα target genes. Differential expression of foam cells compared to macrophages is presented. Red circles indicate upregulated genes, and blue circles show downregulated genes. Grey lines represent the STRING network confidence view. Asterisk indicates LXR knockdown-sensitive genes. Bold circles represent genes with an enriched LXRα-binding site close by.

Atheroprotective potential of 32 novel LXRα target genes in T0901317-treated foam cells

The atheroprotective potential of the synthetic LXRα ligand T0901317 has been previously shown, but the underlying gene networks were only partly understood (49). First, we analysed the direct effects of LXRα binding and their further subsequent impact on foam cell development. Therefore, we selected differentially expressed and LXR knockdown-validated target genes in foam cells and generated a network with interaction data from the FANTOM and STRING database (Figure 5A). The introduced oxLDL stimulus led to activation of 160 genes and repression of 93 genes. The involvement of various transcription factors such as CCAAT-enhancer-binding proteins CEBPB and CEBPG or NF-κB1 illustrated once more the broad impact of ligand-induced LXRα activation on multiple metabolic processes regulated by LXRα itself, interacting transcription factors and other genes. For functional classification, we applied Ingenuity Pathway Analysis (IPA) and could annotate molecular functions, such as lipid metabolism, inflammation and cell death. The apoptotic function of these cells was already observed in the stringent set of LXRα target genes (Supplementary Table S2). In the next step, we performed differential expression analysis of T0901317-treated foam cells compared with untreated, diseased foam cells. We selected LXR knockdown-sensitive genes with LXRα binding site in this model and obtained 32 novel LXRα target genes and six previously described genes. We subjected all 38 genes to IPA and obtained the most enriched network functions in carbohydrate metabolism, molecular transport and lipid metabolism (Figure 5B and Supplementary Figure S6A). To gain more insight in the involvement of the underlying genes in classical LXR/RXR activation pathway (IPA top canonical pathway P-value 1.65e-3), we performed an overlap analysis using IPA (Figure 5C). We found two genes directly connected with LXRα (SCD, MYLIP, Supplementary Figure S6B) and 10 indirect interactions in this pathway, the latter are tightly bound on interactions with the nuclear factor of kappa light polypeptide gene enhancer in B-cells (NF-κB) complex and tumour necrosis factor-α (TNFA) signalling. Our previous analysis revealed all 38 genes to be direct LXRα target genes, indicating a strong bond between LXRα, NF-κB and TNFA signalling. Interestingly, arresten [COL4A1, collagen chain of basement membranes, (50)] and chemokine orphan receptor 1 [CXCR7, member of the G-protein coupled receptor family, (51)] were associated with cardiovascular disease (IPA top disorder, P-value 1.07e-3). Both genes were significantly upregulated after T0901317 treatment in foam cells. The highest upregulation (∼6-fold) showed the interferon-induced guanylate-binding protein 2 (GBP2, Supplementary Figure S6B). This GTPase has not been previously related to LXRα but was found to be associated with diseases such as psoriasis, rheumatoid arthritis and experimentally induced diabetes. One of the highest downregulated (∼2-fold) genes was protein inhibitor of activated STAT 4 (PIAS4, Supplementary Figure S6B), which interacts with the NF-κB complex (52).

Figure 5.

Beneficial effects of LXRα modulation with synthetic ligand T0901317 in foam cells. (A) LXRα target gene network in foam cells. The differential expression of LXR knockdown-sensitive genes is displayed in foam cells compared with macrophages. Well-known and physiologically relevant target genes in lipid metabolism as well as in inflammation and apoptosis are emphasized. Red circles indicate upregulated genes, and blue circles show downregulated genes. Grey lines represent the FANTOM and STRING network. Bold circles represent genes with an enriched LXRα-binding site close by. (B) T0901317-specific LXRα target genes in foam cells. The differential expression of LXR knockdown-sensitive genes in T0901317 foam cells is shown compared to foam cells. Red circles display upregulated genes, and blue circles indicate downregulated genes. Bold circles represent genes with an enriched LXRα-binding site close by. Grey lines represent LXR knockdown-validated interaction with LXRα. The top associated network functions assigned by IPA are carbohydrate metabolism, molecular transport and lipid metabolism. (C) LXR/RXR activation pathway with direct involvement of T0901317-specific LXRα target genes from associated IPA network in foam cells. The differential gene expression of knockdown-sensitive LXRα target genes on T0901317 treatment in foam cells compared to untreated foam cells are visualized in red for significantly upregulated genes and in blue for significantly downregulated genes. Unregulated gene expression of central gene products is displayed in white. Bold circles represent an enriched LXRα-binding site near the gene.

In summary, this analysis revealed a number of so far unknown factors of the artheroprotective network that were found to be enriched in functions, such as carbohydrate metabolism, molecular transport and lipid metabolism. These new genes can be of important value for thorough understanding and future investigations of transcriptional regulation in cellular atherosclerosis processes and ligand-specific activation of LXRα.

DISCUSSION

Comprehensive analyses are important to decode the complex molecular networks that underlie cell physiology and common diseases (53). In contrast to conventional molecular biological studies that usually focus on individual target genes of a nuclear receptor, we aimed to identify gene networks of LXRα that coherently respond to natural ligands or to pharmacological intervention in foam cells and macrophages. To make the analysis most stringent, and to largely exclude any potential off-target effects of LXR ligands, we considered in our network analyses only LXRα loci that showed significant effects on gene expression and were additionally sensitive to LXR knockdown.

In this study, we analysed human macrophages, which in stark contrast to mouse macrophages show a characteristic auto-regulatory loop of LXRα activation. In general, studies on atherosclerotic plaque development in mouse can only be partly adapted for understanding human atherosclerosis involving mainly smaller arteries instead of the large vessels that are analysed in mice (54). The human THP1 cell line applied here is a widely used macrophage model for human foam cell development (55). Consistent with numerous reports, we could largely observe similar effects in THP1 and primary macrophages on LXRα-binding or gene expression patterns. The standardised THP1 cell line allowed our comprehensive analyses, which would have been difficult with primary cells.

Nuclear receptors, such as LXRα, react to various environmental factors to adapt metabolic and other cellular pathways. With the rise of genome-wide binding studies, the classic model of nuclear receptor transactivation was challenged by recent findings (19,39,56,57). The alternative activation mode, observed in these studies, included a complex interplay between ligands, pioneer factors, co-regulators and histone modifications. According to the alternative activation mode, in absence of a specific ligand, the chromatin has a more closed structure and is less accessible for nuclear receptor binding, which changes rapidly upon ligand activation (58). In contrast to a mixed-binding pattern of basal and induced LXRα binding, we observed ligand requirement for all significant binding sites. The few observed ligand-independent, basal-binding sites were not significantly enriched and eliminated during filtering. This observation was confirmed by the overall low protein amount in macrophages in absence of a ligand. For example, using ChIP-seq we did not detect LXRα binding in absence of a ligand in the ABCA1 locus. To rule out sensitivity issues, we tested among others the ABCA1 and ABCG1 loci with ChIP–qPCR and observed that these loci were clearly increasingly occupied by LXRα upon T0901317 treatment. In absence of a ligand, we detected only spurious LXRα binding, which was mostly in the range of negative control binding sites. This result differs from recent findings from Jakobsson et al. (59), who applied a subtype unspecific LXR antibody and performed endpoint ChIP–PCR analyses, which may lead to different results in terms of accurate relative quantification compared with ChIP–qPCR using real-time detection as applied in this study.

As shown here, the evolutionary rather new nuclear receptor LXRα features comparably poor defined sequence specificity but increased flexibility with regard to natural ligand activation including compounds derived from food. Consistently, we detected differential LXRα-binding patterns in foam cells, in T09011317-treated foam cells that contained also oxysterols, and in normal macrophages. Dependent on the individual requirements of the macrophage or the foam cell, LXRα seemed to bind with individual intensity to target gene sites. In the case of foam cell-specific–binding sites, we detected less pronounced, almost DR4 motif-independent and promoter-distant–binding patterns. There are multiple mechanisms that could explain this specific recruitment, including stabilizing interactions with other transcription factors (39) via looping or ‘piggyback’ binding, or assisted binding through specifically modified histones with an open chromatin environment (60).

In contrast to the shared and T0901317-specific peaks at known lipid metabolism relevant target genes LXRα and APOC1, we discovered for the foam cell-specific peaks unexpected gene loci: for example the novel LXRα target gene PARP1, which represents an enzyme that is known to be involved in DNA damage repair. Recently, it has been shown that inhibition of PARP1 leads in particular in brown adipose tissue and skeletal muscle to enhanced mitochondrial metabolism via activation of the histone deacetylase SIRT1, which culminates in protection against metabolic disease (61). Interestingly, in diseased foam cells, we detected increased oxLDL triggered LXRα binding at the PARP1 locus but striking absence of LXRα binding at this locus under T0901317 therapy. Thus, LXRα modulation by an efficient synthetic LXR ligand may play a so far unexplored role in counteracting the effects of (basal) PARP1 gene binding effects of oxLDL-activated LXRα.

Recent changes in human life style and nutrition may overstrain the natural balance in foam cells to reverse cholesterol transport. The synthetic LXRα ligand was more efficient in triggering sequence-specific binding and thereby controlling more stringently LXRα transcriptional regulation. This genome-wide effect form increased atheroprotective effects. In summary, we identified disease-causing and pharmacologically treatable gene networks that can be modulated by specific intervention. Altogether, we found 32 novel LXRα target genes with atheroprotective potential. In contrast to natural oxysterols, T0901317 treatment resulted in striking binding of LXRα at the APOE/C1/C4/C2-gene cluster and its surrounding genes. This genomic region has been found to be highly significantly associated with the disease terms “LDL-cholesterol” and “C-reactive protein”, which suggests that specific ligand-induced LXRα binding—in conjunction with other transcriptional regulators—at this locus is important for establishing a key network of cholesterol transport and anti-inflammatory genes.

Quantitative LXRα-binding studies using clinical samples of diseased patients and control subjects may reveal sequence-dependent variation in LXRα binding efficiency, and thus explain the functional impact of genetic variation in the APOE/C1/C4/C2-gene locus for cardiovascular or other diseases.

We have generated LXRα genomic loci and related expression data in human foam cells, which may also provide a resource for developing further mechanistic studies, including the use of alternative small molecule activators of LXR, such as GW3965 (62) to analyse ligand-dependent fine-tuning of LXR-regulated gene expression via differential co-factor recruitment.

In summary, our data revealed an LXRα ligand-dependent network of transcriptional regulation, which is vulnerable to molecules that efficiently activate LXRα and shape gene expression patterns. Single gene analysis approaches may hit the limits of reductionism when biological function emerges on the interaction network level of genes or gene products (63). Although genome-wide DNA binding of LXRα may be deduced from binding or reporter gene assay efficiencies, expression can vary because of ligand-dependent differential interaction with further transcriptional co-factors leading to varying response on the gene network level (24).

Integration of gene, protein or metabolic networks (64,65), will potentially generate an additional layer of complementary information to provide an even more comprehensive understanding of complex interconnected molecular pathways in cell metabolism during atherogenesis and to optimize treatment of cardiovascular disease.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online: Supplementary Tables 1–3, Supplementary Figures 1–6, Supplementary Methods, Supplementary Datasets 1–4 and Supplementary References [66-70].

FUNDING

German Ministry for Education and Research [BMBF, 0315082]; National Genome Research Net [NGFN, 01 GS 0828]; the European Union [FP7/2007-2013], under grant agreement n° 262055 (ESGI). Funding for open access charge: Max Planck Society.

Conflict of interest statement. None declared.