Characterisation of Genome-Wide PLZF/RARA Target Genes

The PLZF/RARA fusion protein generated by the t(11;17)(q23;q21) translocation in acute promyelocytic leukaemia (APL) is believed to act as an oncogenic transcriptional regulator recruiting epigenetic factors to genes important for its transforming potential. However, molecular mechanisms associated with PLZF/RARA-dependent leukaemogenesis still remain unclear. We searched for specific PLZF/RARA target genes by ChIP-on-chip in the haematopoietic cell line U937 conditionally expressing PLZF/RARA. By comparing bound regions found in U937 cells expressing endogenous PLZF with PLZF/RARA-induced U937 cells, we isolated specific PLZF/RARA target gene promoters. We next analysed gene expression profiles of our identified target genes in PLZF/RARA APL patients and analysed DNA sequences and epigenetic modification at PLZF/RARA binding sites. We identify 413 specific PLZF/RARA target genes including a number encoding transcription factors involved in the regulation of haematopoiesis. Among these genes, 22 were significantly down regulated in primary PLZF/RARA APL cells. In addition, repressed PLZF/RARA target genes were associated with increased levels of H3K27me3 and decreased levels of H3K9K14ac. Finally, sequence analysis of PLZF/RARA bound sequences reveals the presence of both consensus and degenerated RAREs as well as enrichment for tissue-specific transcription factor motifs, highlighting the complexity of targeting fusion protein to chromatin. Our study suggests that PLZF/RARA directly targets genes important for haematopoietic development and supports the notion that PLZF/RARA acts mainly as an epigenetic regulator of its direct target genes.


Introduction
Retinoic acid receptors (RARs) belong to a family of nuclear receptors that can activate or repress transcription of target genes by recruiting co-activator or co-repressor complexes. It is the binding of its natural ligand, retinoic acid (RA), which transforms the activity of RAR from a repressor to an activator of transcription by inducing a conformational change in the ligand binding domain structure [1,2]. In acute promyelocytic leukaemia (APL), the RAR alpha (RARA) gene on chromosome 17, finds itself involved in chromosomal recombination with one of several different potential partner genes (designated ''X''), i.e. PML, PLZF, NuMA, NPM1, FIP1L1, PRKAR1A, STAT5b or BCOR, generating X/RARA fusion proteins [3,4]. Evidence to date suggests that X/ RARA induces a characteristic block in differentiation, associated with the APL phenotype through repression of downstream target genes [1][2][3][4]. The capacity of X/RARAs to form homo-oligomers and hetero-oligomers has been proven critical for inducing transformation and defines their oncogenic propensity [5][6][7][8]. However, it is still not well understood the extent to which the X/ RARA protein, as part of an oligomeric complex, has acquired altered DNA-binding capacity that may result in aberrant expression of genes normally (or not) regulated by wild-type RARA. Several studies have suggested that the presence of the partner X enforces a different conformation of the DNA recognition site for X/RARA versus RARA and that X/RARA would recognise degenerate RA response elements [8][9][10]. This phenomenon could be amplified by recruitment of co-factors that may also impose and/or change the DNA recognition specificity [11].
The nature of the fusion partner X has a decisive impact on the disease state and sensitivity to the therapeutic effects of all-trans RA (ATRA) and arsenic trioxide (ATO) [12]. Importantly, PLZF/ RARA-associated APL is resistant to ATO and exhibits impaired sensitivity to ATRA, giving rise to a significantly poorer clinical outcome compared to patients with classical PML-RARA+ disease, which typically responds to both of these molecularlytargeted therapies [13,14]. Studies conducted over a decade ago provided insights into the distinct natural history of these subtypes of APL, showing major differences in the capacity of the PML/ RARA and PLZF/RARA fusion proteins to bind corepressor complexes according to the level of retinoic acid. While corepressors are displaced from the PML/RARA oncoprotein in the presence of pharmacological doses of retinoid, binding of SMRT/NCoR-HDAC [15,16] and epigenetic factors such as Polycomb group (PcG) complexes [17] to PLZF/RARA persists under such conditions through interaction with the PLZF moiety of the fusion protein [18]. These specific recruitments provide an explanation for the ATRA-resistant phenotype and suggest there may be mechanistic differences in transcriptional repression mediated by PLZF/RARA as compared to other X/RARA fusions. Due to the PLZF moiety, PLZF/RARA could potentially control two different sets of genes: those genes normally regulated by RARA-RXR and those that are de novo target genes due to the sequence-specific DNA binding activity of PLZF. In addition, the situation is made more complicated by expression of the reciprocal RARA/PLZF that can potentially bind to PLZF binding sites where it may function as a transcriptional activator [19]. The development of genome-wide approaches such as ChIP-on-chip has favoured the identification of genomic targets that have widened our understanding of oncogenic transcription factor activities.
Here, we report a genome-wide ChIP-on-chip study of PLZF/ RARA gene targets using an inducible cell system. We identify 413 high-confidence specific target genes that provide clues as to how PLZF/RARA contributes to APL pathogenesis.

