The active enhancer network operated by liganded RXR supports angiogenic activity in macrophages

Here, Nagy and colleagues use genome-wide approaches to uncover the activity of RXR, an enigmatic member of the nuclear receptor superfamily. RXR signaling is predicted to have a major impact in macrophages, but neither the biological consequence nor the genomic basis of its ligand activation is known. Integrating RNA-seq, ChIP-seq, GRO-seq, and 3C-seq, the authors unravel the mechanism of RXR-induced transcriptional events in mouse bone marrow-derived macrophages. Importantly, this study uncovers a novel biological activity—angiogenesis—that is promoted by the receptor.

RXR is an unusual and somewhat neglected member of the nuclear receptor superfamily because it is still not known how this receptor interacts with the genome and regulates gene expression upon ligand activation. It is expressed in every cell type and is required for postnatal life in mice (Mangelsdorf et al. 1992; for review, see Szanto et al. 2004). Its presumed main molecular function is to regulate the activity of a dozen or so nuclear receptors. There is also evidence that it can form homodimers and/or have heterodimer-independent signaling capacity (Zhang et al. 1992;Szeles et al. 2010). A key concept regarding RXR signaling is the permissiveness/nonpermissiveness mutually exclusive dual paradigm. According to this, in certain heterodimers such as RXR:PPAR and RXR:LXR, ligand activation of RXR results in transcriptional activation; hence, these are permissive heterodimers, while in other heterodimers such as RAR:RXR, TR:RXR, and VDR:RXR, RXR is suppressed or ''subordinated,'' and therefore these so-called nonpermissive heterodimers cannot be activated from the RXR side (Germain et al. 2002). Therefore, the activation of all permissive heterodimers present in a particular cell type might lead to pleiotropic gene activation and engagement of potentially conflicting pathways.
The existence of pleiotropy, the role of RXR activation of permissive heterodimers, and the presence and activity of RXR homodimers have been debated and remain largely unresolved.
The fact that certain natural lipids-such as 9-cis retinoic acid (RA), docosahexanoic acid, and phytanic acid-are able to activate RXR gives further support to the biological role of RXR activation in vivo (for a review, see Szanto et al. 2004). There are also potent and selective synthetic compounds such as bexarotone (LG10069) and LG100268 (LG268) (Boehm et al. 1994(Boehm et al. , 1995 that have been used to dissect the role of the receptor in various biological systems and/or used in therapies. In macrophages, there are several heterodimeric receptors with key cellular roles, such as PPARg regulating oxLDL uptake and processing, LXR regulating cholesterol efflux and immune function, and NR4A1 (NUR77) regulating inflammatory response. These heterodimeric receptors have been linked to the development of atherosclerosis and also immune function and provide means to reprogram macrophages (for reviews, see Calkin and Tontonoz 2012;Nagy et al. 2012). Therefore, it is important to understand how activation of RXR contributes to these pathways and potentially to novel ones and regulates gene expression in this cell type.
A key issue in understanding any signal-specific transcription factor is the determination of the genomic sites to which it binds and to link those to the target transcripts. Recent advances in genome-wide approaches allow global assessment of histone modifications and transcription factor genomic binding sites (Ostuni et al. 2013). It is done using microarrays or RNA sequencing (RNA-seq) (Mortazavi et al. 2008). However, the temporal changes in steady-state mRNA levels diverge from changes in transcription; hence, simply determining steady-state mRNA levels will not allow easy identification of primary transcriptional events and will not separate those from secondary and tertiary ones (Bhatt et al. 2012). The development of genome-wide localization studies, primarily ChIP-seq (chromatin immunoprecipitation [ChIP] combined with deep sequencing) approaches, aided the determination of transcription factor-binding sites within the context of the chromatin structure on the genomic scale. Such genome-wide studies suggest that enhancers and promoters exhibit distinct chromatin ''signatures.'' The characteristic signature for enhancers consists of monomethylation of histone H3 Lys4 (H3K4me1), acetylation of histones (H3K9ac and H3K27ac), and binding of the acetyltransferase P300 (Heintzman et al. 2007;Koch et al. 2007;Visel et al. 2009).
In spite of all these developments, the reliable linking of activated enhancers to the regulated gene is still not solved, partly due to the lack of markers linking the genomic binding site to the regulated transcript and partly due to the potentially large distances between the enhancer and the regulated gene. Moreover, the identification of enhancers is even more cumbersome, if not impossible, in the case of signal-dependent transcription factors (such as RXR), which may constitutively occupy genomic sites. In our studies, we tried to solve these issues by combining RNA-seq, ChIP-seq, GRO-seq (global run-on sequencing), and 3C-seq (chromosome conformation capture [3C] combined with sequencing) in a highly integrated way to unravel the mechanism of RXR-induced transcriptional events in mouse bone marrow-derived macrophages (BMDMs) and, as the result of the process, discovered and validated a novel biological activity promoted by the receptor.

The transcriptional consequences of RXR activation in murine macrophages
We set out to systematically determine the genomic events following RXR ligand activation in BMDMs (Fig.  1A). The synthetic and selective RXR ligand LG268 was applied throughout the studies in a 100 nM concentration. We determined the changing transcripts by RNA-seq (analyzed by the pipeline shown in Supplemental Fig. S1A) and found that selective activation of RXR affects the steady-state mRNA levels of hundreds of genes, as expected. A hierarchical clustering of the top 200 changing genes as a time course is shown in Figure 1B. In order to identify primary RXR targets, we embarked on analyses to determine both the genomic binding events of RXR (RXR cistrome) and the enhancers responsible for the transcriptional changes.
Determination of the RXR, PU.1, and P300 cistromes and mapping active histone marks in BMDMs, the effect of ligand By carrying out ChIP-seq experiments, we determined the cistromes of RXR in the absence and presence of its ligand, LG268, and its relationship to the lineage-specific transcription factor PU.1 and the cofactor P300 along with a marker of transcription initiation, H3K4me3, and active histone marks, H3K27ac, H4ac, and H3K4me2, using the algorithm depicted in Supplemental Figure S1B. We determined ;5200 RXR genomic binding regions. Importantly, this cistrome is not impacted greatly by ligand treatment in 60 min. PU.1 occupied the highest number of peaks (;30,000), with only a minimal rearrangement upon ligand exposure (Supplemental Fig. S1C). Although the number of RXR-binding regions does not change much upon ligand exposure, peaks gain ;30% more reads, suggesting that RXR enrichment on the genomic regions is enhanced (Fig. 1C). We confirmed these observations by ChIP and real-time quantitative PCR (RT-qPCR) on binding regions of known directly regulated genes ( Fig. 1D; Nagy et al. 1996;Kersten et al. 2000;Repa et al. 2000). Intriguingly, P300 binding follows the genome-wide RXR enrichment, suggesting that P300 is likely to be recruited to the genome upon RXR activation (Supplemental Fig. S1C).
The motif rank order under the detected RXR peaks reported that PU.1 had the most enriched motif, followed by various combinations of repeats, including DR1 and DR4, of the nuclear receptor-binding (half) site. In addition, with somewhat lower abundance, two other macrophage-associated motifs were also detected: C/EBP and AP-1 (Fig. 1E). As far as the genomic distribution of the detected peaks is concerned, ;90% were found outside of the potentially directly regulated 823 genes identified by RNA-seq experiments (i.e., not within 10 kb of their transcription start sites [TSSs]) (Fig. 1F), suggesting that linking binding sites to regulated genes based simply on proximity will be difficult.
Next, we continued to interrogate the cistromes of PU.1, RXR, and P300 in the ligand-activated state. The RXR peaks overlap with both PU.1 and P300 to a large extent and show increased acetylation ( Fig. 2A,B). Global analyses show that ligand activation of RXR leads to P300 recruitment (Fig. 2B,C). These changes could be validated by ChIP-RT-qPCR (Fig. 2C). Importantly, global P300 binding and recruitment could be detected and confirmed at >400 sites using three biological replicates, and no decrease could be observed (Supplemental Fig. S2A). As far as active histone marks are concerned, H4 acetylation and H3K4 dimethylation are increasing, while H3K27 acetylation remains largely unchanged upon ligand activation on the RXR-bound genomic regions (Fig. 2D).
Next, we sought to determine whether some of the potential partner ligands had a detectable effect on the RXR cistrome. Therefore, we determined the RXR cis-tromes in the presence of a synthetic LXR (GW3965) or PPARg (RSG) ligand as the activators of the two major suspected heterodimeric partners present in macrophages. We found that neither ligand had a statistically significant effect on RXR's genome-wide distribution and/or affinity, applying a <0.1 false discovery rate (FDR) threshold (Fig.  2E). Finally, we wanted to evaluate the possibility of whether an endogenous RXR ligand is masking a ligandinduced (re)distribution of the receptor. We reasoned that eliminating the effect of a presumed endogenous ligand by an RXR antagonist (LG1208) would reveal the true unliganded RXR cistrome. However, as shown in Figure 2E, the RXR cistrome does not change upon antagonist treatment either. These findings were confirmed using RT-qPCR (Supplemental Fig. S2B). These data suggest that the LG268. Differentiated macrophages were treated with 100 nM LG268 for the indicated time, and total RNA was subjected to RNA-seq gene expression analysis (Supplemental Fig. S1A). Biological replicates had at least 70 million reads each. The presented gene set was filtered at FDR < 0.1 and a more than twofold change. (C) Average read distribution of RXR peaks in the presence or absence of LG268. ChIP-seq data were analyzed as shown in Supplemental Figure S1B (peaks were predicted by MACS2, and the consensus of the 2 3 2 replicates was determined by DiffBind); the histogram was generated by HOMER. (#) Normalized throughout the study. (D) RXR binding on the indicated individual enhancers confirmed using ChIP-RT-qPCR. Macrophages were treated with 100 nM LG268 for 1 h. The mean and 6SD of triplicate determinations are shown. (E) De novo and targeted (asterisk) identification of motifs under RXR peaks from ChIP-seq data using HOMER. ''Target %'' refers to the ratio of the peaks having the given motif, and ''Bg %'' refers to the ratio of a random background as described in the Supplemental Material. (F) The genome-wide distribution of RXR peaks relative to the TSS of the closest regulated genes identified by RNA-seq (P < 0.05; n = 823). RXR cistrome is not changing upon short-term ligand activation and that the RXR selective ligand only increases the receptor's enrichment at its predetermined binding sites, which typically bind PU.1 and recruit P300. Importantly, the dominant ligands of relevant permissive heterodimers such as PPARg or LXR do not appear to have an effect on either distribution or affinity.

Determination of liganded RXR regulated nascent RNA production
The merging of the ChIP-seq and RNA-seq data sets proved to be insufficient to reliably identify the network of primary regulated genes and their enhancers. This is due to the fact that the physical location of the bound transcription factor cannot be linked to the regulated gene. Almost 90% of the RXR peaks are located outside of the 10-kb region upstream of the TSS of the closest regulated gene (Fig. 1F), which is in good correlation with others' observations (Thurman et al. 2012).
We reasoned that in order to link a liganded RXRoccupied enhancer to the corresponding regulated gene, one needs a high-resolution method, which can provide dynamic temporal information about transcription. The recently recognized fact that nascent RNA production can be detected on both the regulated gene , and H4ac occupancies in 3-kb windows around the summit of the RXR peaks in the presence or absence of LG268. Read distribution was determined by HOMER, clustering was done by Gene Cluster 3.0 using centered correlation similarity metric with single linkage clustering method, and heat maps were created by Java TreeView in log 2 scale. (C) Read distribution of P300 on RXR peaks in the presence or absence of LG268. P300 binding was confirmed on the indicated individual enhancers using ChIP-RT-qPCR. Macrophages were treated with 100 nM LG268 for 1 h. The mean and 6SD of three biological replicates are shown. Asterisk represents significant difference at P < 0.05; n = 3. (D) Average read distribution of the indicated active histone marks on RXR peaks in the presence or absence of LG268 determined by HOMER. (E) RXR enrichments of the significantly changing peaks in the presence of the indicated ligands. Cistrome was determined in the presence or absence of LG268, and the log 2 -normalized read numbers of the significantly changing peaks were plotted by DiffBind. The number of changing peaks and the statistical stringency applied is indicated below each plot. and the enhancer offered a plausible solution .
Therefore, we decided to detect the dynamics of nascent RNA production using GRO-seq. As shown in Figure 3A, we carried out a time-course experiment to determine the sites and dynamics of nascent RNA production upon LG268 activation, a typical example of which is shown for Abca1 induction in Figure 3B. Importantly, we noticed the presence and induction of nascent RNA transcripts at sites of enhancers, eTranscripts, and also short divergent transcripts (Core et al. 2008). We collectively call these enhancer RNAs (eRNAs). Notably, if one aligns the detected GRO-seq activities with cistromic data, it is easy to recognize that short transcripts overlap with transcription factor-bound regions-in our case, PU.1, RXR, and P300 ( Fig. 3C; Supplemental Fig. S3A). A comprehensive analysis of GRO-seq active regions, including gene transcript and eRNA calling, was carried out by a unique algorithm (Supplemental Fig. S3B). We determined that ;20% of the mouse genome is transcribed (with at least 0.006 RPKM [reads per kilobase per million mapped reads] expression) in this macrophage cell type and identified 10,586 known genes also marked by H3K4me3 and, all together, 25,560 transcripts, which include eRNAs and other noncoding RNAs as well (Supplemental Fig. S3C). Importantly, we found 51,657 GRO-seq divergent sites (>0.2 RPKM expressed) characterized by divergent transcription with a distance <300 bases between the 59 ends of the transcripts on the two strands. These sites do not necessarily show elongation, just initiation irrespective of transcribed genes and other longer transcripts. One could reasonably assume that this set of genomic regions contains all active promoters and enhancers. Using this approach, we determined changing (induced/repressed) gene transcripts and eRNAs as well. If one plots the changing levels of nascent RNA production of already established direct target genes, the dynamics of the changes are indicative of an immediate induction (Fig.  3D). This also allows easy identification of additional direct targets.
These data further reinforced the notion that the dynamically changing nascent RNA landscape provides clues about direct transcriptional regulation specific for RXR liganding. Furthermore, this allowed us to generate a list of genes most likely being directly regulated by liganded RXR (Supplemental Table S1).
As far as the active enhancers are concerned, we classified the identified ;5200 RXR-binding sites into two categories using an algorithm (pipeline depicted in Supplemental Fig. S4A): (1) enhancers that overlap with a GRO-seq-positive region (divergent site; 2781) and (2) enhancers that do not overlap with GRO-seq activity (2425), as shown in Figure 4, A and B). This classification showed an enrichment of RXR peaks within the proximity of the closest TSS (<10 kb from the TSS), suggesting that our assumption was likely correct and that we were identifying the functional/active binding regions and efficiently separating those from the silent/nonactive/ parking ones. The distribution of GRO-seq-positive and GRO-seq-negative RXR-bound regions relative to the nearest TSS shows a bias toward upstream regions in cases of the positive ones. Most interestingly, the enriched motif distribution under these peaks also showed characteristic differences. The GRO-seq-positive RXR sites showed an enrichment for DR1 and DR4 (binding sites of RXR:PPAR and RXR:LXR heterodimers, respectively) as well as AP-1 sites when compared with the GRO-seqnegative sites (Fig. 4B).
GRO-seq-positive RXR-binding sites are characterized by an increased H3K27ac, H4ac, and also H3K4me2 histone marks when compared with negative ones (Fig.  4C,D). These provided further support to the notion that these sites are indeed functionally distinct and likely represent active enhancers.
A closer look at the up-regulated GRO-seq-positive regions revealed that those that overlap with RXR peaks overlap with a significant fraction of PU.1 regions as well (Fig. 5A, top). The regions of down-regulated GRO-seq activity have much fewer RXR peaks but show overlap with PU.1 regions to a larger degree as well (Fig. 5A, bottom). These latter analyses showed that there are very few liganded RXR-occupied negative binding sites or silencers, suggesting also that liganded RXR is predominantly a transcriptional activator. Next, we matched up the divergent GRO-seq sites with RXR sites (Fig. 5B;  Supplemental Fig. S5A). The identified 51,657 divergent GRO-seq sites likely contain all TSSs, short divergent transcripts, and eRNAs. More than 3300 of these change upon RXR ligand activation, and 718 overlap with RXR-binding sites as well (Fig. 5A,C). We identified 252 regulated genes to which we could assign 414 RXR-bound regions, 387 enhancers, and 27 silencers using the criteria set by us (Fig. 5B,C; Supplemental Fig. S4A).
The distribution of the various motifs within these potentially active enhancers show that ;50% can be classified as a nuclear receptor-binding site repeat, another 20% is PU.1 only, and a quarter is other or unknown (Fig. 5D). The number of enhancers per gene ranges from one to nine (Fig. 5E). The distribution of these putative enhancers relative to the TSS is fairly symmetrical, and they span long distances, with only 16.54% found within 10 kb upstream of the TSS (Fig. 5F). The identified 387 enhancers have a pronounced increase in RXR occupancy, P300 recruitment, and H4ac, while H3K4me2 and H3K27ac are not different upon ligand treatment (Supplemental Fig. S4B).
Taking these together, we (1) uncovered the core active enhancer network operated by liganded RXR; (2) separated these sites from nonfunctional, silent/parking RXR-binding sites; and (3) paired the active enhancer network with the regulated genes using a set of criteria.

Functional validation of novel distant and long-range enhancers
Before we embarked on molecular validation of the enhancers, we intended to assess the contribution of the activated partner receptors to the induction of various RXR-regulated genes and eRNAs also. Therefore, we determined the induction of steady-state RNA levels of selected genes. We found that Vegfa and Hbegf were induced primarily by RXR liganding with some activity by RAR, Tgm2 and Ccl6 were both induced by RXR and RAR liganding, Abcg1 was induced by RXR and LXR, and Angptl4 was induced by the RXR and PPARg ligands (Supplemental Fig. S5).
Interestingly, primarily RXR-induced eRNA production could be detected on an enhancer assigned to Vegfa or Tgm2, while an enhancer of Abcg1 also showed robust LXR ligand activation, as expected (Fig. 5G). These data suggested that the eRNAs can be easily validated and show ligand induction similar to the regulated genes and therefore most likely are linked.
For functional validation, we chose 45 newly identified regulatory regions assigned to 23 genes (Supplemental Table S1). We selected 30 distant enhancers (>10 kb from the TSS), which earlier studies performing traditional enhancer analyses (Schwartz et al. 2000) would likely have missed, including one of ours (Nagy et al. 1996). In addition, we sought to validate 15 more proximal enhancers (<10 kb from the TSS).
Next, we cloned the 45 putative cis-regulatory elements (35 enhancers and 10 silencers) in front of a luciferase reporter gene and measured their transcriptional activity in COS1 cells in the presence of combinations of receptors and ligands. We clustered the enhancers based on the induction patterns obtained in the transient transfection-based reporter system, which is detailed in the Supplemental Material (Supplemental Fig. S6). These analyses show that the identified and tested enhancers are RXR-specific enhancers but also complex and very versatile, allowing combinatorial and context-dependent regulation of the genes by distinct dimers. The context is determined by the receptor expression profile and probably also lineage-specific transcription factors and cofactors. In addition, a single gene is likely to have distinct types of enhancers.  A). Columns represent the peak number of the 1-kb distance bins. (Insert) De novo and targeted (black asterisk) identification of motifs under GRO-seq-positive (left side) and GRO-seq-negative (right side) RXR peaks determined using HOMER. Differentially represented motifs are marked by red asterisks. ''Target %'' refers to the ratio of the peaks having the given motif, and ''Bg %'' refers to the ratio of a random background as described in the Supplemental Material. (C) Heat map representation of RXR, P300, H3K4me2, H3K27ac, and H4ac occupancies in 3-kb windows around the summit of the GRO-seq-positive and GRO-seq-negative RXR peaks in the presence or absence of LG268. Read distribution was determined by HOMER, clustering was done by Gene Cluster 3.0 using centered correlation similarity metric with single linkage clustering method, and heat maps were created by Java TreeView in log 2 scale. (D) Read distribution of the indicated histone marks relative to the GRO-seq-positive and GRO-seq-negative RXR cistromes determined by HOMER.

