Dissecting the retinoid-induced differentiation of F9 embryonal stem cells by integrative genomics

doi:10.1038/msb.2011.73

Mol Syst Biol. 2011; 7: 538.

Published online 2011 October 11. doi: 10.1038/msb.2011.73

PMCID: PMC3261707

Dissecting the retinoid-induced differentiation of F9 embryonal stem cells by integrative genomics

Marco A Mendoza-Parra,¹ Mannu Walia,¹ Martial Sankar,¹^* and Hinrich Gronemeyer^a¹

¹Department of Cancer Biology, Institut de Génétique et de Biologie Moléculaire et Cellulaire (IGBMC)/CNRS/INSERM/Université de Strasbourg, Illkirch Cedex, France

^aDepartment of Cancer Biology, IGBMC, 1, rue Laurent Fries, BP10142, Illkirch 67404, France. Tel.: +33 3 88 65 34 73; Fax: +33 3 88 65 34 37; Email: hg/at/igbmc.u-strasbg.fr

^*Present address: Department of Plant Molecular Biology, University of Lausanne, Biophore Building, CH-1015 Lausanne, Switzerland

Received March 3, 2011; Accepted August 20, 2011.

This is an open-access article distributed under the terms of the Creative Commons Attribution Noncommercial Share Alike 3.0 Unported License, which allows readers to alter, transform, or build upon the article and then distribute the resulting work under the same or similar license to this one. The work must be attributed back to the original author and commercial use is not permitted without specific permission.

This article has been cited by other articles in PMC.

Abstract

Retinoic acid (RA) triggers physiological processes by activating heterodimeric transcription factors (TFs) comprising retinoic acid receptor (RARα, β, γ) and retinoid X receptor (RXRα, β, γ). How a single signal induces highly complex temporally controlled networks that ultimately orchestrate physiological processes is unclear. Using an RA-inducible differentiation model, we defined the temporal changes in the genome-wide binding patterns of RARγ and RXRα and correlated them with transcription regulation. Unexpectedly, both receptors displayed a highly dynamic binding, with different RXRα heterodimers targeting identical loci. Comparison of RARγ and RXRα co-binding at RA-regulated genes identified putative RXRα–RARγ target genes that were validated with subtype-selective agonists. Gene-regulatory decisions during differentiation were inferred from TF-target gene information and temporal gene expression. This analysis revealed six distinct co-expression paths of which RXRα–RARγ is associated with transcription activation, while Sox2 and Egr1 were predicted to regulate repression. Finally, RXRα–RARγ regulatory networks were reconstructed through integration of functional co-citations. Our analysis provides a dynamic view of RA signalling during cell differentiation, reveals RAR heterodimer dynamics and promiscuity, and predicts decisions that diversify the RA signal into distinct gene-regulatory programs.

Keywords: ChIP-seq, retinoic acid-induced differentiation, RXR–RAR heterodimers, temporal control of gene networks, transcriptomics

Introduction

Retinoic acid receptors (RARs) and retinoid X receptors (RXRs) are members of the nuclear receptor (NR) gene family of ligand-regulated transcription factors (TFs). RARs and RXRs form heterodimers that act as master regulators for multiple physiological processes, including embryogenesis, organogenesis, immune functions, reproduction, and organ homeostasis (Mark et al, 2006). Apart from their impact on physiology, RARs and RXRs have major promise for therapy and prevention of cancer and other diseases, and several therapeutic paradigms have been established (Altucci et al, 2007; Liby et al, 2007; Shankaranarayanan et al, 2009; de The and Chen, 2010; Zhang et al, 2010).

The biological importance of the retinoid signalling system and its cancer therapeutic potential has inspired intense research that provided detailed insight in the structural basis of, and molecular events at the early steps of retinoid action. Mechanistically, the binding of a ligand facilitates the exchange between corepressor (CoR) and co-activator (CoA) complexes by allosterically altering receptor surfaces involved in these interactions. The recruitment of such epigenetically active and/or chromatin modifying complexes leads to chromatin structure alterations and post-translational modifications that ultimately regulate cognate gene programs (Gronemeyer et al, 2004; Rosenfeld et al, 2006).

The retinoid signalling system is highly complex, as it comprises three RXR (RARα, β and γ) and three RAR (RARα, β and γ) subtypes expressed from distinct genes as multiple isoforms which act as heterodimers; in addition, RXRs can form heterodimers with a plethora of other NRs (Laudet and Gronemeyer, 2002). While insight into (some of) the physiological functions of the various RAR and RXR subtypes has been obtained by exploiting mouse genetics (Mark et al, 2006) our understanding of the cell physiological functions of these various subtypes is rather limited. The generation of subtype-selective ligands has provided important tools (de Lera et al, 2007), while the study of RAR subtype-deficient F9 embryonal carcinoma (EC) cells (Su and Gudas, 2008), despite its values, has been hampered by the observation of artifactual ligand responsiveness of the expressed RAR subtypes. Thus, we are presently facing a situation in which significant knowledge has been accumulated about the very early steps in retinoid action and the (patho)physiological impact of RAR and RXR signalling. However, what has remained entirely enigmatic is how a single compound upon activating subtype-specific RXR–RAR heterodimers can set up the temporal order of complex signalling networks that are at the basis of (patho)physiological phenomena.

Knowledge about the early events in retinoid signalling has been derived mainly from in vitro models like F9 EC cells, which differentiate into primary endodermal-like cells upon exposure to all-trans retinoic acid (ATRA); this differentiation is well characterized by morphological changes and marker expression. F9 cells display a very low rate of spontaneous differentiation, such that homogeneous cell populations can be generated during ATRA-induced differentiation. Previous studies demonstrated that, while different RXR–RAR isotype combinations control the expression of different target genes, the RXRα–RARγ heterodimer is essential for inducing differentiation (Taneja et al, 1996; Chiba et al, 1997a, 1997b). Together, these data support a model in which various RXR–RAR heterodimers regulate subtype-selective gene programs, of which RXR–RARγ establishes a path that leads to the changes which specify a differentiated F9 cell.

Here, we have addressed the question of how RXRα−RARγ upon activation by ATRA sets up a sequence of temporally controlled events that generate different subsets of primary and secondarily induced gene networks. We hypothesized that these networks required temporally defined step(s) of diversification, thereby forming separable gene cohorts that constitute the various facets of differentiation, such as altered proliferation, cell physiology, signalling, and finally terminal apoptogenic differentiation. To this aim, we performed RARγ and RXRα chromatin immunoprecipitation (ChIP) analyses coupled with massive parallel sequencing (ChIP-seq) together with the corresponding microarray transcriptomics at five time points during differentiation (Box 1). To understand the dynamics of ATRA-regulated gene expression during differentiation, gene-regulatory decisions were inferred in silico from characterized targets of RXRα−RARγ and other annotated TFs (Ernst et al, 2007). This dynamic regulatory map was used to reconstruct RXRα–RARγ signalling networks by integration of functional co-citation. Altogether, we present a genome-wide view of the temporal gene-regulatory events elicited by the RXRα–RARγ during F9 cell differentiation.

Results

Genome-wide characterization of RXRα-RARγ binding sites during ATRA-induced F9 cell differentiation

We first confirmed the induction of markers (Rarβ, Hoxa1, and Col4a1) for F9 cell differentiation by RT–PCR (Supplementary Figure S1A) and the detection of binding at previously described RAREs in the Cyp26a1 promoter (Loudig et al, 2000, 2005) using anti-RXRα antibodies (R1 and R2 in Supplementary Figure S1B and C). As expected, these sites were empty in F9 cells lacking RXRα (Rxrα^−/−).

We reasoned that combining uniquely aligned reads from all ChIP-seq time points (0, 2, 6, 24, and 48 h) would generate a valuable meta binding site profile for subsequent analyses, as it (i) cumulates all stable and transient binding events over the 48-h period and (ii) increases the peak calling confidence due to the combination of five data sets. Therefore, uniquely aligned reads from the RXRα and RARγ ChIP-seqs at different time points were combined and processed (see Materials and methods) to generate the corresponding metaprofiles.

To identify chromatin sites occupied by RXRα–RARγ heterodimers, binding sites for the two receptors in the metaprofiles were compared at different P-value thresholds and the percentage of co-occupancy was plotted for each receptor (Figure 1A). This analysis identified an optimal confidence threshold (CT40; P-value 10⁻⁴) for which all 4281 identified RARγ meta sites were co-occupied by RXRα. For the same CT RXRα bound to 9065 additional sites, most likely as heterodimer with partner(s) other than RARγ. Note that the implication of other RXRα heterodimers in ATRA-induced F9 cell differentiation has been reported (Chiba et al, 1997a).