Results and Discussion
Genome-wide identification of specific PLZF/RARA target genes To identify gene targets directly regulated by PLZF/RARA, we took advantage of a zinc-inducible U937 cell system either expressing the PLZF/RARA fusion protein or harbouring an empty vector [20]. In the U937-B412 cell line, PLZF/RARA expression can be induced upon zinc induction ( Figure S1). We performed ChIP experiments using an anti-PLZF antibody that recognises both wild-type PLZF and the PLZF/RARA fusion protein. Immunoprecipitated DNA samples were hybridised to a custom promoter array containing ,18,000 human promoters and replicate experiments were merged (see Materials and Methods). Using a peak detection algorithm [21] we identified a total of 1545 significantly bound promoter regions in the U937-PLZF/RARA cell line. Among these targets, we found previously characterised PLZF targets, such as several genes from the HOXD cluster [17,22] (Figure 1a). Indeed, we detected high enrichment signals at the promoter regions of the HOXD4 and HOXD9 genes in both U937-PLZF/RARA and control cell lines. These data were corroborated by quantitative PCR analyses of independent PLZF ChIP samples (ChIP-qPCR) (Figure 1b). To identify PLZF/ RARA-specific target genes we compared the targets found in the zinc-stimulated U937-PLZF/RARA cell line to those found in the zinc-stimulated control U937-empty vector cell line and identified 412 promoter regions that were significantly bound only upon PLZF/RARA induction (Table S1). Among these genes, we identified RARB2 (the RARA fusion model target gene) as well as C/EBPE, ASB2, PRAM1 and IL8, that have all previously been reported to be regulated by X/RARA fusion proteins [23][24][25][26]. To further validate the specificity of our approach, we performed quantitative PCR analyses on ChIP samples (ChIP-qPCR) on a selection of ten PLZF/RARA specific target genes and one negative control (Fig. 1c). All tested PLZF/RARA-bound promoter regions were significantly amplified in the induced U937-PLZF/RARA cell line only (Fig. 1d), underscoring the specificity of the generated PLZF/RARA-binding profiles.
In a recent study, Rice and colleagues [26] identified more than 4900 genes as direct PLZF and/or PLZF/RARA targets. However, in this work the PLZF and PLZF/RARA targets were not discriminated, making it difficult to compare with the specific PLZF/RARA targets identified in the present study. Nevertheless, 168 specific PLZF/RARA target genes from our study were also identify as potential targets in the Rice study. The discrepancies between the two studies might be explained by the different experimental procedures and by the methods and the thresholds used to select bona-fide binding regions. Functional classification of specific PLZF/RARA target genes according to gene ontology (GO) analyses showed that PLZF/RARA binds to genes implicated in either basic cellular processes or developmental related processes (Table 1). Strikingly, a significant number of the gene products are involved in the regulation of transcription, highlighting the oncogenic potential that could arise from their deregulation, including transcription factors implicated in haematopoietic development, such as C/EBPA, C/EBPE, GATA1 and BCL6, as well as other general chromatin modifiers such as BRD2, MBD6 and SMARCD2 (Table 2).