The genome architectural context of RXR-regulated gene expression
Once we determined the RXR-regulated genes and linked at least a subset of them to their enhancers (some of them are long-range ones), we wanted to get some insight into whether these regulatory units are confined to known features of genome architecture. Therefore, we decided to explore the genome architectural context of RXR signaling, asking the question of whether RXR signaling is confined by known structural features of genome organization and/or contributes to reorganization of these. We determined the cistromes of CTCF and a member of the cohesin complex (RAD21) in the absence or presence of ligand-activated RXR. We found ;30,000 peaks for CTCF and ;24,500 for RAD21 (Supplemental Fig. S1C); the latter showed some increase upon ligand activation. Next, we determined an overlap of ;12,660 peaks representing similarly high CTCF and RAD21 occupancy (with peaks having a MACS2 score >15 and having less than threefold difference). We considered these as boundaries of functional domains, as suggested by others (Merkenschlager and Odom 2013;Sofueva et al. 2013). We identified 10,204 such functional domains by pairing the CTCF and RAD21 copeaks with the closest ones with similar scores (Supplemental Material).
The median length of such domains was 81.15 kb, while the average was 149.98 kb. The neighboring domains then were merged into active domains with <100-kb distances between them. We could designate almost 700 such domains, which have a median length of 1.1 Mb and an average length of 1.88 Mb. If one overlays this architectural domain structure with active RXRbinding sites and regulated genes that we identified (Fig.  5C), one can find that 203 out of the 252 (80%) identified RXR-regulated genes along with their enhancers fall on such functional domains. In addition, one can identify ;1113 CTCF genomic binding regions on which RAD21 (cohesin) binding is enhanced and 128 regions on which RAD21 binding is reduced upon RXR ligand treatment (Fig. 6A, left). Moreover, RAD21 (cohesin) binding is more enhanced on active RXR enhancers (GRO-seq + ones) when compared with nonactive ones (GRO-seq À ones) (Fig. 6A, middle and right). Importantly, 40% (84 out of 203) of regulated genes have induced RAD21 on their enhancer and/or CTCF-binding sites. Collectively, these data suggest that most if not all RXR enhancers act inside functional domains and that activation of the receptor contributes to formation of such domains by stabilizing the genomic architecture and, in some cases, even inducing the binding of enhancers to promoters. Next, we validated molecularly the interactions in such functional domains and the impact of ligand on these for long-range enhancers of Abcg1, Vegfa, and Tgm2 using RT-qPCRbased 3C (Fig. 6B). These data clearly documented that the enhancers identified by the combination of RNA-seq, ChIP-seq, and GRO-seq act in functional domains, loop to the promoter, and can be readily functionally validated using transient transfection and 3C as well.
Finally, we asked the question of whether any of these enhancers communicate with other cis-regulatory elements or functional domains in the genome. Thus, we carried out 3C-seq using pairs of baits located in or close to these regions. We could detect proximal and also longrange interactions (for details of the analyses, see the Supplemental Material). At the Vegfa locus, we could detect interactions between the distant enhancer and the neighboring enhancers as well as the intronic region with remarkable specificity (Fig. 6C). Similar results were obtained in the case of the Tgm2 and Abcg1 loci (Supplemental Fig. S7A,B). Importantly, we also detected interchromosomal interactions with much less frequency (at least 50-fold less) though. In order to provide statistical context to these findings, we compared the interaction frequency of a given bait with inactive topological domains, active domains devoid of RXR-regulated regions, and active domains with RXR-regulated regions. As shown in Figure 6D in the case of Abcg1 and Vegfa, there is significant difference between the frequencies of such interactions. Similar results were obtained with additional enhancer of these and the Tgm2 gene. This suggests that active RXR enhancers interact with other active genomic regions and with an even higher likelihood with other RXR-regulated ones. These formally suggest that the active RXR enhancers form an interchromosomal hub or network.