Figure 1

RXRα and RARγ nuclear receptors present a highly dynamic binding to chromatin during ATRA-induced F9 differentiation. (A) Uniquely aligned reads sequenced from samples associated with the different time points were combined and processed (more ...)

Highly dynamic binding of RXRα–RARγ during differentiation

Temporal analysis of RXRα and RARγ at its 4281 meta binding sites revealed a highly dynamic binding (Supplementary Figure S2). In absence of ATRA, 2158 of the meta binding sites were co-occupied by RXRα and RARγ. Two hours later, 1124 additional meta sites were occupied by the heterodimer, thus increasing the number of co-occupied sites; a similar addition of new heterodimer binding sites was observed at later time points, albeit with decreasing tendency (Figure 1B). Importantly, the number of RARγ–RXRα binding sites decreased when cells moved through the differentiation program from initially ~2000 sites at 0 h to <1000 sites at 48 h. At 2 and 6 h, the gain in heterodimer binding compensated the loss of sites present at 0 h, while after 6 h there was an overall loss of RXRα–RARγ binding and at 48 h only 814 were observed. A similar loss was observed for the number of sites that were newly added at a given time point and decreased thereafter.

The observed decrease of RARγ–RXRα binding sites during differentiation could be due to (i) dissociation of both heterodimer subunits or (ii) replacement of the RXRα–RARγ by another RXR heterodimer. Monitoring the fraction of RXRα-bound sites to which RARγ is bound revealed that exposure to ATRA significantly decreased RARγ co-binding to RXRα-bound sites over time (Figure 1C). An example is the binding of the RARγ–RXRα heterodimer to the well-known RARE of the Rarβ promoter for which the level of RARγ binding decreases over time while RXRα binding is maintained, if not increased (Figure 1D). Most importantly, reChIP experiments, in which RARγ or RARα is immunoprecipitated from the RXRα ChIP, demonstrated an unexpected strong increase of RARα co-occupancy at 48 h which was not observed at earlier time points (Figure 1E and F). Note the Rarγ^−/− and Rarα^−/− F9 cell control ChIPs, which reveal the background of the assay.

Together, the above data give not only a global view of the chromatin binding dynamics of the RXRα–RARγ heterodimer but also provide moreover evidence for its replacement during F9 cell differentiation by RXRα heterodimers with other partners at common response elements. At present, we cannot distinguish between swapping of RXRα partners, i.e., dissociation followed by the formation of a distinct RXRα heterodimer, and the replacement of RXRα–RARγ by other pre-formed RXRα heterodimers.

RXRα–RARγ co-occupancy correlates with gene induction while gene repression is largely independent of this heterodimer

Transcription profiling using microarrays performed at the same time points as ChIP-seqs revealed a biphasic global gene induction with peaks at 2 and 48 h, reminiscent of results obtained by co-exposure to ATRA and cAMP (Harris and Childs, 2002). Indeed, 2 h after ATRA induction 281 genes exhibited an induction of [gt-or-equal, slanted]

1.8-fold relative to 0 h, followed by a progressive decline until 24 h (6 h, 189 genes; 24 h, 128 genes; Figure 2A). In contrast, a strong ‘wave' of gene induction was apparent at 48 h, with 926 genes getting induced.

Figure 2

Temporal correlation between RXRα–RARγ heterodimer binding and transcriptional regulation of putative target genes. (A) Genes exhibiting ATRA-induced or repressed mRNA levels at the indicated time points during F9 cell differentiation (more ...)

When comparing the differential gene expression with the location of RXRα or RARγ inferred from the metaprofiles we found that >50% of the genes induced during the first 24 h presented an RXRα or of RXRα–RARγ site within 10 kb distance (referred to as ‘putative target genes'). Similarly as for the oestrogen receptor (Carroll et al, 2005, 2006), ~70% of RXRα (heterodimer) binding sites are beyond this distance at all time points and may regulate non-annotated transcripts, such as ncRNAs, or cognate targets through chromosomal looping (Supplementary Figure S1D and E). At 48 h, the fraction of genes with RXRα/RXRα–RARγ sites dropped to 34% of all induced genes. This reveals that the majority of gene inductions at this time are due to secondary responses. Less than 10% of the downregulated genes presented a proximal RXRα or RXRα–RARγ binding site, suggesting that this heterodimer functions predominantly as positive regulator of transcription in this context.

A comparison of induced mRNA levels and gene-proximal temporal binding of RXRα–RARγ indicated a significant correlation between binding and transcription activation. Indeed, sorting of putative RXRα–RARγ target genes by induction levels revealed that at 2 h RXRα–RARγ is bound predominantly to strongly induced genes (Figure 2B). At 6 h, RXRα–RARγ binding is more prevalent at moderately induced genes, while at 24 and 48 h the number of binding events in gene-proximal RXRα–RARγ sites has dramatically decreased and the remaining subset is progressively associated with weakly induced genes.

To further assess the connection between RXRα–RARγ binding and transcription regulation of putative target genes, we mapped RNA Polymerase II (PolII) recruitment during ATRA-induced F9 cell differentiation by ChIP-seq. This analysis provided information about binding of PolII at both Transcription Start Sites (TSSs) and gene bodies (Supplementary Figure S3). For this, the PolII binding profiles were processed with POLYPHEMUS (Mendoza et al, submitted), which entails non-linear normalization of PolII enrichment of multiple ChIP-seq data sets. Genes presenting proximal binding sites for RXRα–RARγ were subsequently ranked by their PolII recruitment to TSSs at a given time point relative to 0 h. Interestingly, most of the top 50 genes (Figure 2C) presented significant PolII enrichment in both gene body and at the TSSs, indicative of active transcription. Furthermore, except Cyp26a1 and Prr14 the top 10 genes are TFs, supporting a hierarchical model of ATRA-regulated gene networks in which RXRα–RARγ induces TFs, which in turn induce their cognate gene programs.

The spatio-temporal binding of RXRα and RARγ and target gene profiling reveal distinct classes of temporally controlled gene induction patterns

To link the binding of RXRα and RARγ to transcription activation, we clustered the putative target genes by their temporal receptor binding and gene expression characteristics using a self-organization tree algorithm (SOTA; Figure 3). This classification revealed the existence of four classes of genes, which differ in the timing of heterodimer binding and gene induction (i) early induced genes with sustained expression over 48 h; (ii) early transiently induced genes; (iii) early-late transiently induced genes and (iv) late induced gene expression (Figure 3A and B). These classes contain several established RXR–RAR targets, such as Cyp26a1, Rarβ or Hoxa1 (Supplementary Table I). Note that we found a third RXRα–RARγ binding site (R3) localized ~2 kb downstream of the Cyp26a1 coding region apart from the distal (R2) and proximal (R1) RAREs and detected binding sites in genes shown to respond to ATRA but for which no RARE is described, such as Stra6, Stra8, Cdx1, Aqp3, Foxa2/HNF-3, and Nostrin/mDaIP2.

Figure 3

Temporal (transcription) regulation defines distinct classes of RXRα–RARγ target genes. (A) SOTA classification of putative RXRα–RARγ target genes according to the indicated criteria for RXRα and (more ...)

For each of the four classes the timing of coordinate binding and gene activation was the distinctive feature, while no common feature could be defined for the binding of the two receptors before or after this phase. Indeed, the ~2.5-kb distal RXRα–RARγ binding site of Aquaporin (Aqp3) (Bellemere et al, 2008; Cao et al, 2008) was co-occupied by both receptors already in absence of ATRA, while binding of RARγ was strongly reduced at 24 h and no binding of either receptor was apparent at 48 h (Figure 3C and D). In addition, co-activator components like RAC3 and p300 were recruited to this site at 2 h and were progressively reduced at later time points (Supplementary Figure S4). Notably, Aqp3 expression increased even after receptors/co-activators disappeared from the locus (Figure 3B and C; Supplementary Figure S4). As for Apq3, RXRα–RARγ occupied the putative RARE of Notch4 (Uyttendaele et al, 1998) in absence of the cognate ligand and induced transcription from 2 h on, but the loss of RARγ correlated with termination of Notch4 induction and decreasing mRNA levels. In the case of Ksr1 (Wang et al, 2006), binding of RXRα–RARγ was detected at 2–6 h, followed by a short pulse of transcriptional induction around 6 h, which ceased before 24 h together with the loss of receptors from the binding site. The late induced Nostrin (Cho et al, 1999; Cho and Park, 2000) exhibited a strongly delayed binding of RXRα and RARγ at 48 h which correlated with late RAC3 and p300 co-activator recruitment and late gene induction. The RXRα–RARγ co-occupancy of these binding sites at different time points was confirmed by reChIP assays (Figure 3E). In summary, the spatio-temporal cross-comparison between RXRα–RARγ binding and transcriptional activation revealed the existence of at least four different gene classes with distinct temporal inductions.

The putative RXRα–RARγ target genes contain a subset of promiscuously regulated genes that respond to other RAR isotypes

To assess the selectivity and promiscuity of RAR isotype signalling the use of isotype-selective ligands (de Lera et al, 2007) in the context of wild-type cells appeared to us superior to the use of RAR isotype-deficient cells, as such cells may exhibit artifactual ligand responses (Chiba et al, 1997a, 1997b). To reveal RAR isoform-selective transcription of putative RXRα–RARγ target genes, we thus used the RARγ-selective ligand BMS961. Notably BMS961, which suffices to drive F9 cells into differentiation (Taneja et al, 1996; see Supplementary Figure S5A and B), activated 62% of the ATRA-induced putative RXRα–RARγ targets (Figure 4). The RARα or RARβ-selective BMS753 and BMS641, which do not induce F9 differentiation (Taneja et al, 1996 and our unpublished results), still activated 40 and 10%, respectively, of the ATRA-induced transcriptome, thus providing evidence for both RARγ selectivity and RAR isotype promiscuity of RXRα–RARγ target genes in the context of F9 wild-type cells. That 38% of the ATRA-induced RARγ–RXRα target genes were not activated by BMS961 indicates that they are not required for F9 cell differentiation according to generally used criteria (Supplementary Figure S5). Mechanistically, these genes may be activated through direct or indirect action of RARα and/or RARβ isotypes. Possible scenarios are that both RARγ and RARα or RARβ heterodimers sequentially or coordinately bind to their regulatory regions, or that RARα or RARβ activate factors that synergize with RARγ action.

Figure 4

RXRα–RARγ putative target genes activated by specific RAR agonists. mRNA expression heatmaps of putative RXRα-RARγ target genes illustrate their induction in presence of ATRA or the indicated RAR isotype-selective (more ...)

A dynamic regulatory map for ATRA-induced F9 cell differentiation

The above results reveal that the putative RXRα–RARγ gene program suffices to trigger primitive endodermal F9 cell differentiation. It is reasonable to assume a hierarchical architecture of this program in that a few key genes coordinate cascades of gene-regulatory events thus establishing subprogram networks. Indeed, the induction of multiple TFs supports a concept in which regulatory decisions are taken, albeit not exclusively, through TF action at defined time points. To identify these decisions, we used ATRA-induced temporal gene expression, TF-target gene annotations (NCBI database annotations and/or MatInspector predictions; Cartharius et al, 2005) and the identified putative RXRα–RARγ target genes as input into the Dynamic Regulatory Events Miner (DREM; Ernst et al, 2007). DREM models bifurcation points (BPs) from the expression of a subset of genes that diverges from the co-expression pattern shared with a larger population in the previous time frame. In addition, DREM evaluates if a co-expression path is enriched for genes regulated by particular TFs whose action may account for, or contribute to the predicted bifurcation. DREM predicted six different co-expression paths from three BPs (Figure 5A). The first BP occurs between 0 and 2 h and results in the establishment of three distinct programs generating induced (orange), constitutive (grey or path (iv); this class gets induced late) and repressed (red) cohorts. The second BP subdivides the repressed path between 2 and 6 h. It separates one cohort that is progressively induced between 24 and 48 h (path (v)) from a permanently repressed gene set (path (vi)). A third BP between 6 and 24 h derives three cohorts from the induced path; one that gets repressed (path (iii)) and two others that are induced with different kinetics and mean intensities (paths (i) and (ii)). To support the validity of the predicted co-expression paths, the three gene sets originating from the first BP were classified by hierarchical clustering. As shown in Figure 5B, each of these subsets clustered into cohorts predicted by the second and third BP, with the exception of related paths (i) and (ii) which appear as one class.

Figure 5

Dynamic regulatory map of ATRA-induced transcriptome. (A) DREM co-expression analysis is represented by colour-coded paths that summarize common characteristics. The number of genes per co-expression path is indicated. Diamonds indicate three predicted (more ...)

One of the advantages of DREM is the possibility to derive associations between TFs and predicted BPs. In agreement with results described above (Figures 2A and and3),3), DREM preferentially associates RXRα–RARγ with induced paths (i) and (ii). In addition, target genes of TF-like members of the Homeobox family (e.g., Hoxa1, Hoxb2, Hoxb4, Hoxb5), Myc, Rara, Rarb, Runx1, Jun, Foxa2, Gata4, Pbx1 were also predicted to be enriched in these cohorts (see Supplementary Figure S6 for TF enrichment scores). Note that the repressed path (vi) is associated with TFs like Egr1 (Min et al, 2008) and Sox2 (Orkin et al, 2008), which are involved in regulating cell proliferation and stem cell pluripotency, respectively.