Gene expression changes at PLZF/RARA target genes
To determine the transcriptional effect of PLZF/RARA on its identified target genes, gene expression profiling was performed in U937-MT (control) and U937-B412 (PLZF/RARA) cell lines. Expression profiling following PLZF/RARA induction was evaluated by microarray analysis and 4321 genes were identified with altered expression levels (.1.5 fold; P,0.05). A highly significant 20% of the genes targeted by PLZF/RARA have altered expression patterns (P = 4.99610 29 , Fisher exact) with 11% becoming down regulated and 9% being up regulated upon PLZF/RARA binding ( Figure 2a). This result is consistent with previous microarray analyses that have revealed that PLZF/ RARA induces a signature characterised by both down regulated and up regulated genes [24,26]. This suggests that while PLZF/ RARA can act as a transcriptional regulator, effects of PLZF/ RARA binding to chromatin cannot be measured only by transcriptional analysis. In addition, this shows that the direct transcriptional activity of PLZF/RARA is difficult to assess in an established cell line that has already shut down transcriptional pathways. To analyze the significance of PLZF/RARA targets identified by ChIP-on-chip with the APL phenotype we performed analyses of microarray gene-expression profiles from diagnostic samples from PLZF/RARA+ patients. These analyses identified 22 direct PLZF/RARA target genes that were significantly down regulated in PLZF/RARA APL blasts implicating their deregulation in tumour pathogenesis ( Figure 2b). The analysis confirmed that the ITGB2 gene encoding CD18, which is characteristically under expressed in primary APL cells (a feature used to help identify patients with PML-RARA+ disease) [27], was also down regulated in the presence of the PLZF/RARA fusion. Moreover, we demonstrated that several PLZF/RARA-target genes relevant for malignant haematopoiesis were significantly down regulated in PLZF/RARA patient samples, including IDH2, PAX5 and DAZAP1. IDH2 is frequently mutated in acute myeloid leukaemia and its mutations confer adverse prognosis in cytogenetically normal AML [28,29]. PAX5 and DAZAP1 are both involved in translocations associated with acute lymphoblastic leukaemia [30,31]. In addition, PLZF/RARA significantly down regulated   genes involved in self-renewal. This is the case for transforming growth factor beta-1 (TGFbeta-1), a pluripotent cytokine that controls key tumour suppressive functions. Interestingly, cytoplasmic PML has been described as a critical TGF-beta regulator and PML/RARA APL cells have defects in TGF-beta signalling similar to those observed in PML-null cells [32,33]. This is also the case for DUSP6, an inhibitor of the ERK signalling that has been described to be under expressed in different types of cancer, including pancreatic [34] and lung cancer [35], and which was reported as a potential PLZF-RARA target gene by Rice and colleagues [26].
The expression profiles obtained with the patient samples correlates with the expression profile obtained in B412 cell system as the majority of the genes down regulated in PLZF/RARA+ patients were also down regulated in B412 cells expressing PLZF/ RARA (Figure 2c). This integrated genome wide analysis thus highlights the importance of a subset of PLZF/RARA target genes and reveals new pathways that could be disrupted in the PLZF/ RARA APL.
Chromatin modification at PLZF/RARA target genes PLZF/RARA was shown to act as a dominant-negative inhibitor of RARA, mediating repression of retinoid target genes important for myeloid differentiation: PLZF/RARA binding to RARB2 has previously been shown to recruit PcG complexes, leading to an increase in trimethylation of Lysine 27 of histone H3 (H3K27me3) and RARB2 repression (ref. [17] and data not shown). To investigate whether PLZF/RARA can induce chromatin modification at its target genes, we analysed the trimethylation status of lysines 4 and 27 of histone H3 (H3K4me3 and H3K27me3) and the acetylation status of histone H3 (H3K9K14Ac) in response to PLZF/RARA-induction by ChIPon-chip. Indeed, PLZF/RARA induction provokes epigenetic variations represented as an increase in H3K27me3, a decrease in H3K4me3 or a decrease in H3K9K14Ac on 267 of its 413 target genes (Figure 3a). To evaluate the consequences of epigenetic changes on gene expression, we calculated the percentages of deregulated genes in function of individual or combinatorial changes for each histone modification (Figure 3b). We observed that the correlation between histone modifications and gene expression changes was higher when taking into account at least 2 of the analysed marks. Strikingly, genes displaying epigenetic changes for all three of the histone modifications displayed the higher percentage of repressed genes. These results strongly support the idea that PLZF/RARA protein could recruit different co-repressor complexes to its target gene promoters, resulting in a multiplicity of epigenetic changes.
To validate the transcriptional and epigenetic effects of PLZF/ RARA on its identified target genes, we analysed the mRNA expression level as well as changes in epigenetic histone modifications of a subset of PLZF/RARA targets. The five selected genes were repressed by direct PLZF/RARA binding, while their epigenetic changes were showing similar characteristics ( Figure 4). With the exception of ASB2, the genes down regulated by PLZF/RARA recruitment display an increase in H3K27me3 at their promoters with a concomitant decrease in H3K9K14Ac (Figure 4). H3K4me3 was either decreased or unmodified but was not found increased at PLZF/RARA repressed loci (Figure 4). Interestingly, ITGB2 and IL8 that are strongly repressed in B412 cell line and that display epigenetic changes for the 3 histone modifications (Figure 4), were down regulated in PLZF/RARA+ patient (Figure 2). These observations emphasize the importance of these two genes in APL pathology and highlight the relevance of epigenetic changes consistent with a PcG-repressed chromatin state that is mediated by PLZF/RARA binding.
In this study we show that PLZF/RARA expression induces epigenetic changes in a subset of genes of which expression is not necessarily modified. This observation is interesting in the context of the different nature of the fusion partner X that impact on the therapeutic effects of ATRA and ATO. Previous studies have demonstrated that X/RARA proteins recruit different partners that exert influences upon disease characteristics with respect to sensitivity to RA treatment [15][16][17][18]. In this study we reveal an epigenetic profile that is characteristic to PLZF/RARA bound loci. While a study has shown that high levels of H3K27me3 was associated with the PML/RARA binding region in the RARB2 promoter [36], genome-wide analysis of H3K27me3 at PML/ RARA binding sites reveals no significant increase of this repressive mark related to PML/RARA induction [3]. Thus, characterisation of epigenetic changes induced by X/RARA, even if not associated with detectable changes in gene expression, could provide an explanation for the difference in sensitivity to RA.