Determination of the impact of RXR activation on the angiogenic capacity of macrophages
Finally, we wanted to see whether some of the identified novel transcriptional pathways could be validated biologically. The functional annotation of the genes controlled by liganded RXR enhancers assigned a large number of genes into the angiogenesis category (Supplemental Fig.  S4C). These include Vegfa, Hbegf, Litaf, and Hipk2. Therefore, we decided to functionally test this activity.
First, we determined the excreted VEGFa protein levels from cell supernatants (Fig. 7A) using different ligands and could show that RXR induced this protein, as was suggested by our gene expression measurements as well (Supplemental Fig. S5). Using RT-qPCRs, we could show the RXR requirement using RXRa/b double-knockout macrophages (Fig. 7B). To test whether activation of RXR promoted angiogenesis in an in vivo relevant setting, we performed a chorioallantoic membrane (CAM) assay using macrophages pretreated with RXR agonist. These cells showed significantly increased angiogenic activity, which was not detectable in RXRa/b double-knockout macrophages (Fig. 7C). These data suggest that RXR activation can program macrophages toward a distinct cellautonomous angiogenic phenotype (Fig. 7D,E).

Discussion
The ultimate goal of studying the function of a particular transcription factor is to discern its primary biological action in the entire genome comprehensively. In principle, this should be done by defining the genomic binding  B1 enhancer (chr17: 31,172,998-31,173,656) and Vegfa B1 enhancer (chr17: 45,890,890,829) determined by 3C-seq. The interchromosomal interactions of these enhancers were determined as described in the Supplemental Material. The detected interactions were mapped onto 1-Mb fragments covering the mouse genome. Interaction frequency was determined by expressing the read number per 1 Mb normalized to 1000 reads. Genome regions (inactive, active, and RXR target) were qualified as described in the Supplemental Material. Unpaired two-tailed t-test analysis was used to determine significant differences. Asterisk represents significant difference at P < 0.0001. Figure 7. RXR activation leads to increased angiogenic activity in macrophages. (A) VEGFa protein levels determined from supernatants of macrophages treated with LG268 for the indicated periods of time determined by ELISA. The mean and 6SD of triplicate determinations are shown. Asterisk represents significant difference at P < 0.05; n = 3. (B) Vegfa mRNA levels determined using RT-qPCR in wild-type (WT) and RXRa/b knockout (KO) macrophages. The mean and 6SD of duplicate determinations are shown. (C) Macrophage angiogenesis activity was determined using CAM assay in the presence or absence of LG268 using wild-type and RXRa/b knockout macrophages. Representative images (left) and bar chart of quantification (right) are shown. The mean and 6SD of three biological replicates are shown. Asterisk represents significant difference at P < 0.05; n = 3. (D) The molecular mechanism through which Vegfa is regulated in macrophages. Vegfa harbors a set of very-long-range enhancers. The enhancers are marked by the lineage-determining transcription factor PU.1 and are able to recruit RXR presumably as a heterodimer with RAR or in a homodimeric form. Upon ligand stimulation, P300 is recruited, and eRNA production is highly increased. These elements can get into the close proximity of the gene promoter by looping, which is enhanced in the presence of the RXR activator (LG268) and leads to increased gene expression. This regulatory unit is bordered by CTCF/cohesin (RAD21) cobinding, which is thought to contribute to the topological domain structure of mammalian genomes and stabilizes chromatin loops. (E) Schematic representation of the various nuclear receptor dimers and their biological significance in macrophages. PPAR/RXR are known to inhibit the inflammatory response upon activation by either side as a permissive heterodimer, while the LXR/RXR heterodimer regulates inflammatory response, cholesterol metabolism, and triglyceride biosynthesis. Our results shed light on the effect of putative, permissive, RAR/RXR, or possibly RXR/RXR dimers in the angiogenesis program of macrophages, most probably through the regulation of Vegfa and/or Hbegf. regions and determining the genes proximally regulated in a given cell type. Using computational tools, one can identify potential RXR-binding sites on the scale of several times 100,000 in the mouse genome. When all of the ChIP-seq data available from fat, from liver, and now in macrophages are combined, the combined number of genomic regions bound by RXR is a few times 10,000. Still, in macrophages, the number of genomic regions is ;5200, but the number of regulated genes by liganding the receptor is a few times 100. This is further complicated by the fact that primary targets are hard to reliably identify by simply measuring steady-state mRNA levels. So the question remains: How can one find the regulated genes and the linked receptor-bound enhancers? We attempted to do this by integrating genome-wide analyses.