The dynamics of TF-mediated subprogramming of the RXRα–RARγ regulon is further illustrated by the temporally regulated expression of TFs themselves (Figure 5C; Supplementary Figure S6B). Indeed, with the exception of genes like Sox2, the majority of TFs are generally induced. Of interest is the biphasic response of Egr1 and Myc, which together with Sox2, is associated with class (vi) genes. Egr1 and Myc are induced when paths (v) and (vi) separate and get silenced between 6 and 24 h. This suggests that not only enhanced transcriptional activity but also temporally regulated expression of TFs contributes to the formation of temporal gene programs.

To validate the role of DREM-predicted TFs involved in BPs, we performed small interference RNA (siRNA) knockdown assays using as readout the mRNA expression of differentiation markers Laminin α1 (Lama1), Laminin β1 (Lamb1), Laminin γ1 (Lamc1), type IV collagen α1 (Col4a1); in addition, we monitored siRNA effects on the morphological changes associated with differentiation (Figure 5E–G). We also knocked down expression of Foxa1, a TF that is not predicted by DREM but is strongly and exclusively induced by ATRA and BMS961 (Figure 2C; Supplementary Figure S8; class I). Knockdown of Hoxb2, Hoxb5, Foxa1 or Foxa2 (see Supplementary Figure S7A for silencing efficiencies) reduced significantly the differentiation marker expression levels (Figure 5E). Notably, the expression levels of Nostrin, a late induced direct RXRα–RARγ target, Bmp2, an established RA target or GAPDH were not, or only marginally affected (Supplementary Figure S7B). Tracking transfection with fluorescent 6-FAM revealed that transfected cells were generally delayed (or arrested) in differentiation, while non-transfected cells within the same population exhibited a differentiated morphology (Figure 5F). Counting of blinded samples by two independent persons provided a semiquantitative analysis (Figure 5G), which fully supports the notion that these TFs have important roles in the (temporal) regulation of gene networks that are at the basis of ATRA-induced cell differentiation.

The dynamic map derived by DREM classified the differentially regulated genes during cell differentiation in six major paths, which can be distinguished by the relative enrichment of their components according to Gene ontology (GO) terms (Figure 5D; Supplementary Figure S8). Indeed, while the early and sustained induced paths (i) and (ii) are enriched for genes related to embryonic morphogenesis and actin cytoskeleton organization, respectively, the early temporally induced path (iii) is enriched for genes involved in steroid/cholesterol metabolic processes. The late induced path (iv) is associated with cell adhesion, positive regulation in response to external stimuli while path (v) is linked to cell-cycle regulation. Interestingly, the repressive path (vi) is enriched for genes that negatively regulate cell differentiation.

A comprehensive ATRA-induced RXRα–RARγ signalling network

With the aim of enhancing the dynamic landscape of the RXRα–RARγ regulome inferred by DREM, we reconstructed the corresponding gene networks on the basis of functional co-citation (Genomatix Bibiosphere PathwayEdition) and the identification of essential nodes by topology-based scoring methods (cytoHubba; Lin et al, 2008). The illustration of the resulting RXRα–RARγ regulome (Figure 6; Supplementary File S1) depicts the relevant components of the six co-expression classes (compare Figure 5) and specifies their intraclass and interclass co-citation interactions.

Figure 6

A comprehensive ATRA-RXRα/RARγ signalling network. Genes associated with the different co-expression paths illustrated in Figure 5 are represented in the context of their functional gene co-citation interactions. For simplicity, only the (more ...)