Sequence analysis of the PLZF/RARA-bound DNA
PLZF/RARA has conserved the DNA binding domain of RARA that binds to RA response elements (RAREs) defined by a direct repeat of the core motif 59-RGKTCA-39 separated by 1, 2 or 5 nucleotides (DR1, DR2, DR5) [9]. To determine whether the core motif RGKTCA (RARE half-site or RAREh) was the determining factor in recruiting PLZF/RARA, the sequences overlapping each peak were analyzed for the presence of this motif. We found that approximately half (233 out of 412) of the PLZF/RARA-bound DNA sequences contained at least one RAREh and that among these DNA sequences, almost 50% of them contained two or more repeats of the RAREh (Figure 5a). Spacing and orientation of core motifs dictate the specificity of DNA recognition by the nuclear receptors and the X/RARA fusion proteins [37]. We analysed for the presence of the possible combination of the core motif sequence; direct repeat (DR), inverse repeat (IR) or everted repeat (ER) in PLZF/RARA promoters where we could detect more than one RAREh. This showed that binding more frequently occurs at DR elements than at IR or ER motifs, with DR1, DR2 and DR5 accounting for 35% of promoter sites targeted by PLZF/RARA overall (Figure 5b). Thus, PLZF/RARA appears to bind to promoter regions harbouring either consensus or highly degenerated RARE, but the presence of RAREs is not absolutely required for PLZF/RARA to bind to its target promoters. However, we could not find a correlation between the presence of RARE in PLZF/RARA binding sites and PLZF/RARA target gene expression (Table S1). This result confirms previous observations suggesting that the binding of X/RARA to chromatin does not require direct consensus RARE [37]. This is due to the presence of the fusion  showed that even in the case of wild-type RARA, DR1, DR2 or DR5 elements do not account for the majority of bound loci, proposing that chromatin topology contributes to RAR occupancy [38]. A potential difference between PLZF as a fusion partner and most of the other APL fusion partners (e.g. PML) is that PLZF can bind DNA in its own right. Although the PLZF DNA binding domain is disrupted by the fusion, through heterodimerization PLZF/RARA could be targeted to PLZF binding sites. In order to find out whether significant binding sites can be identified in PLZF/RARA bound regions we searched for transcription factor (TF) binding sites using Genomatix (Table 3 and Figure S2). We found that PLZF/RARA target regions were significantly enriched for several TF binding sites, including the KLF and ETS family of TFs. However, we did not find that PLZF/RARA target regions were specifically enriched for PLZF binding sites. Interestingly, PML/RARA was recently shown to be targeted to genes regulated by PU.1, a member of the ETS family [4]. Thus, it is possible that like PML/RARA, PLZF/RARA collaborates with other tissuespecific TFs for DNA binding.

Conclusion
By comparing PLZF bound regions between U937 and PLZF/ RARA-induced U937 cells we identified 413 specific PLZF/RARA target genes. This selection of bona-fide PLZF/RARA targets was significantly enriched in genes important for hematopoietic development and included genes effectively downregulated in primary PLZF/RARA+ APL cells. Moreover, PLZF/RARA binding was preferentially associated with gene repression at sites enriched for the repressive H3K27me3 histone modification mark. Analysis of PLZF/RARA bound sequences suggests PLZF/RARA binds canonical and degenerate RARE motifs, and recruitment of PLZF/RARA may also occur via other hematopoietic specific TFs. This study helps to clarify the role of the PLZF/RARA fusion protein as an oncogenic transcriptional regulator during leukaemogenesis and provides further insights into the pathogenesis of the disease.