Integration of genome-wide localization and nascent RNA production data identifies active enhancers
A key premise of the work presented here is that genomewide localization studies, predominantly ChIP-seq experiments, can be intersected with data from nascent RNA determinations such as GRO-seq time courses, which then enables one to comprehensively annotate cistromes and identify active enhancers (Hah et al. 2013;Lam et al. 2013). This relies on the assumption that in vitro (nuclear run-on) determination of the activity of RNA polymerases (RNAPs) is a reliable indicator of transcriptional activation. This has been proven in other systems and leads to the identification of active enhancers in other cell types (Koch and Andrau 2011;Wang et al. 2011;Bonn et al. 2012), although systematic comparisons or crossvalidation with other methods such as subcellular RNA fractionation (Bhatt et al. 2012) has not been done. The usage of GRO-seq for the determination of direct transcriptional responses has two major advantages: (1) The dynamics of the nascent RNA levels depends only on the rate of RNAP activity, and therefore it is matching the expected time course of a directly regulated gene. Our data are clearly demonstrating this because known and established direct target genes (Abca1, Abcg1, Angptl4, and Tgm2) as well as newly identified ones (i.e., Vegfa and Hbegf) show an immediate induction when assessed by GRO-seq rather than a complex, often delayed one determined by RNA-seq or RT-qPCR (Fig.  3D). (2) Enhancers contain engaged RNAPs, and their activity can be revealed by GRO-seq in the form of a typically few-hundred-base-pair-long divergent transcript, termed eRNA. The existence of such transcripts has been shown before (Kim et al. 2010;Wang et al. 2011). Although the mechanism of their production and role in transcription is far from being clarified, it can be reasonably assumed that they represent active enhancers, as the presence of an engaged RNAP strongly indicates it. Furthermore, if the changes of nascent RNA production (increase or decrease) at these sites show correlation with the stimulating signal, it is likely to be associated with or the consequence of the activity of the signal-specific transcription factor. A key limitation of today's technology is the lack of high-throughput validation of the identified enhancers. Nonetheless, our high rate of success in validating such enhancers with low-throughput enhancer trap approaches, 3C-seq and 3C-RT-qPCR, suggests that signalspecific changes in eRNA production can be used to filter active enhancers from the general pool of genomic binding regions. Altogether, this also means that this combined approach can serve as proof of concept and a model to tackle similar problems with other signal-specific transcription factors. A remarkable feature of the approach is that it allowed the identification of active enhancers in spite of the fact that these represent a fraction of the binding regions (not more than 10%-15% by using our stringent criteria), including very-long-range enhancers. It is also interesting that the enhancers' distribution relative to the TSS appears to be rather symmetrical instead of being much biased to the upstream regions, suggesting that one needs to look in both directions and also far away to identify the true enhancer controlling a particular gene in a given cellular context. A cautionary conclusion of this study is that it is very likely that many of the enhancers identified with less comprehensive methods might need to be revisited.
RXR colocalizes with PU.1 and recruits P300 as a cofactor to its binding sites Our data showed that RXR has >5200 binding regions in macrophages. These peaks are likely to be RXRa-bound because this is the dominant subtype present in macrophages. However the used antibody is pan-RXR-specific and would recognize all three subtypes. Most of these sites contain nuclear receptor half sites (AGGTCA), suggesting that the localization of RXR uses nuclear receptor-mediated direct DNA binding. In addition, two independent lines of evidence suggest that there is an intimate relationship between the lineage-determining factor PU.1 and RXR: (1) 45% of the RXR peaks contain PU.1 sites, and (2) there is a significant overlap between RXR and PU.1 cistromes, as more than two out of three of RXR-binding regions overlap with PU.1 peaks. This is in agreement with our anticipation based on previous reports by Lazar and colleagues (Lefterova et al. 2010) and Liu and colleagues (Pott et al. 2012), who suggested that binding regions for PPARg, a heterodimerization partner of RXR, colocalizes with PU.1. This finding supports the pioneering or bookmarking factor concept put forward by several investigators (for a review see, Zentner and Scacheri 2012) to explain various interactions and genome-wide colocalization between lineage-specific and signal-specific transcription factors, suggesting that the lineage-specific factor opens up or marks particular regions in the genome, which then allows or even facilitates the binding of the signal-specific factors. Although the sequence of events or the mechanisms of such interactions are unknown, our data are fully compatible with such a scenario. An additional line of evidence supporting this is that there is only a 14% overlap between our RXR-binding regions and the ones found in 3T3L1 cells differentiating into adipocytes (Nielsen et al. 2008). Our work further extended these studies and suggests that many if not all RXR heterodimers at least in part colocalize with PU.1 in macrophages and that upon ligand treatment, PU.1 is at least partially released and therefore might not be required for further binding or transcriptional activity.
Another unexpected aspect of the cistromic interactions identified is that P300 binding is dynamically recruited upon the RXR ligand treatment to most of the active RXR sites. Two conclusions can be drawn from these findings: (1) P300 is major cofactor of mediating transcriptional activity by RXR in BMDMs in the steady state.
(2) Liganded RXR-regulated transcription complexes use P300 as one and potentially the major cofactor, and thus P300 is likely to be responsible for the acetylation of histone H4K5/8. This is consistent with the initially proposed role for this protein in nuclear receptor signaling (Chakravarti et al. 1996) and with the more recent findings (Jin et al. 2011) but also suggests that P300 is likely to be specific for this signaling pathway in this cell type.
A key issue with signal-specific transcription factors is whether the signal contributes to the (re)distribution of the particular transcription factor. In case of the RXR-containing receptor dimers, the accepted view is that RXRs bind to the genome in both the absence and presence of ligand (Boergesen et al. 2012). Our findings further support this view now by adding a genome-wide perspective and show that the only effect RXR ligand has on the genome-wide distribution of the receptor is that the enrichment on the preformed binding sites increases. These data obviously have to be interpreted within the context of ChIP experiments, meaning that the sum of all binding in all cells is detected in the time frame of the 40-min cross-linking used. This method does not allow the construction of a more dynamic and/or higherresolution picture of the receptors' activation and mobility. Methods using shorter time resolution show a much more dynamic behavior of the receptor though (Brazda et al. 2014). Therefore, our interpretation of the data is that RXRs' genome-wide localization is determined by its own DNA-binding capacity, which allows it to find preformed sites in the genome to which it can bind, and this is facilitated by additional factors such as a lineagespecific factor; i.e., PU.1. The increased affinity that is observed is either the reflection of more cells being involved in the response or a higher affinity for the binding site. However, the experimental approaches used are not able to determine the relative contribution of these two mechanisms. We also excluded the possibility that an endogenous ligand plays a major role in directing RXRs to its genomic binding sites by using an antagonist. This behavior is in stark contrast to that of steroid hormone receptors such as GR or estrogen receptor (ER), whose genomic binding is dictated by the addition of ligand (Carroll et al. 2005;Welboren et al. 2009). As far as the genome architectural context of the RXR-regulated enhancers and the regulated genes are concerned, we could show that the vast majority of enhancers and their regulated genes are confined to functional domains bordered by CTCF/cohesin (RAD21) complexes, suggesting that their activity localized to these loops (Supplemental Fig. S8).
In addition, the rearrangement of such loops is moderately though but impacted by ligand activation. Curiously though, our 3C-seq analyses also revealed that the functional domains and the enhancers within them interact with each other even on different chromosomes. The functional significance of such interactions is not clear but might suggest the existence of RXR-or partner-specific transcription foci or factories inside the nucleus.