Several general features can be extracted from this dynamic network of co-expression classes. First, each class is unique in expressing a particular set of genes with similar general functionality, such as the TF-rich class (i). Second, genes regulating complex biological phenomena may appear in different classes with distinct expression profiles, as the subsequent inductions of cyclins and cyclin-dependent kinase inhibitors. Third, the present ChIP-seq data identify putative RAREs in a great number of genes, some of which are known to respond to retinoids (see Supplementary Table I). Fourth, the described F9 RXRα–RARγ regulome integrates several factors with important roles in other cell systems, such as Egr1 and Notch4. Fifth, comparing regulation of the putative target genes by subtype-selective ligands reveals RAR subtype selectivity and promiscuity; moreover, the subset of genes commonly regulated by ATRA and BMS961 which are divergently regulated by RARα and RARβ ligands is likely constitute the bona fide differentiation program.

Within class (i), topology-based scoring identified Jun, Myc, Rara or Rarb as most important nodes. While the positive regulation of Jun and Myc expression by ATRA has been described (Supplementary Table I) the biphasic expression seen upon ATRA exposure is not maintained with the RARγ-selective BMS961 (Supplementary Figure S6B). Indeed, BMS961 only recapitulates the early and late downregulation of the expression of Jun and Myc, respectively. Thus, the temporally regulated repression but not the induced expression of these TFs correlates with cell differentiation. Multiple other TFs contribute to the definition of class (i). Apart from two other RAR isotypes, there is a strong representation of members of the homeobox TF family, including Cdx1, Meis2, some of which have well-characterized RAREs (Supplementary Table I) and served as validation marks for our ChIP-seqs. Finally, Foxa1 and two NR co-regulators (Ncoa7 and Nrip1) are putative regulatory factors of class (i). In addition to TFs, this class contains also RA-target genes involved in retinoid homeostasis, including Cyp26a1, Crabp2 or Rbp1. Importantly, all of these genes are similarly regulated by ATRA and BMS961 but not by BMS753 or BMS641 (Supplementary Figure S9), thus supporting a functional role in F9 cell differentiation.

According to GO terms, class (i) is predicted to trigger positive regulation of transcription, cell differentiation and responses to vitamin A. Class (ii), which shares a common ancestor with classes (i) and (iii), is characterized by the enrichment of genes involved in actin cytoskeleton organization (Supplementary Figure S8). This cohort contains also several apoptogenic factors, including Casp3, Casp8, Bcl2l11 and Mcl1, and the signalling factors Jak2, Rhob and Pim; several of these genes are known to respond to retinoids (Supplementary Table I). Comparing the induction profiles of these genes by the three RAR subtype-selective agonists indicates that their ATRA regulation may not be directly linked to F9 differentiation; examples for this notion are Id2, Casp 3 or Pim1 (see class (ii) in Supplementary Figure S9).

Several genes that are components of a similar biological process are found in different classes and it is tempting to speculate that this may be linked to their distinct temporal role during the differentiation process. For instance, the temporally controlled expression of cell cycle-regulatory genes during differentiation can be rationalized from their functionalities; while Ccnd2 (class (ii)) expression increases, the class (iii) genes Ccnd1, Ccnd3 and Ccne1 are only transiently induced and regress with the expression of ‘late' class (v) genes p27/Cdkn1b, p21/Cdkn1a and p57/Cdkn1c, which were repressed at earlier time points. Note that these observations corroborate previous reports showing that also the protein levels of cyclin D2 and p27/Cdkn1b increased during F9 cell differentiation, whereas that of cyclin D1, D3 and cyclin E decreased (Li et al, 2004).

In addition to cell cycle-regulatory components, class (iii) contains also genes like Egr1 and Notch4. The ATRA-responsive Egr1 (Edwards et al, 1991) is sufficient to inhibit cell proliferation of haematopoietic stem cells (Min et al, 2008) where Notch4 negatively regulates cell differentiation (Vercauteren and Sutherland, 2004; Ye et al, 2004). Note the temporal binding of RXRα–RARγ to the Notch4 promoter region in F9 cells, which correlates with its transient mRNA expression (Figure 3).

The ‘late' class (iv) comprises genes involved in processes like extracellular matrix organization (Col18a1, Sparc, Timp3) or cell adhesion (Cd9, Cd47, Nedd9, Cdh2, Thbs1, Sirpa, Cdh5, Cyr61, Cdh11). These gene regulations are likely readouts of the morphological changes associated with differentiation, supported by the fact that RARα or RARβ-specific agonists do not induce a similar expression as ATRA or RARγ-specific agonists (Supplementary Figure S9). Also, late induced TFs like Foxa2/HNF-3, Runx1 and the TF-related factor Zmiz1 are members of this class (Supplementary Table I). Finally, germline-specific differentiation markers, like Otx2 and Fgf4 for ectodermal or Bmp2, Gata4 and Hnf1b for endodermal differentiation were retrieved in classes (iv) and (vi) for which their temporal transcription pattern correlates with differentiation into primitive endoderm (Supplementary Figure S9).

Taken together, the comprehensive co-expression network, reconstructed from an integrative analysis of RXRα–RARγ binding and transcription regulation, illustrates the complex temporal coordination of a plethora of molecular processes during differentiation, including cell-cycle regulation, the activation of apoptotic/cell survival or repression of self-renewal/pluripotency.

Discussion