Cell lines
Human myelomonoblastic cell lines U937-MT and U937-B412 [20] were maintained at exponential growth in RPMI supplemented with 10% foetal calf serum. U937-MT is the empty vector control, U937-B412 contains PLZF/RARA cDNA under the control of the zinc inducible human-metallothionein promoter [20]. For PLZF/RARA induction cells were prestimulated with 0.1 mM ZnSO 4 for at least 24 h.

RNA isolation and RT-PCR
Total RNA was extracted using the RNeasy Mini kit according to the manufacturer's protocol (Qiagen), For cDNA synthesis, reverse transcription was performed with 2 mg of the total RNA, oligo dT, dNTPs, DTT, buffer and Superscript Retrotranscriptase II (Invitrogen). cDNA was analyzed by qPCR using Brilliant

ChIP-on-chip
Total DNA and ChIP DNA were amplified using the whole genome amplification kit (WGA-2, Sigma-Aldrich, Taufkirchen, Germany) and fluorescently labelled using BioPrime Array CGH Genomic labelling kit (Invitrogen). Labeled samples were hybridized following the manufacturer's instructions to a custom promoter array (Agilent, Santa Clara, USA) containing 236.992 probes spanning the promoter regions (22,5 kb to +1 kb) of 18.000 refseq genes (UCSC assembly HG17). Experiments were performed in duplicate and showed high correlation (R$0.769). All data is MIAME compliant. Median-normalized log2 enrichment ratios (ChIP/Input) of merged replicates were calculated using CoCAS software [21]. For PLZF ChIP-on-chip experiments, significantly enriched promoters were isolated using the peak detection algorithm in CoCAS and the following settings: 26 and 16 standard deviations for main and extended thresholds, respectively. Enrichment scores were calculated by measuring the effective peak area [21]. For histone modification experiments, enrichment scores were calculated by computing the average signals at each promoter region. Next, fold changes between MT and B412 chromatin at PLZF/RARA target promoters were calculated for each epigenetic mark. Gene promoters displaying a differential enrichment of .0,135 (P val,0.05) for a given histone modification, were considered to be differentially marked.

Functional annotation of genes
Target genes were classified using DAVID [39] and according to PANTHER molecular functions and biological process categories. Only the highest significant categories were selected (P,0,001).

Sequence analysis
We looked for the consensus RARE motif using Perl regular expressions (http://www.perl.org) on ChIP-on-chip peak sequences. The expression code was built according to the known repeat pattern (A/G)G(G/T)TCA. We allowed a gap ranging from 0 to 100 bases and performed the calculation for finding DR, IR and ER on both strands. Analyses of transcription factor binding site enrichments were performed using MatInspector tool from GENOMATIX software (http://www.genomatix.de).

Microarray analysis
MT and B412 cells were cultured in duplicate in the presence of zinc for 48 hours. RNA was extracted using QIAGEN RNeasy and quality was tested using an Agilent Bioanalyzer. 1 mg of RNA was hybridized to U133 Plus 2.0 arrays (Affymetrix). Affymetrix U133A.CEL data files from B412 cells and MT cells were imported into R and probe intensities were normalized using GC-RMA. A total of 4321 genes were found to be differentially expressed (.1,5 fold; P,0.05). For patient samples, microarray raw data was taken from Gene Expression Omnibus (http://www.  ncbi.nlm.nih.gov/projects/geo/), under the data set GSE8510 for primary PLZF/RARA+ APL patient cases and under the data set GSE15061 for normal bone marrow samples (NBM). Details of the patient samples analysed have been reported previously [19]. Expression data were Robust Multichip Average (RMA) [40] normalised using the 'affy' R/Bioconductor library. The differential gene expression study was performed by Significance Analysis of Microarrays using MeV software [41]. Figure S1 Expression of PLZF/RARA, PLZF and RARA proteins in MT and B412 cell lines. Nuclear extracts were prepared from U937-MT (1) and U937-B412 (2) treated with zinc and were immunobloted for PLZF/RARA and PLZF (a-PLZF) or RARA (a-RARA).

Supporting Information
(TIF) Figure S2 Transcription factor binding sites overrepresented in PLZF/RARA target promoters. MatInspector tool from GENOMATIX was used to identified transcription binding sites in PLZF/RARA bound promoters.