Characteristics of the liganded RXR-operated enhancer network
Our data show, as expected, that the two major permissive heterodimers in BMDMs are RXR:LXR and, to a lesser extent, RXR:PPAR, as shown by the ready induction of their established target genes, Angptl4, Abca1, and Abcg1, and the motifs identified within their RXR peaks being RXR:PPAR-specific (DR1) and RXR:LXR-specific ones (DR4). However, our analyses also uncovered a set of regulated genes, including Tgm2, and novel ones, such as Vegfa and Hbegf, which could not be induced efficiently by ligands for permissive heterodimers, only the RXR-selective LG268 or ligands activating the RAR receptor. This raises the possibility of the existence of permissive RXR:RAR heterodimers or a complex regulatory mechanism/enhancer allowing activation by either RAR or RXR-activating ligands (Fig. 7D).
Underpinning this complexity, the enhancers identified for the regulated genes show a large degree of functional versatility and can be grouped into four broad categories. Cluster A requires the presence of liganded RXR; another (cluster B) is best activated if both RXR and one of its partners is expressed; a third (cluster C) is mediating RXR-specific signaling, provided RXR is expressed at a significantly higher level than its partners; and the enhancers in the fourth (cluster D) work with multiple combinations of dimers (Supplemental Fig. S4). A further level of complexity is that a given gene (i.e., Vegfa or Abca1) has enhancers from more than one category.
All of these support that RXR has a unique genomic effect that cannot be recapitulated by any ligand or combination of ligands. This is mechanistically served by activation of permissive heterodimers such as LXR:RXR and PPAR:RXR and activation of enhancers inducible by RAR as well. In addition, some of the eTranscript inductions appear to be specific for RXR.
This and the facts that, using our most stringent criteria, we identified only 226 induced genes, and only 7.4% of the identified binding sites are active enhancers suggest that activation of RXR is tightly controlled in this cell type and leads to distinct and selective and not pleiotropic gene expression. Unexpectedly, a network of enhancers that could be linked to a set of genes-including Vegfa, Hbegf, Hipk2, Litaf, Cxcl2, and Foxo3, which have been clearly linked to angiogenesis-was also revealed.
The regulation of Vegfa by RXR is intriguing because it appears to use a set of very-long-range enhancers (270 kb downstream) and shows RXR/RAR specificity, including RXR ligand-induced looping (Fig. 7D). We presented here a novel integrated approach to identify functional enhancers and link those to primarily regulated genes for a signal-dependent transcriptional factor, liganded RXR. This integrated approach revealed that ligand stimulation of RXR activates only a small fraction of the DNA-bound molecules confined by CTCF/cohesin-delimited functional domains (Supplemental Fig. S8) and leads to a distinct gene expression program, which results in increased angiogenic potential and might be a valid macrophage reprograming/therapeutic target (Fig. 7E).