NR-regulated cell physiological phenomena, such as differentiation, growth, survival or death, are at the basis of complex physiological processes. Despite an enormous gain of knowledge about molecular and structural features of NR function (Perissi and Rosenfeld, 2005; Altucci et al, 2007; O'Malley and Kumar, 2009), we are still far from understanding of how a single compound, such ATRA, can induce temporal patterns of coordinate regulation of gene networks that finally lead to the observed changes of cell, organ or animal physiology. With the availability of a number of genome-wide profiling technologies, it is now feasible to initiate integrative analyses to decipher the gene-regulatory phenomena and describe the dynamic regulation of gene networks. To this aim, we have used the ATRA-induced F9 EC cell differentiation model and describe to our knowledge the first dynamic analysis of the RXRα–RARγ regulon by an integrative analysis of genome-wide temporal binding of the RXRα–RARγ heterodimer, the corresponding temporal gene regulation patterns and the response to subtype-selective RAR agonists. Our aim was to understand signal diversification and specification at different functional levels, ranging from the specific action of RAR subtypes to decisions within the process of differentiation that result in the formation of separate gene programs and finally to the description of the key components therein which may provide clues to the program function(s).

One of the unexpected outcomes of this analysis is the astounding dynamics of RXRα–RARγ binding. It was recently shown that ATRA-induced differentiation of mouse ES cells to motor neurons results in widespread changes in RAR genomic binding (Mahony et al, 2011). This documented that in contrast to the apo oestrogen receptor, which binds to a limited number of target sites and acquires extensive chromatin binding capacity in presence of agonists (Carroll et al, 2005, 2006; Ceschin et al, 2011), apoRAR already binds to a large number of sites (some of) which it may silence by recruiting corepressor complexes. On the other hand, ligand exposure was shown to generate de novo RAR binding to certain target genes (Chen et al, 1996). Our present data not only confirm those studies at the genome-wide level but also demonstrate, moreover, a precise temporal order of gain and loss of binding sites for a particular RXR–RAR heterodimer during the first 48 h of differentiation. We observed an unexpected fluctuation of different RXRα heterodimers at sites initially occupied by RXRα–RARγ. Two interesting scenarios could account for this phenomenon, replacement of RXRα–RARγ by another RXRα heterodimer or swapping of RXRα partners. This temporal order of binding of different RAR heterodimers to the same site is supported by the observations (i) that reChIP experiments reveal a temporally distinct co-binding of RARγ and RARα together with RXRα to the Rarb RARE (Figure 1E and F) and (ii) that a significant number of putative RXRα–RARγ target genes respond not only to RARγ-selective but also to RARβ and/or RARα-selective ligands (Figure 4; Supplementary Figure S9). Taken together, our analysis of RXR–RAR binding reveals a highly dynamic occupancy of its target loci suggesting the existence of mechanism(s) that orchestrate in a clock-like manner the sequential recruitment and release of RXR heterodimers. It is tempting to speculate that site-selective chromatin modifications have a role in this process, which may involve previously observed histone methyltransferase and demethylase gatekeepers (Garcia-Bassets et al, 2007).

Such mechanism(s) may also account for the absence of temporal correlation between heterodimer binding and gene activation. Indeed, the presence of RXRα–RARγ localization preceded in certain cases by several hours the transcription induction, indicating that holoRXRα–RARγ binding is not per se sufficient for transcription induction; inversely, the dissociation of RXRα–RARγ does not necessarily correlate with attenuated transcription. Due to this lack of temporal correlation the generation of metaprofiles from different time points facilitated the comparison with temporal gene regulation.

The observation that distinct RAR subtype heterodimers can bind at identical target sites complicates the definition of subtype-selective target gene repertoires. As it has been shown that gene deletion can generate artifactual RAR subtype responses (Chiba et al, 1997a, 1997b), and affect global gene expression profiling even in absence of any treatment (Su and Gudas, 2008), we used RAR subtype-selective agonists (Gehin et al, 1999; Germain et al, 2004; de Lera et al, 2007) to define selective and promiscuous functionalities of RAR subtypes. Notably, the RARγ-specific agonist BMS961, which is competent to induce F9 cell differentiation, activated >60% of the putative RXRα–RARγ target genes, thus confirming a causal link between binding and transcription activation for this cohort. While our analysis identified the RARγ-selective gene program for F9 cell differentiation and revealed RAR subtype promiscuity, definition of the complete program will require the target site identification of the great number of binding sites that could not be linked to annotated genes.

To approach an understanding of the decisions that lead to the establishment of subprograms, we used DREM for the in silico reconstruction of a dynamic regulatory map. This approach was supported by the observation that multiple TF were among the genes regulated by ATRA and BMS961. DREM analysis predicted six different gene co-expression paths, which result from three distinct bifurcations during the first 6 h of ATRA treatment. Note that this is the time span required for differentiation commitment (Levine et al, 1984). Some of the TFs associated with the DREM-predicted gene co-expression paths were shown here to have an important role during F9 cell differentiation (Figure 5E–G; Supplementary Figure S7), others were described previously (Gata4, Futaki et al, 2004; Sox2, Chew et al, 2005; Boer et al, 2007).

To gain insight into the functionality of the co-expression gene cohorts and identify key players, we established a dynamic regulatory map and constructed a differentiation-related gene network from functional co-citation. While not exhaustive, this RXRα–RARγ driven dynamic network provides a global view on the regulatory events during ATRA-induced F9 cell differentiation. Indeed, genes associated with each class were shown to share similar general functions and genes of distinct classes displayed a coordinate temporal order of expression that is requisite for regulating complex biological phenomena. That an important number of RXRα–RARγ target genes did not show up in the reconstructed network is due to their original identification in the present study, accounting for the lack of co-citation. It is worth pointing out that a comparison of gene regulation by subtype-selective ligands, as illustrated in the reconstructed network, was useful to validate their importance for differentiation and discover RAR subtype promiscuity.

Overall, the current study illustrates for the F9 model the different levels of ATRA-induced signalling pathway diversification due to regulatory decisions at different levels and time points. Whereas other regulatory principles, like the chromatin modification status and co-regulator function and modification (Ceschin et al, 2011) may be overlaid at any of these levels, the current gene network identifies potential nodes that act as key regulators of various subprograms of RA-initiated signal transduction during differentiation. It will be interesting to integrate other RXR–RAR components and compare other RA-regulated cell models to retrieve common and cell-specific features. Ultimately, it may be possible to predict critical nodes from applying computational models to reconstructed gene networks (Fisher and Piterman, 2010) and extrapolate this information towards other model systems that mimic RAR-dysfunctional human diseases (Gronemeyer and Zelent, 2009; Mikesch et al, 2010), as targetable key factors for therapy.

Materials and methods

F9 mouse EC cell culture

EC cells were cultured in Dulbecco's modified Eagle's medium supplemented with 10% fetal calf serum and 40 μg/ml Gentamicin. Cells were seeded in gelatine-coated tissue culture plates (0.1%) and ATRA was added to the plates in a final concentration of 1 × 10⁻⁶ M at different time points. For assays involving RAR subtype-specific agonists, cells were incubated with BMS961 (RARγ specific; final concentration 10⁻⁷ M), BMS753 (RARα specific; final concentration 10⁻⁶ M) or BMS641 (RARβ specific; final concentration 10⁻⁶ M).

ChIP and reChIP assays

Cells were fixed with 1% para-formaldehyde (Electron Microscopy Sciences) for 30 min at room temperature. ChIP assays were performed following standard conditions: chromatin sonication and immunoprecipitation in lysis buffer (50 mM Tris–Cl pH=8, 140 mM NaCl, 1 mM EDTA, 1% Triton, 0.1% Na-deoxycholate) complemented with protease inhibitor cocktail (Roche 11873580001); 2 × washes with lysis buffer; 2 × washes with lysis buffer containing 360 mM NaCl; 2 × washes with washing buffer (10 mM Tris–Cl pH=8, 250 mM LiCl, 0.5% NP-40, 1 mM EDTA, 0.5% Na-deoxycholate); 2 × washes with 1 × TE; elution at 65°C; 15 min in elution buffer (50 mM Tris–Cl pH=8, 10 mM EDTA, 1% SDS). RXRα and RARγ have been immunoprecipitated with antibodies generated by immunization of rabbits with the following peptides:

mRXRα: PB105 (MDTKHFLPLDFSTQVNSSSLNSPTGRGC),

mRARγ: PB288 (CSKPGPHPKASSEDEAPGGQGKRGQSPQPD).

Polyclonal anti-RXRα and anti-RARγ were purified from the crude serum by affinity chromatography. RNA Polymerase II (sc-9001 H-224), p300 (sc-584 N-15) and Rac3/NCoA-3 (sc-9119) antibodies were purchased from Santa Cruz biotechnology (Santa Cruz, CA). For RXRα, RARγ, p300 and RAC3 6 × 10⁶ cells were used per ChIP, whereas for RNA polII 2 × 10⁶ cells were used. For reChIPs, at least four ChIP assays of 6 × 10⁶ cells were used for the first IP. For reChIPs, the first antibody (anti-RXRα) was covalently linked to the sepharose protein A (Sigma P92424) using disuccinimidyl suberate (DSS). The covalently linked Ab beads were washed with ethanolamine (0.1 M), followed by glycin at pH=2.8 to remove non-covalently linked antibodies from the beads. Beads were washed with 50 mM sodium borate at pH=8.2 and PBS, and were incubated overnight at 4°C with the corresponding whole cell extract as in a regular ChIP assay. Following standard washing, elution was performed with 10 mM DTT (30 min, 37°C). The eluates from four ChIPs were combined, diluted at least 30 times with lysis buffer (containing protease inhibitors like in a regular ChIP assay), and incubated overnight with the second antibody (anti-RARγ) and sepharose Protein A beads at 4°C. The subsequent steps were performed as for regular ChIPs.

The immunoprecipitated chromatin was validated using positive (recruitment to known targets) and negative (‘cold' region) controls and the binding was expressed as enrichment relative to the whole cell extract input control (% input) and/or relative to a ‘cold' region (fold occupancy); validation assays were performed by quantitative real-time PCR (qPCR, Roche LC480 light cycler device) using Quantitect (Qiagen).

Massive parallel sequencing

After qPCR validation, immunoprecipitated chromatin was quantified using Qubit (Quant-It dsDNA HS Assay Kit; Invitrogen). In all, 10 ng of the ChIPed material was used for preparing the sequencing library (ChIP-seq DNA sample preparation kit, Illumina). In all, 5 pmol of the library was used per flow cell in the Solexa 2G Genome analyzer (Illumina). The Illumina Pipeline v1.4.0 was used for primary data analysis (image processing: Firecrest; base calling: Bustard; alignment: Gerald) from which uniquely aligned reads with up to two mismatches relative to the mouse mm9 reference genome were kept.

Peak detection approach

MACS v. 1.3.6 (Zhang et al, 2008) was used as peak caller at the first level of data treatment to obtain signal intensity wiggle files and annotated peak regions using a Poisson model distribution (P-value confidence cutoff: 1 × 10⁻⁵). Considering that a certain number of false positive regions were observed during this treatment, which in our hands cannot be decreases without compromising true positives by playing with the P-value confidence parameter, we used subsequently an approach based on a machine-learning algorithm able to define the model distribution associated with the ChIP-seq data set under study. This approach, implemented in the R package MeDiChI, has been previously described (Reiss et al, 2008) for analysis of ChIP-chip assays; its implementation for ChIP-seq assays will be described elsewhere (manuscript in preparation).

RT–PCR and microarrays

Total RNA was extracted from cells, treated with ATRA during different period of times, using the GENEElute^TM RNA extraction kit (Sigma). In all, 2 μg of the eluted RNA was used for reverse transcription (AMV-RTase, Roche; oligo(dT) New England Biolabs; 1 h; 42°C). The cDNA was diluted 10-fold and used for real-time qPCR (Roche LC480). Expression of the following marker genes was assessed to follow the process of F9 cell differentiation:

For early response to the differentiation inductor (ATRA):

Hoxa1 (CCCAGACGGCTACTTACCAG; CATGGGAGTCGAGAGGTTTC);

Rarb (GATCCTGGATTTCTACACCG; CACTGACGCCATAGTGGTA).

For late response to the differentiation inductor (ATRA):

Col4a1 (ATGCCCTTTCTCTTCTGCAA; ATCCACAGTGAGGACCAACC);

Lama1 (CCGACAACCTCCTCTTCTACC; TCTCCACTGCGAGAAAGTCA);

Lamb1 (TCTATGCTCGGCAGTGTGAC; CAGTGGTCTCCTGACCCAAT);

Lamc1 (GGCCGAGTGCCTACAACTT; CAGTGGCAGTTACCCATTCC).

During data analysis, all qPCR values were normalized relative to the constitutively expressed 36B4 mRNA (AATCTCCAGAGGCACCATTG; CCGATCTGCAGACACACACT).

Using 250 ng of starting total RNA, biotin-labelled cDNAs were synthesized and hybridized on Affymetrix GeneChip® Mouse Gene 1.0 ST Array (Affymetrix, Santa Clara, CA, USA) according to the manufacturer's recommendations as described in ‘GeneChip® whole transcript sense target labelling assay manual' (P/N 701880 Rev.4). The arrays were further washed and stained with streptavidine-phycoerythrin in an Affymetrix GeneChip® Fluidics station 450 using the script protocol F450-0007, then scanned with an Affymetrix Gene Chip Scanner 3000 7G. Expression values were generated with Affymetrix software Expression Console version 1.1, using sketch quantile normalization and median polish summarization as in Robust Multiarray Analysis.

RNA polymerase II ChIP-seq data processing

RNA Polymerase II genome-wide intensity profiles were generated using MACS v. 1.3.6 (Zhang et al, 2008) as peak caller and signal intensity profiles were processed by MeDiChI to identify chromatin regions presenting significant RNA Polymerase II enrichment (P-value cutoff: 1 × 10⁻²). MeDiChI-predicted enriched regions were further filtered according to their genomic localization by comparing them with the mm9 coding region Ref-seq annotation. Thereby coding regions presenting significant levels of RNA PolII at the TSS were identified and selected for further analysis. A Quantile normalization procedure was used to enable a comparison between different samples. Normalized profiles were compared in the context of their signal intensities to identify local changes in the Pol II occupancy relative to the control sample. Note that this approach differs from previously described linear corrections between samples based on total numbers of reads; it has been implemented in the R package POLYPHEMUS (Mendoza et al, submitted). The final comparisons between the PolII levels at the TSS relative to the average behaviour through the gene body have been used as readouts to identify and describe genes presenting an enhancement of their transcriptional activity at the different analysed time points (Figure 2).

Data integration

The global RXRα and RARγ localization at different time points during ATRA-induced differentiation has been followed by ChIP-seq. To increase the confidence in binding localization assessment, we collected the mappable reads from all five time points in a single file called metaprofile. The metaprofile associated with the localization of RXRα and RARγ was processed by MACS and MeDiChI to identify significant enriched regions. Interestingly, this approach generates a higher signal/background contrast, thus increasing the sensitivity of peak calling. Taking in consideration the heterodimeric nature between RXRα and RARγ, we compared the two metaprofiles to identify the optimal confidence parameter at which the coexistence between RXRα and RARγ is maximal (see Figure 1A). Indeed, for a P-value threshold of 1 × 10⁻⁴ (CT40; where CT=−10 × log (P-value)) all identified RARγ sites in the metaprofile colocalize with RXRα. These co-occupied sites were further evaluated in the context of their temporal RXRα and RARγ co-occupancy. Whereas a CT40 has been used for metaprofiles binding sites selection, the same parameters are too stringent for the analysis of time-point profiles. We have taken a CT25 per time-point profile for evaluating the RXRα and/or RARγ occupancy at the pre-identified sites (Figure 1B).

RXRα–RARγ co-occupied sites (metaprofiles comparison) were annotated based on their proximity (<10 kb) to transcriptionally active coding regions (upregulation and downregulation levels relative to control sample: ratio cutoff [gt-or-equal, slanted]

1.8 and

0.5, respectively); annotated genes are referred to as ‘putative' target genes. In order to gain temporal information about the correlation between RXRα–RARγ binding and putative target gene transcription, the RXRα and RARγ binding at a given time point was compared with the temporal mRNA levels of putative target genes. To further remove possible false positive annotations, genes presenting coexistence between RXRα–RARγ binding and differential mRNA expression levels at least at one time point were retained during selection. For classification purposes, the coexistence of several events was scored per time point following a hierarchical order of importance in the context of gene regulation:

An external file that holds a picture, illustration, etc.
Object name is msb201173-m2.jpg

Finally, the temporal incidence of the evaluated events was classified using an SOTA approach (Euclidean distance, Max. cycles=7, cell variability P-value=0.01) under the open access multiExperiment Viewer (Saeed et al, 2003, 2006; see Figure 3).

Dynamic regulatory maps and Network reconstruction

RXRα–RARγ putative target gene information as well as further TF-gene regulatory interaction annotations extracted from the NCBI database and/or predicted by MatInspector (Cartharius et al, 2005) were combined into a binary matrix (‘1' for TF-target gene association; otherwise ‘0'). Furthermore, time-point transcriptomics data were expressed in log2 ratios relative to the 0-h control. These two data sets were integrated into the DREM (Ernst et al, 2007) to establish dynamic regulatory maps. In addition to classifying temporal sets of gene expression data into co-expression paths, DREM predicts BPs, defined as a time point at which the co-expression behaviour of a subset of genes diverges from the common path present in a previous time point. Moreover, the association of a given TF with the corresponding BPs is estimated by evaluating the enrichment of its target genes using hypergeometric distributions relative to the genes associated with the common path.

To construct gene networks, genes associated with each predicted co-expression path were evaluated in the context of their functional co-citation relationships (Genomatix Bibiosphere PathwayEdition). To identify functional co-citation interactions between different co-expression paths, all genes used for the initial analysis were evaluated for their functional co-citations relationships as explained before. Finally, the intra- and inter-co-expression paths co-citation interactions were integrated into Cytoscape (Shannon et al, 2003; Cline et al, 2007). To further increase the confidence of the reconstructed network, topological information concerning the number of edges per node (Hubba; Lin et al, 2008) were used as filter in the final representation. Briefly, we have used Hubba's Double screening scheme (DSS) approach, which combines a Maximum Neighbourhood Component (MNC) with a Density of Maximum Neighbourhood Component (DMNC) to score nodes based on the number of their interconnections. The final gene-network representation corresponds to the top 100 ranked nodes (red to green colour code) as well as their corresponding first level neighbours. Additional information, like TF annotations, BMS961 and ATRA responsiveness, number of functional co-citations between nodes (edge's broadness) and co-expression path belonging (colour coded) are also included in the final display of this RXRα–RARγ induced genes network, which is also available in Cytoscape format (Supplementary File S1). Transcriptomics data associated with ATRA and RAR subtype-specific agonist treatments were also included as attributes for the reconstructed network (Supplementary Files S2–S6).

siRNA transfections

F9 cells were transfected with siRNA oligomers directed against RNA target sequences for the following TFs: Hoxb2 (QIAGEN; FlexiTube GeneSolution: Cat. No. SI01069089, SI01069096, SI01069075, SI01069082; working concentration: 50 nM); Hoxb5 (QIAGEN; FlexiTube GeneSolution: Cat. No. SI01069173, SI01069180, SI01069159, SI01069166; working concentration: 50 nM); Foxa1 (QIAGEN; FlexiTube GeneSolution: Cat. No. SI01004493, SI01004500, SI01004479, SI01004486; working concentration: 50 nM); Foxa2 (QIAGEN; FlexiTube GeneSolution: Cat. No. SI01004528, SI02737182, SI01004514, SI01004521; working concentration: 50 nM); Gata4 (QIAGEN; FlexiTube GeneSolution: Cat. No. SI01009813, SI01009820, SI01009799, SI01009806; working concentration: 50 nM). In addition, F9 cells were transfected with an siRNA oligomer directed against the RNA target sequence of GFP (Dharmacon; P-002048-01-20; working concentration: 50 nM) as a control for the specificity of the assay. The transfection efficiency has been followed by co-transfecting the previous oligomers with the 6-FAM labelled siGLO transfection indicator (Dharmacon; D-001630-01-05; working concentration: 50 nM). F9 cells were transfected with Lipofectamine 2000 (Invitrogen; Cat. No. 11668-027) during 18 h, then medium has been changed and cells were treated either with ATRA or with ethanol during 48 h as previously described.

Data access

Affymetrix microarrays as well as Illumina platform ChIP-seq data described in this study are available from the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo) under accession number GSE30539.

