Identification of p53-target genes in Danio rerio

To orchestrate the genomic response to cellular stress signals, p53 recognizes and binds to DNA containing specific and well-characterized p53-responsive elements (REs). Differences in RE sequences can strongly affect the p53 transactivation capacity and occur even between closely related species. Therefore, the identification and characterization of a species-specific p53 Binding sistes (BS) consensus sequence and of the associated target genes may help to provide new insights into the evolution of the p53 regulatory networks across different species. Although p53 functions were studied in a wide range of species, little is known about the p53-mediated transcriptional signature in Danio rerio. Here, we designed and biochemically validated a computational approach to identify novel p53 target genes in Danio rerio genome. Screening all the Danio rerio genome by pattern-matching-based analysis, we found p53 RE-like patterns proximal to 979 annotated Danio rerio genes. Prioritization analysis identified a subset of 134 candidate pattern-related genes, 31 of which have been investigated in further biochemical assays. Our study identified runx1, axin1, traf4a, hspa8, col4a5, necab2, and dnajc9 genes as novel direct p53 targets and 12 additional p53-controlled genes in Danio rerio genome. The proposed combinatorial approach resulted to be highly sensitive and robust for identifying new p53 target genes also in additional animal species.

Scientific RepoRts | 6:32474 | DOI: 10.1038/srep32474 network 15,17,18 . Therefore, the identification and characterization of a species-specific p53 BS consensus sequence and of the associated target genes may help to provide new insights into the evolution of the p53 regulatory networks across different species. Here, we designed and biochemically validated a computational tool to identify p53-target genes in Danio rerio (zebrafish). Our approach expanded the number of direct p53 target genes in this organism.