Materials
The following ligands were used: LG268 and LG1208 were gifts from M. Leibowitz (Ligand Pharmaceuticals), RSG and AM580 were from Sigma, and GW3965 was a gift from T.M. Wilson (GlaxoSmithKline).

Differentiation of BMDMs
Isolation and differentiation were completed as described earlier (Barish et al. 2005).

RNA-seq
The RNA-seq library was prepared from two biological replicates by using a TruSeq RNA sample preparation kit (Illumina) according to the manufacturer's protocol. Analysis was carried out as described in the Supplemental Material.

ChIP
ChIP was performed as previously described (Barish et al. 2010), with minor modifications. Libraries were prepared by Ovation Ultralow Library Systems (Nugen) from two biological replicates according to the manufacturer's instructions. ChIP-seq analysis and domain prediction were carried out as described in the Supplemental Material.

GRO-seq
GRO-seq and library preparation were performed as described earlier (Core et al. 2008;Hah et al. 2011), with limited modifications. Libraries were generated from two biological replicates of BMDMs treated with 100 nM LG268. Analysis was carried out as described in the Supplemental Material.

RT-qPCR
RNA was isolated with Trizol reagent (Molecular Research Center). RNA was reverse-transcribed with Tetro reverse transcriptase (Bioline). Transcript quantification was performed by RT-qPCR reaction using SYBR Green dye (Diagenode). Transcript levels were normalized to Ppia or Rplp0.

3C
3C experiments were completed as described previously with minor modifications (Miele et al. 2006). The detailed protocol is available in the Supplemental Material.

3C-seq
Experiments were carried out as previously described (Stadhouders et al. 2013). For details, see the Supplemental Material.

Transient transfection
Enhancer sequences were PCR-amplified from a BAC or genomic DNA and cloned into pUC18 HSV TK-LUC. Transient transfections were carried out as previously described .

ELISA
ELISA experiments were carried out according to the manufacturer's instructions (R&D Systems).

CAM assay
CAM assays were performed as described previously (Movahedi et al. 2008). The detailed protocol is presented in the Supplemental Material.

Data access
Sequencing data were submitted to Sequence Read Archive under accession number SRP019970.