Supplementary Material

Supplementary File S1

Cytoscape file for ATRA-induced RXRα-RARγ functional co-citation network.

Click here to view.^{(247K, zip)}

Supplementary File S2

Cytoscape file for ATRA-induced RXRα-RARγ functional co-citation network displaying temporal transcriptomics information per node associated to ATRA and RAR-specific agonists treatment.

Click here to view.^{(511K, zip)}

Supplementary File S3

RXRα-RARγ functional co-citation network illustrating the temporal transcriptomics information per node associated to ATRA treatment

Click here to view.^{(20M, ppt)}

Supplementary File S4

RXRα-RARγ functional co-citation network illustrating the temporal transcriptomics information per node associated to RARγ-specific agonist (BMS961) treatment.

Click here to view.^{(12M, ppt)}

Supplementary File S5

RXRα-RARγ functional co-citation network illustrating the temporal transcriptomics information per node associated to RARα-specific agonist (BMS753) treatment.

Click here to view.^{(12M, ppt)}

Supplementary File S6

RXRα-RARγ functional co-citation network illustrating the temporal transcriptomics information per node associated to RARβ-specific agonist (BMS641) treatment.

Click here to view.^{(11M, pps)}

Suppl Information

Click here to view.^{(1.9M, pdf)}

Review Process File