Results
In silico similarity search of putative p53-binding sites. A total of 3,033 hits for the RRRCWWGY YYN{0, 15}RRRCWWGYYY motif were found in the Danio rerio genome (Table S1), according to the PatMatch scanning analysis, with a substantially uniform pattern distribution among the 25 chromosomes (2.25 patterns per Mb, on the average). This analysis detected the 47% (1,420 out of 3,033) of the retrieved p53 DNA binding sites within repeat or transposable elements of Danio rerio genome. We ranged the spacer size from 0 to 15 bp and found a bi-modal frequency distribution, with peaks at 0 and 11 bp. A total of 1,277 gene-related patterns were found located in the ± 5 Kb genomic neighbourhood of Danio rerio genes. In some cases, multiple genes resulted to be close to the same pattern or multiple patterns were found in the proximity of the same gene (Table 1 and  S2). Particularly, 1,099 patterns were uniquely mapped to distinct genes, while 178 were proximal to more than one gene. Totally, this analysis identified 984 pattern-related genes including 979 protein-coding and 5 miRNAs. 140 out of 984 Danio rerio genes were associated to more than one pattern. Such genes as arl13a, cetn3, rnf150a, runx1, zgc:136572, zgc:158766 have 10 or more motifs within their neighbourhoods. 472 out of 984 genes were positioned on the same genomic strand of their corresponding patterns (Table 1). All the genomic features of the found patterns are described in Tables S1 and S2.
Thus using the prioritization criteria (see Material and Methods section) (Table S3), along with the information about tissue specificity and the position of the patterns relative to the TSS, we selected a subset of 134 candidate motifs (Table S4).
As a preliminary investigation, we searched human orthologs for 134 prioritized zebrafish genes (Table S4) and checked their expression level into p53-related public datasets (Gene Expression Atlas 3.0, access: June 2016, https://www.ebi.ac.uk/gxa/home) (Table S5). A considerable proportion of Dr-retrieved gene/human-gene sets (53%, 73%, and 15,7% for p53, p53 and cancer, doxorubicin condition queries, respectively) was reported to be under-or over-expressed according to different experimental conditions suggesting their possible role as putative p53-controlled genes.
p53-binding sites in vitro validation. We randomly selected 10 out of 134 predicted Danio rerio BSs, located at varying distances from the TSS of the related gene, for further in vitro validation ( Table 2 and Table S6).
First, to assess whether in silico candidate patterns may drive the Drp53-mediated transactivation or repression, we performed a luciferase-based assay. Short genomic regions harbouring selected Drp53 BSs associated to 10 genes ( Table 2 and Table S6) were cloned in luciferase reporter plasmids, and transfected into the human p53-deficient cell line H1299 19 along with vectors expressing wild-type Drp53, or an empty vector as a negative control. All amplified genomic fragments contain only one predicted p53 BS, but the region associated to runx1 gene that includes four p53 BSs. A vector containing trim8a p53 REs was used as a positive control 19 . We found that all putative BSs mediate p53-responsiveness (Fig. 1a), although with different strength. Particularly, reporter gene activation mediated by Drp53 was stronger for the p53 BSs potentially associated to runx1 (60 folds), axin1 (25 folds), traf4a (50 folds), and hspa8 (9 folds) genes, when compared to the control. To test whether the activation was dependent on functional p53, we used an N-terminal truncated p53 isoform, Δ 1-113 Drp53 (pCS2-Drp53Mut), which interferes with the p53 functions 20 . Consistent with the effect of this truncated form of p53, co-transfection of both full-length Drp53 and Δ 1-113 Drp53 constructs in H1299 cells result in a luciferase activity decrease for all analysed BSs, when compared to transfection with full-length Drp53 protein only (Fig. 1a).
Next, we generated mutant vectors carrying the deletion of BSs for those genomic fragments that exhibited strong reporter gene activation. As showed in Fig. 1b for all constructs, except for p53-delRE-runx1, which was still responsive to the action of p53, likely because this construct still maintains three out of four predicted BSs. To define whether the conserved p53 RE CWWG core motif is necessary for reporter activation in Danio rerio as well as in human, we performed luciferase assays by using reporter vectors carrying nucleotide changes in the CWWG consensus (Fig. 1c). The results showed that replacing the core CWWG of runx1, axin1, traf4a, and hspa8 p53 REs with ACCT sequence the responsiveness to p53 was dramatically reduced (Fig. 1d). Taken together, these results revealed that the analysed Danio rerio p53 BSs are functionally active, being each of them capable of conferring p53-dependent transcriptional activation on a luciferase reporter gene. Accordingly to the human p53 RE consensus, the core CWWG motif has a fundamental role in determining p53 activation function. In order to confirm the functionality of these p53 BSs we investigated the native binding of endogenous p53 to the genomic locations where we have identified putative p53 BS. Chromatin immunoprecipitation in wild type 24 hour Danio rerio embryos (Fig. 1e) confirms that p53 binds to the runx1, traf4a, and hspa8 p53 REs (ChIP). We also observed a slight enrichment in the binding to axin1 p53 RE with respect to the control, probably reflecting the low level of axin1 expression at this developmental stage. Additionally, we detect clear enrichment of p53 in the col4a5, dnajc9, and necab2 p53 REs. This result indicates that the described p53 REs can act as enhancers through direct binding of p53 in the predicted binding sites, confirming that these seven genes are direct targets of p53. As a positive control for ChIP we analysed the binding of p53 to the p53-responsive element of a previously reported Drp53 target gene, puma 21 . Conversely, we did not find any direct p53 binding in the stat3, vcp, and ttyh2l p53 RE.
To verify if the genes embedding the p53 BSs are per se responsive to the action of p53, we profiled in vivo the expression of ten p53 RE associated genes on 54 hpf zebrafish embryos incubated for 16 hours in presence of 50 μM R-roscovitine, a cyclin-dependent kinase inhibitor that can efficiently stabilize and activate the nuclear p53 in human and zebrafish cells 22,23 . We found that the transcript levels of the genes surrounding the ten p53 BSs were increased in embryos upon exposure to R-roscovitine, compared to untreated embryos (Fig. 1f). Consistently, we observed a reduced mRNA level for all 10 genes in Danio rerio embryos injected with p53 antisense morpholino 24 , when compared to uninjected embryos (Fig. 1g).
Moreover we evaluated the responsiveness to human p53 of tested Dr BSs. We observed that the human p53 protein is able to efficiently activate the reporter gene expression (Fig. 1h), while the transactivation defective mutant p53-R175H (pcDNA3-MutHsp53) was unable to bind and activate the reporter constructs. Overall, these results confirm that p53 REs were responsive to both human and Danio rerio p53 proteins.
PWM-driven similarity search of Drp53-target genes and in vitro validation. Subsequently, by using an additional similarity-based motif-search computational strategy we identified a second list of candidate p53 pattern-related genes for in vitro validation.
We aligned the nucleotide sequences of the ten responsive elements initially validated by luciferase reporter assays (Fig. 2a,b, Table 2) calculating a new PWM (Table S7), which was finally inputted to two popular motif analysis resources, MEME and RSAT, and used as seed for an additional pattern similarity search. This strategy identified a group of 52 putative p53 pattern-related genes (Table S8), 21 of which were randomly selected for further in vitro validation (Table 3).
To verify whether such 21 genes are per se responsive to p53, we performed a gene expression profiling by qPCR in 54 hpf Danio rerio embryos incubated for 16 hours in presence of 50 μM R-roscovitine (Fig. 3a) or in embryos injected with p53 antisense morpholino (Fig. 3b), finding a significant increase of transcript levels in 12 out of 21 analysed genes or an mRNA decrease of all genes, respectively, compared to untreated embryos (Fig. 3a,b).    We also investigated the effects of our computational strategy for pattern prioritization by constructing two additional matrices and comparing results obtained. We randomly chose other two sets of sequences from the list of 134 candidate genes (Table S4) and created two further PWMs ( Fig. 4 and Table S7). The selected patterns (Table S4, columns "PWM-II" and "PWM-III") were analysed with the same combined strategy of motif similarity search, MEME and RSAT, applied to the first PWM, which was obtained from the first 10 validated motifs. Comparing the results obtained with the three PWMs, we found that 12 and 17 out of the 52 selected patterns (which are those with at least one significant similarity score calculated by the first PWM) were obtained when applying the two methodologies on the PWM-II, respectively. The analysis of the third PWM yielded, instead, 22 and 21 patterns from MEME and RSAT implementation, respectively (Table S8). The whole computational and biochemical analysis workflow is summarized in Fig. 5.
The functional role of all the 31 investigated genes is summarized in Table S9, in which, for each gene, Gene Ontology Accession numbers and Terms have been reported.

Discussion
The capacity of transcriptional regulatory machinery to quickly bring about changes in the gene expression pattern allows organisms to adapt to environmental and genetic changes. Transcription factor binding sites are generally known to evolve more rapidly than transcription factors themselves, playing a crucial part in the phenotypic divergence between and within closely related species 25,26 .
The binding of p53 with the DNA represents the crucial step by which the so-called "guardian of the genome" coordinates and regulates a wide set of downstream target genes. Because of the degenerate nature of the consensus p53 RE, changes in the level of p53 responsiveness have been observed and the rapid evolution of these regulatory elements determines variations in the p53 network between and within species 13 .
Therefore, the characterization of species-specific p53 BSs and target genes may help to provide insights into the evolution of p53 regulatory networks and to understand species-specific differences in several biological processes.
Here, we designed and biochemically validated a computational analysis that allowed us to identify novel p53 target genes in Danio rerio genome, an organism model widely used in human diseases and cancer, for whom little is known about the Drp53-mediated transcriptional signature. Consensus sequence logo has been created by using the WebLogo tool: http://weblogo.berkeley.edu/logo. cgi 63 . A total of 52 patterns were retrieved by MEME and RSAT motif similarity search (details shown in Supplementary Table S4).
Screening all the Danio rerio genomic contexts by pattern-matching-based analysis, we found a total of 3,033 matching p53 RE-like patterns with 1,099 of them proximal to the annotated Danio rerio genes. In this analysis 47% (1,420 out of 3,033) of the retrieved patterns were detected within repeat or transposable elements in Danio rerio genome, as previously reported in ref. 19. In particular, the query pattern partially matched the elements belonging to the hAT-AC (30% of the total pool) and DNA (14%) repeat families or to the TE-X-4_DR transposable element (14%). These elements can promote the spread of p53 BSs, thereby contributing to the diversity of the p53 response between and within species, including p53-influenced processes, such as tumour formation, longevity, angiogenesis and fertility [27][28][29] .
Among the predicted gene-related patterns, we selected 10 binding sites and 10 matching pattern-related genes for in vitro validation using luciferase and qPCR assays, respectively, confirming their responsiveness to both Drp53 and Hsp53. ChIP assays, performed in wild type 24 hours Danio rerio embryos, confirm that p53 directly binds to the runx1, traf4a, axin1, hspa8, col4a5, dnajc9, and necab p53 RE indicating that these genes can be considered as bona fide direct targets of p53. Drp53 was not detected to stat3, vcp, and tty2l p53 RE that exhibited the lowest level of gene reporter activity for both Drp53 and Hsp53.
The failure to detect direct binding of p53 to some of the p53 RE such as stat3, vcp, and tty2l might be due to technical issues such as low affinity sites, absence of cofactors or developmental stage as is shown by the low level of reporter activity for both Drp53 and Hsp53.
Notably, we observed a difference in reporter gene activation between human and zebrafish p53. This could be due to several reasons: i) Drp53 REs were searched by using the human motif as starting query. It is thus conceivable that the selected Drp53REs could be more specific for human p53 binding with respect to Drp53; ii) the Drp53 and Hsp53 are temperature sensitive proteins exhibiting their optimal activities at 28 °C and 37 °C, respectively, with Drp53 displaying no transactivation or very weak activation at 37 °C 19 ; iii) we carried out the luciferase assays in human cells cultured at 28 °C, which is not the standard condition for cell growth.
The validated pattern-gene pairs were subsequently used as seed for an additional similarity-based motif-search strategy that identified a second list of putative p53 targeted genes, among which we randomly chose 21 for further validation. We profiled their expression by qPCR both in Danio rerio embryos either silenced or enhanced for p53 expression, funding 12 additional potential p53 regulated genes.    Therefore, we cannot exclude the possibility that the other 9 genes, which are significantly deregulated in Danio rerio embryos that were injected with p53 antisense morpholino, but not in embryos treated with R-roscovitine, could be regulated by Drp53.
Many factors could contribute to the differential transactivation to p53 observed in our assayed genes. Post-transcriptional and post-translational modifications of p53, recruitment of co-factors (co-repressors/ chromatin-modifying factors), competition with transcription activator for DNA-binding site 30 , cellular and genomic context 31 , chromatin state, complexity of transcriptional regulation in living cells, features of the RE sequences, which take into account, i.e., of the variability of the spacer lengths, and the response variability upon each different p53 stimuli, as well as the R-roscovitine that interferes with the preferential activation of some p53-target genes, are possible causes. Regarding the latter point, we observed that the expression of dnai2a and mfge8b was down regulated in both R-roscovitine-treated and p53-antisense-injected embryos (Fig. 3a,b).  Probably, a different p53-dependent transactivation of some genes could be observed in response to several DNA-damaging reagents activating specific p53 transcriptional programs.
Collectively, our analysis identified 7 novel p53 target genes: runx1, axin1, traf4a, hspa8, col4a5, necab2, and dnajc9. Interestingly, the p53 target genes identified in this study appear to involve broad functional classes and multiple gene families that are related, directly or indirectly, to the p53 pathway in human.
Human RUNX1 is a member of the RUNX transcription factor family that play essential roles in the vertebrate embryo by directly regulating transcription in a variety of developmental pathways including neurogenesis, bone development, and segmentation 32 . Runx1 has also an important role during hematopoiesis in fish as in mouse and human 33 . Interestingly, Masse et al. reported that a functional interplay between p63 and p53 regulates RUNX1 expression in the control of the transition from proliferation to early differentiation in human keratinocytes 34 .
TRAF4 is a member of the TRAF family of adaptor proteins that mediate cellular signalling by binding to various members of the tumour necrosis family receptor superfamily and interleukin-1/Toll-like receptor super-family. The analysis of traf4-deficient mice revealed embryonic lethality underlying its potential importance during embryogenesis 35 . TRAF4 is specifically regulated by p53 in response to overexpression of p53 by use of an adenovirus and stabilization of p53 in response to exposure to DNA-damaging agents. The murine Traf4 genomic locus contains a functional p53 DNA-binding element located approximately 1 kilobase upstream of the translation start site. Overexpression of TRAF4 induces apoptosis and suppresses colony formation in multiple tumour cell lines suggesting a role for TRAF4 in determining cell fate in response to stabilization of p53 36 . Danio rerio ortholog traf4a shows a highly restricted expression pattern mainly in sensorial and neural cells, and in somites of Danio rerio embryos 37 .
Axin1 is a negative regulator of axis formation in the development of mouse embryos. Its deficiency leads to axis duplication 38 and its overexpression blocked embryo axis formation in Xenopus and caused apoptosis in transgenic mice 39,40 . Accumulating data showed that Axin1 controls multiple important pathways, including the canonical Wnt pathway, JNK signaling, TGF-β signalling, and p53 activation cascade [41][42][43] .
A large variety of cellular functions have been attributed to heat shock protein HSPA8, most of them through its cooperation with cochaperones like its role as a clathrin-uncoating ATPase during clathrin-mediated endocytosis. HSPA8 is also a major actor in chaperone-mediated autophagy and its involvement in protein import into organelles or cellular compartments has been the focus of numerous studies. HSPA8 has been identified as a p53-repressed target 44 . p53 binds a direct specific sequence on the HSPA8 promoter followed by p53-dependent recruitment of Sin3B/HDAC1 co-repressor complex as well as hyper-methylation. On the other hand, HSPA8 repression is critical for the functional activation of p53 because Hspa8 protein is known to antagonize the p53 nuclear localization by masking the nuclear localization signal sequence of p53 45 .
In conclusion, the approaches described here illustrate the great advantage of combining global p53 DNA-binding site in silico analysis and molecular assays, providing an efficient strategy to identify new p53 Danio rerio REs and target genes.
The identification and characterization of a Danio rerio-specific p53 BS consensus sequence and of the associated target genes may help to provide new insights into the evolution of the p53 regulatory networks in determining phenotypic diversity among species.