Click here to view.^{(220K, pdf)}

Acknowledgments

We would like to thank Bernard Jost and Serge Vicaire for sequencing library preparation and Solexa sequencing; Stephanie Le Gras for computational Illumina pipeline treatment; Ingrid Colas and Christelle Thibault for Microarray gene expression samples preparation and primary data processing, all the members of the IGBMC Microarray and Sequencing Platform and all the members of our laboratory for discussions. MW is supported by a fellowship of the Association de Recherche Contre le Cancer and the Fondation pour la Recherche Médicale. This work was supported by the European Community (contracts HEALTH-F4-2007-200767 ‘APO-SYS', LSHC-CT-2005-518417 ‘EPITRON', HEALTH-F4-2009-221952 ‘ATLAS' and LSHG-CT-2005-018882 ‘X-TRA-NET') and the Ligue National Contre le Cancer (laboratoire labelisé).

Author contributions: MAMP produced the ChIP-seq and transcriptomics data. In addition, he developed/implemented the in silico approaches like the use of metaprofiles, MeDiChI-seq; ChIP-seq/transcriptomics data clustering. RNA Polymerase II ChIP-seq data normalization was developed by MAMP and MS. MW generated most of the ChIP–qPCR and RT–qPCR data and collaborated with MAMP for ChIP. HG planned, coordinated and supervised the project. MAMP and HG wrote the manuscript. All authors read and approved the manuscript.

Footnotes

The authors declare that they have no conflict of interest.

References

Altucci L, Leibowitz MD, Ogilvie KM, de Lera AR, Gronemeyer H (2007) RAR and RXR modulation in cancer and metabolic disease. Nat Rev Drug Discov 6: 793–810. [PubMed]
Ashburner M, Chihara C, Meltzer P, Richards G (1974) Temporal control of puffing activity in polytene chromosomes. Cold Spring Harb Symp Quant Biol 38: 655–662. [PubMed]
Bellemere G, Von Stetten O, Oddos T (2008) Retinoic acid increases aquaporin 3 expression in normal human skin. J Invest Dermatol 128: 542–548. [PubMed]
Boer B, Kopp J, Mallanna S, Desler M, Chakravarthy H, Wilder PJ, Bernadt C, Rizzino A (2007) Elevating the levels of Sox2 in embryonal carcinoma cells and embryonic stem cells inhibits the expression of Sox2:Oct-3/4 target genes. Nucleic Acids Res 35: 1773–1786. [PMC free article] [PubMed]
Cao C, Wan S, Jiang Q, Amaral A, Lu S, Hu G, Bi Z, Kouttab N, Chu W, Wan Y (2008) All-trans retinoic acid attenuates ultraviolet radiation-induced down-regulation of aquaporin-3 and water permeability in human keratinocytes. J Cell Physiol 215: 506–516. [PubMed]
Carroll JS, Liu XS, Brodsky AS, Li W, Meyer CA, Szary AJ, Eeckhoute J, Shao W, Hestermann EV, Geistlinger TR, Fox EA, Silver PA, Brown M (2005) Chromosome-wide mapping of estrogen receptor binding reveals long-range regulation requiring the forkhead protein FoxA1. Cell 122: 33–43. [PubMed]
Carroll JS, Meyer CA, Song J, Li W, Geistlinger TR, Eeckhoute J, Brodsky AS, Keeton EK, Fertuck KC, Hall GF, Wang Q, Bekiranov S, Sementchenko V, Fox EA, Silver PA, Gingeras TR, Liu XS, Brown M (2006) Genome-wide analysis of estrogen receptor binding sites. Nat Genet 38: 1289–1297. [PubMed]
Cartharius K, Frech K, Grote K, Klocke B, Haltmeier M, Klingenhoff A, Frisch M, Bayerlein M, Werner T (2005) MatInspector and beyond: promoter analysis based on transcription factor binding sites. Bioinformatics 21: 2933–2942. [PubMed]
Ceschin DG, Walia M, Wenk SS, Duboe C, Gaudon C, Xiao Y, Fauquier L, Sankar M, Vandel L, Gronemeyer H (2011) Methylation specifies distinct estrogen-induced binding site repertoires of CBP to chromatin. Genes Dev 25: 1132–1146. [PubMed]
Chen JY, Clifford J, Zusi C, Starrett J, Tortolani D, Ostrowski J, Reczek PR, Chambon P, Gronemeyer H (1996) Two distinct actions of retinoid-receptor ligands. Nature 382: 819–822. [PubMed]
Chew JL, Loh YH, Zhang W, Chen X, Tam WL, Yeap LS, Li P, Ang YS, Lim B, Robson P, Ng HH (2005) Reciprocal transcriptional regulation of Pou5f1 and Sox2 via the Oct4/Sox2 complex in embryonic stem cells. Mol Cell Biol 25: 6031–6046. [PMC free article] [PubMed]
Chiba H, Clifford J, Metzger D, Chambon P (1997a) Distinct retinoid X receptor-retinoic acid receptor heterodimers are differentially involved in the control of expression of retinoid target genes in F9 embryonal carcinoma cells. Mol Cell Biol 17: 3013–3020. [PMC free article] [PubMed]
Chiba H, Clifford J, Metzger D, Chambon P (1997b) Specific and redundant functions of retinoid X Receptor/Retinoic acid receptor heterodimers in differentiation, proliferation, and apoptosis of F9 embryonal carcinoma cells. J Cell Biol 139: 735–747. [PMC free article] [PubMed]
Cho SY, Lee SH, Park SS (1999) Differential expression of mouse Disabled 2 gene in retinoic acid-treated F9 embryonal carcinoma cells and early mouse embryos. Mol Cells 9: 179–184. [PubMed]
Cho SY, Park SS (2000) Genomic organization and promoter analysis of mouse disabled 2 gene. Biochem Biophys Res Commun 275: 189–194. [PubMed]
Cline MS, Smoot M, Cerami E, Kuchinsky A, Landys N, Workman C, Christmas R, Avila-Campilo I, Creech M, Gross B, Hanspers K, Isserlin R, Kelley R, Killcoyne S, Lotia S, Maere S, Morris J, Ono K, Pavlovic V, Pico AR et al. (2007) Integration of biological networks and gene expression data using Cytoscape. Nat Protoc 2: 2366–2382. [PubMed]
de Lera AR, Bourguet W, Altucci L, Gronemeyer H (2007) Design of selective nuclear receptor modulators: RAR and RXR as a case study. Nat Rev Drug Discov 6: 811–820. [PubMed]
de The H, Chen Z (2010) Acute promyelocytic leukaemia: novel insights into the mechanisms of cure. Nat Rev Cancer 10: 775–783. [PubMed]
Edwards SA, Darland T, Sosnowski R, Samuels M, Adamson ED (1991) The transcription factor, Egr-1, is rapidly modulated in response to retinoic acid in P19 embryonal carcinoma cells. Dev Biol 148: 165–173. [PubMed]
Ernst J, Vainas O, Harbison CT, Simon I, Bar-Joseph Z (2007) Reconstructing dynamic regulatory maps. Mol Syst Biol 3: 74. [PMC free article] [PubMed]
Fisher J, Piterman N (2010) The executable pathway to biological networks. Brief Funct Genomics 9: 79–92. [PMC free article] [PubMed]
Futaki S, Hayashi Y, Emoto T, Weber CN, Sekiguchi K (2004) Sox7 plays crucial roles in parietal endoderm differentiation in F9 embryonal carcinoma cells through regulating Gata-4 and Gata-6 expression. Mol Cell Biol 24: 10492–10503. [PMC free article] [PubMed]
Garcia-Bassets I, Kwon YS, Telese F, Prefontaine GG, Hutt KR, Cheng CS, Ju BG, Ohgi KA, Wang J, Escoubet-Lozach L, Rose DW, Glass CK, Fu XD, Rosenfeld MG (2007) Histone methylation-dependent mechanisms impose ligand dependency for gene activation by nuclear receptors. Cell 128: 505–518. [PMC free article] [PubMed]
Gehin M, Vivat V, Wurtz JM, Losson R, Chambon P, Moras D, Gronemeyer H (1999) Structural basis for engineering of retinoic acid receptor isotype-selective agonists and antagonists. Chem Biol 6: 519–529. [PubMed]
Germain P, Kammerer S, Perez E, Peluso-Iltis C, Tortolani D, Zusi FC, Starrett J, Lapointe P, Daris JP, Marinier A, de Lera AR, Rochel N, Gronemeyer H (2004) Rational design of RAR-selective ligands revealed by RARbeta crystal structure. EMBO Rep 5: 877–882. [PubMed]
Gronemeyer H, Gustafsson JA, Laudet V (2004) Principles for modulation of the nuclear receptor superfamily. Nat Rev Drug Discov 3: 950–964. [PubMed]
Gronemeyer H, Zelent A (2009) Fingering modulators of retinoic acid signaling identifies new prognostic marker for neuroblastoma. Cancer Cell 15: 249–251. [PubMed]
Harris TM, Childs G (2002) Global gene expression patterns during differentiation of F9 embryonal carcinoma cells into parietal endoderm. Funct Integr Genomics 2: 105–119. [PubMed]
Laudet V, Gronemeyer H (2002) The Nuclear Receptor Factsbook. San Diego: Academic Press.
Levine RA, Campisi J, Wang SY, Gudas LJ (1984) Butyrate inhibits the retinoic acid-induced differentiation of F9 teratocarcinoma stem cells. Dev Biol 105: 443–450. [PubMed]
Li R, Faria TN, Boehm M, Nabel EG, Gudas LJ (2004) Retinoic acid causes cell growth arrest and an increase in p27 in F9 wild type but not in F9 retinoic acid receptor beta2 knockout cells. Exp Cell Res 294: 290–300. [PubMed]
Liby KT, Yore MM, Sporn MB (2007) Triterpenoids and rexinoids as multifunctional agents for the prevention and treatment of cancer. Nat Rev Cancer 7: 357–369. [PubMed]
Lin CY, Chin CH, Wu HH, Chen SH, Ho CW, Ko MT (2008) Hubba: hub objects analyzer—a framework of interactome hubs identification for network biology. Nucleic Acids Res 36: W438–W443. [PMC free article] [PubMed]
Loudig O, Babichuk C, White J, Abu-Abed S, Mueller C, Petkovich M (2000) Cytochrome P450RAI(CYP26) promoter: a distinct composite retinoic acid response element underlies the complex regulation of retinoic acid metabolism. Mol Endocrinol 14: 1483–1497. [PubMed]
Loudig O, Maclean GA, Dore NL, Luu L, Petkovich M (2005) Transcriptional co-operativity between distant retinoic acid response elements in regulation of Cyp26A1 inducibility. Biochem J 392: 241–248. [PubMed]
Mahony S, Mazzoni EO, McCuine S, Young RA, Wichterle H, Gifford DK (2011) Ligand-dependent dynamics of retinoic acid receptor binding during early neurogenesis. Genome Biol 12: R2. [PMC free article] [PubMed]
Mark M, Ghyselinck NB, Chambon P (2006) Function of retinoid nuclear receptors: lessons from genetic and pharmacological dissections of the retinoic acid signaling pathway during mouse embryogenesis. Annu Rev Pharmacol Toxicol 46: 451–480. [PubMed]
Mikesch JH, Gronemeyer H, So CW (2010) Discovery of novel transcriptional and epigenetic targets in APL by global ChIP analyses: emerging opportunity and challenge. Cancer Cell 17: 112–114. [PubMed]
Min IM, Pietramaggiori G, Kim FS, Passegue E, Stevenson KE, Wagers AJ (2008) The transcription factor EGR1 controls both the proliferation and localization of hematopoietic stem cells. Cell Stem Cell 2: 380–391. [PubMed]
O'Malley BW, Kumar R (2009) Nuclear receptor coregulators in cancer biology. Cancer Res 69: 8217–8222. [PMC free article] [PubMed]
Orkin SH, Wang J, Kim J, Chu J, Rao S, Theunissen TW, Shen X, Levasseur DN (2008) The transcriptional network controlling pluripotency in ES cells. Cold Spring Harb Symp Quant Biol 73: 195–202. [PubMed]
Perissi V, Rosenfeld MG (2005) Controlling nuclear receptors: the circular logic of cofactor cycles. Nat Rev Mol Cell Biol 6: 542–554. [PubMed]
Reiss DJ, Facciotti MT, Baliga NS (2008) Model-based deconvolution of genome-wide DNA binding. Bioinformatics 24: 396–403. [PubMed]
Rosenfeld MG, Lunyak VV, Glass CK (2006) Sensors and signals: a coactivator/corepressor/epigenetic code for integrating signal-dependent programs of transcriptional response. Genes Dev 20: 1405–1428. [PubMed]
Saeed AI, Bhagabati NK, Braisted JC, Liang W, Sharov V, Howe EA, Li J, Thiagarajan M, White JA, Quackenbush J (2006) TM4 microarray software suite. Methods Enzymol 411: 134–193. [PubMed]
Saeed AI, Sharov V, White J, Li J, Liang W, Bhagabati N, Braisted J, Klapa M, Currier T, Thiagarajan M, Sturn A, Snuffin M, Rezantsev A, Popov D, Ryltsov A, Kostukovich E, Borisovsky I, Liu Z, Vinsavich A, Trush V et al. (2003) TM4: a free, open-source system for microarray data management and analysis. Biotechniques 34: 374–378. [PubMed]
Shankaranarayanan P, Rossin A, Khanwalkar H, Alvarez S, Alvarez R, Jacobson A, Nebbioso A, de Lera AR, Altucci L, Gronemeyer H (2009) Growth factor-antagonized rexinoid apoptosis involves permissive PPARgamma/RXR heterodimers to activate the intrinsic death pathway by NO. Cancer Cell 16: 220–231. [PubMed]
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T (2003) Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 13: 2498–2504. [PubMed]
Su D, Gudas LJ (2008) Gene expression profiling elucidates a specific role for RARgamma in the retinoic acid-induced differentiation of F9 teratocarcinoma stem cells. Biochem Pharmacol 75: 1129–1160. [PMC free article] [PubMed]
Taneja R, Roy B, Plassat JL, Zusi CF, Ostrowski J, Reczek PR, Chambon P (1996) Cell-type and promoter-context dependent retinoic acid receptor (RAR) redundancies for RAR beta 2 and Hoxa-1 activation in F9 and P19 cells can be artefactually generated by gene knockouts. Proc Natl Acad Sci USA 93: 6197–6202. [PubMed]
Uyttendaele H, Soriano JV, Montesano R, Kitajewski J (1998) Notch4 and Wnt-1 proteins function to regulate branching morphogenesis of mammary epithelial cells in an opposing fashion. Dev Biol 196: 204–217. [PubMed]
Vercauteren SM, Sutherland HJ (2004) Constitutively active Notch4 promotes early human hematopoietic progenitor cell maintenance while inhibiting differentiation and causes lymphoid abnormalities in vivo. Blood 104: 2315–2322. [PubMed]
Wang X, Wang TT, White JH, Studzinski GP (2006) Induction of kinase suppressor of RAS-1(KSR-1) gene by 1, alpha25-dihydroxyvitamin D3 in human leukemia HL60 cells through a vitamin D response element in the 5′-flanking region. Oncogene 25: 7078–7085. [PMC free article] [PubMed]
Ye Q, Shieh JH, Morrone G, Moore MA (2004) Expression of constitutively active Notch4 (Int-3) modulates myeloid proliferation and differentiation and promotes expansion of hematopoietic progenitors. Leukemia 18: 777–787. [PubMed]
Zhang L, Ren X, Alt E, Bai X, Huang S, Xu Z, Lynch PM, Moyer MP, Wen XF, Wu X (2010) Chemoprevention of colorectal cancer by targeting APC-deficient cells for apoptosis. Nature 464: 1058–1061. [PMC free article] [PubMed]
Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W, Liu XS (2008) Model-based analysis of ChIP-Seq (MACS). Genome Biol 9: R137. [PMC free article] [PubMed]

Articles from Molecular Systems Biology are provided here courtesy of
The European Molecular Biology Organization and Nature Publishing Group