Materials and Methods
Pattern searching in Danio rerio genome. The first step of our computational strategy, which we have set up to identify the potential motifs of the p53 responsive elements in Danio rerio, encompasses the search of the human motif RRRCWWGYYYN{0, 15}RRRCWWGYYY in the Danio rerio genome. The build Zv9/danRer 7 version, July 2010, was considered as the reference genome (available at ftp://ftp.ensembl.org/pub/release-79/ fasta/danio_rerio/dna/) 46 , from which we retrieved the genomic coordinates of the genes and of the repeat/transposable elements.
The pattern-matching tool initially used was PatMatch 1.2 (available at ftp://ftp.arabidopsis.org/home/tair/ Software/) 47 . It implements an exact matching algorithm. It was run with this command line: perl patmatch.pl -n "RRRCWWGYYYN{0, 15}RRRCWWGYYY" Zv9.fasta.prepared > patmatch_Zv9_patterns.out, where the argument "-n" defines the pattern, while the "prepared" fasta file consists of the input genomic sequences (strand ' + '), manipulated as required by the tool. The expression "{0, 15}" indicates the spacer size. PatMatch produced a tabular output, containing all the nucleotide sequences that matched the pattern, together with the associated chromosomal coordinates. Results obtained by PatMatch were confirmed with p53scan 1.05 (available for download at www.rimls.nl/bioinfo/p53scan/) 48 . This tool implements a position weight matrix-based pattern-finding algorithm and uses a probabilistic approach that associates every matching site to a multinomial distribution (site-specific relative nucleotide frequencies, as defined by the p53scan PWM), rather than a single character state (one of the four nucleotides). It was run on the human p53 binding-site position weight matrix (PWM) accompanying the software package, by setting a spacer size ranging from 0 to 15. The command line was: python p53scan.py -i Zv9.fasta -s 0-15 > p53scan_Zv9_patterns.out, where "Zv9.fasta" specifies the Danio rerio reference genome, while the "-s" parameter sets the allowable sizes of the p53 motif spacer.
Then, we explored the genomic context of each found pattern by choosing 5 Kb up and downstream following [48][49][50] and focused only on contexts that included at least one annotated gene, according to the Zv9 reference gene coordinates. For each gene, we collected NCBI and Ensembl accession numbers, gene official symbol, CDS and exons coordinates, associated Gene Ontology terms via the Biomart query tool 51 . We thus annotated them in order to improve the sensitivity of the overall analytical strategy. First, we checked if putative target genes were located on the same genomic strand of the retrieved patterns. Considering also the distance from the TSSs, we labelled the patterns with the gene regions where they fell in (i.e., promoter, UTR, exon and intron). Hence, the Homo sapiens orthologs of all the genes falling within the ± 5 Kb up/down-stream genomic context of each pattern were searched in Ensembl 78, through Biomart, taking care of specifying that these were Danio rerio Ensembl Gene IDs. Homology relationships were obtained by extracting the gene trees of interest from the Ensembl Compara pre-computed phylogenies. The output consisted in a list of pairs of Ensembl Human Gene ID and Homology Type, comprising of only homolog genes with a one-to-one relationship. Among the different options of orthology predictions ("one2one", "one2many", "many2many"), we opted for the former, hypothesizing that such relationship should preserve the functional similarity of orthologs between Danio rerio and Homo sapiens.
Then, the resulting human ortholog genes were prioritized according to their functional proximity to 127 renowned p53-regulated human genes (Table S3), obtained from Riley et al. 9 . This set of 127 genes was basically used as training dataset for two gene prioritization bioinformatics tools: Endeavour 52 and ToppGene 53 , which were run with default settings. We then considered the first 50 genes, which exhibited the most significant overall scores by Endeavour or ToppGene (Table S3).
In addition, information on gene expression was retrieved from the Zfin web resource 54 by searching for tissue expression data for wild type organisms, from zygote to adult samples, through the "Expression" menu of the main page of the tool.
In summary, we first collected genomic and functional data related to each discovered "pattern-gene" pairs and then, we manually curated the selection of a limited number of them for in-vitro/in-vivo validation.

Cloning of Danio rerio p53 binding sites, site-directed mutagenesis and dual-luciferase assay.
Fragments of approximately 300-500 bp containing the predicted p53 BS were amplified by polymerase chain reaction (PCR) using a Danio rerio genomic DNA as template and then cloned into pGL3-basic vector (Promega). All constructs were verified by sequencing.
The human lung adenocarcinoma p53-null cell line H1299 was grown in DMEM medium supplemented with 10% fetal bovine serum (FBS; Life Technologies, USA) and 1% antibiotics mixture in a 37 °C incubator with 5% CO 2 . The reporter construct, p53 expression vector (pCS2-Drp53 or pcDNA3-Hsp53) or mutant p53 expressing vector (pCS2-Drp53Mut or pcDNA3-Hsp53Mut, an N-terminal truncated p53 isoform, Δ 1-113 Drp53, that lacks the transcription activation domain) or empty control vector, and pSV-Renilla (Promega) were co-transfected into p53-null H1299 cells plated in 96-well plates using Lipofectamine ® LTX with Plus ™ Reagent (Life Technologies) according to the manufacturer's instructions. Since it was previously reported that Drp53 is a temperature sensitive protein 18 , cells were grown at 28 °C and 37 °C when transfected with pCS2-Drp53 or pcD-NA3-Hsp53 respectively. 48 h after transfection, cells were washed with PBS, lysed and assayed for both Firefly and Renilla luciferase activity using the Dual-GLO ® Luciferase Assay System (Promega) by a Glomax 96 microplate luminometer. Firefly luciferase activity was normalized to the Renilla luciferase activity for each transfected well, as reported in ref. 19. Constructs changing the p53 BS were generated by QuickChange II site-directed mutagenesis kit (Stratagene). The primer pairs used are listed in the Table S10.

Chromatin immunoprecipitation.
ChIP experiments were carried out following the EZ-Magna ChiP A/G Chromatin Immunoprecipitation Kit (Merck-Millipore) protocol essentially as described in ref. 21. Briefly, 30 dechorionated embryos (24 to 34 hour after fertilisation, AB strain) were fixed for 10 minutes in 1% formaldehyde, washed and sonicated in 300 μ l of lysis buffer. Chromatin was immunoprecipitated with the anti-p53 antibodies kindly provided by Prof. G. Del Sal (University of Trieste, Trieste, Italy). The positive and negative controls were Rabbit anti-H3 (Abcam) and Rabbit anti-IgG (Diagenode) respectively. PCRs were performed with the Kapa 2G Fast HS (Kapa Biosystems) with the oligonucleotides indicated in Table S10. The data presented are representative of three independent experiments. Zebrafish care. Wild-type zebrafish (Danio rerio, AB type) embryos were obtained by natural spawning and were raised under standard conditions at 28 °C, with a 14-hour light/10-hour dark cycle in fish water (0.1 g/L Instant Ocean Sea Salts, 0.1 g/L sodium bicarbonate, 0.19 g/L calcium sulphate, 0.2 mg/L methylen blue, H 2 O) until the desired developmental stage was reached.
Embryos were staged according to Kimmel et al. 55 and developmental stages of zebrafish embryos were expressed as hpf or dpf (hours or days post-fertilization, respectively) at 28 °C. This study has been approved by OPBA (Organismo Preposto al Benessere degli Animali) of the Ethics Committee of the University of Brescia. The use of the teleost fish Danio rerio -zebrafish -for the study of human disease has been approved by Ministry of Health. Animal experiments were carried out according to the relevant guidelines and regulations of the EU Directive 2010/63/EU. Morpholino injection. p53 Morpholino Oligo (MO) (GeneTools LLC, Philomath, OR, USA) targeting p53 zebrafish translation start site was used (GCGCCATTGCTTTGCAAGAATTG). p53 MO was diluted in 1x Danieau buffer (pH 7.6) and 2.5% tetramethylrhodamine dextran (Molecular Probes) as a tracer. 0.5 pmol of p53 expression-inhibiting morpholino were injected into the yolk of 1-2 cell stage embryos as previously reported 24,56 . Embryos were then raised until 30 hpf and processed for RNA extraction as described below. Injections were repeated three times; 100 embryos were injected and processed in every experiment. Non-injected embryos were used as negative control.
Roscovitine treatment, RNA isolation and quantitative PCR. 54  kit (Qiagen) and reverse transcribed using QuantiTect ® Reverse Transcription Kit (Qiagen), according to the manufacturer's instructions. The expression levels of putative genes associated to p53 BS were examined by quantitative PCR (qPCR). The primer pairs were designed using the primer express software package 57 with default parameters (Table S10). top1 was used as housekeeping gene (Table S10). The reactions were run in triplicate in 10 μ l of final volume with 10 ng of sample cDNA, 0.3 mM of each primer, and 1x Power SYBR Green PCR Master Mix (Termo Fisher Scientific-Applied Biosystems). Reactions were set up in a 384-well plate format with a Biomec 2000 (Beckmann Coulter) and run in an ABI Prism 7900HT (Termo Fisher Scientific-Applied Biosystems) with default amplification conditions. Raw Ct values were obtained using SDS 2.4.1 (Termo Fisher Scientific-Applied Biosystems). Calculations were carried out according to the comparative Ct method as reported in ref. 19. Identification of additional potential targets by pattern similarity search. The validated pattern-gene pairs were subsequently used as seeds for an additional similarity-based motif-search computational strategy. In detail, by using the CLC Sequence Viewer 7.5 program (http://www.clcbio.com/products/ clc_sequence_viewer/), we aligned the nucleotide sequences of the responsive elements initially validated by luciferase reporter analysis. Therefore, since the consensus sequence of the p53 RE is degenerate, we used the aligned sequences to construct two new positional weight matrices (PWMs), one for each decamer, which were finally inputted to two popular motif analysis resources: MEME 58 and the Regulatory Sequence Analysis Tools (RSAT) 59 . We run MAST 4.9.1 of the MEME suite with default parameters on both decamers, separately. MAST returned a series of top scoring input sequences, ordered by increasing E-values (i.e., the expected number of sequences in a random database of the same size that would match the motifs as well as the sequence, with 10 as cut-off value). In addition, we used the "matrix-scan QUICK" method of the RSAT with default parameters. The Danio rerio-specific nucleotide frequencies were used as background. As a result, we selected the best significant (i.e., p-value < 0.05) "input matrix vs. input sequence" alignments. The sites involved in the pairwise alignments were the alignment positions "1-10" or "24-34". Finally, we crossed the top-ranked binding sites obtained with the two tools and made a second list of candidate patterns for in-vitro validation. This procedure was applied to two further matrices, built on 20 randomly selected patterns picked from the list of 52 "prioritized" patterns/ genes. With this approach, we evaluated the degree of consistency of the results obtained with MEME and RSAT, when applied on the three PWMs, as well as the influence of the matrix-specific compositional features in the calculation of pattern similarity scores.
Validated genes were subjected to in-silico functional enrichment analysis by Ingenuity Pathway Analysis (IPA; QIAGEN, Redwood City, CA; www.qiagen.com/ingenuity). The procedure was based on the prior calculation of the activation z-scores, which were used to infer the activation states of the predicted biological functions. An enrichment score (Fisher's exact test, p-value) was calculated to measure the overlap between observed and predicted functions 60 . A function was considered to be significantly over-represented by a set of genes if the corresponding p-value < 0.05.