Dissecting the DNA binding landscape and gene regulatory network of p63 and p53

The transcription factor (TF) p53 is the best-known tumor suppressor, but its ancient sibling p63 (ΔNp63) is a master regulator of epidermis development and a key oncogenic driver in squamous cell carcinomas (SCC). Despite multiple gene expression studies becoming available in recent years, the limited overlap of reported p63-dependent genes has made it difficult to decipher the p63 gene regulatory network (GRN). In particular, analyses of p63 response elements differed substantially among the studies. To address this intricate data situation, we provide an integrated resource that enables assessing the p63-dependent regulation of any human gene of interest. Here, we use a novel iterative de novo motif search approach in conjunction with extensive publicly available ChIP-seq data to achieve a precise global distinction between p53 and p63 binding sites, recognition motifs, and potential co-factors. We integrate all these data with enhancer:gene associations to predict p63 target genes and identify those that are commonly de-regulated in SCC and, thus, may represent candidates for therapeutic interventions.


Introduction
In contrast to the tumor suppressor p53 with its extensive set of target genes controlling the cell cycle and apoptosis (Fischer, 2017;Sammons et al., 2020), its phylogenetically ancient sibling p63 governs epidermis development (Mills et al., 1999;Yang et al., 1999) and is an oncogenic driver of squamous cell carcinoma (SCC) (Campbell et al., 2018;Gatti et al., 2019) that is overexpressed or amplified in SCCs, which depend on its expression (Ramsey et al., 2013). Together with p73, p63 and p53 form the p53 transcription factor (TF) family that shares a highly conserved DNA binding domain (DBD) through which they bind to very similar DNA recognition motifs. The mechanisms that enable these sibling TFs to shape their unique gene regulatory network (GRN) leading to the different phenotypic control, however, remain poorly understood.
The TP53 and TP63 genes encode for two major isoform groups that are controlled by distinct promoters leading to transcripts differing in their N-terminus (Murray-Zmijewski et al., 2006). In the case of TP53, the longest isoform, p53α, is ubiquitously expressed while the alternative intronic promoter has little activity across virtually all tissues. Conversely, the usage of the two TP63 promoters is highly cell type-dependent. For instance, the long isoform TAp63 is predominantly expressed in germ cells, while the smaller transcript, ΔNp63, is most copious in stratifying epithelia (Sethi et al., 2015). Similar to p53, alternative splicing leads to α, β, and γ protein isoforms that differ in their C-terminus (Murray-Zmijewski et al., 2006). While both TAp63 and ΔNp63 may bind to DNA through a specific binding domain, ΔNp63 lacks the canonical N-terminal transactivation domain (TAD) (Yang et al., 1998) and has long been thought to be a dominant-negative regulator of other p53 family members or its own isoforms (Gebel et al., 2016;Yang et al., 1998). However, ΔNp63 has also been shown to harbor alternative TADs, that endow transactivation activity (Helton et al., 2006;King et al., 2003;Yang et al., 2006). Notably, the many ΔNp63 binding sites are associated with enhancer regions, where ΔNp63 has been proposed to "bookmark" genes that are expressed in stratifying epithelia (Karsli Uzunbas et al., 2019;Kouwenhoven et al., 2015a;Lin-Shiao et al., 2019;Qu et al., 2018;Somerville et al., 2018). Here, we focus on the most widely expressed isoforms p53α (hereafter p53) and ΔNp63 (hereafter p63).
The p53 TF family shares many binding sites, but all three family members have been shown to bind to substantial sets of unique target genes (Lin et al., 2009;McDade et al., 2014).
Indeed, there are differences in the DBDs, e.g. regarding thermostability, hydrophobic potentials (Enthart et al., 2016), zinc-coordination (Lokshin et al., 2007), and redox sensitivity (Tichý et al., 2013). In addition, the different C-terminal domains (CTD) of p53 family members may also affect their DNA binding specificity (Sauer et al., 2008). p53 binds to a canonical 20 bp response element (RE) made of two decameric half-sites that both contain the sequence RRRCWWGYYY (R = A/G; W= A/T; Y = C/T). p53 has also been shown to bind to decameric half-sites separated by spacers or to single half-sites (Kitayner et al., 2010;Menendez et al., 2013;Vyas et al., 2017). Results from systematic evolution of ligands by exponential enrichment (SELEX) (Ortt and Sinha, 2006;Perez et al., 2007) and high-throughput analyses of chromatin immunoprecipitation (ChIP) McDade et al., 2012;Yang et al., 2006) yielded p63 binding motifs with high similarity to the p53RE but still showed some unique characteristics. These unique characteristics identified for p63REs, however, differed substantially between the studies.
While multiple genome-wide p63 gene expression datasets became available in recent years, our understanding of the p63 GRN remains incomplete. This is in part due to the limited overlap of the p63-dependent genes identified in individual studies (Kouwenhoven et al., 2015b). Also, the frequent binding of p63 to enhancers (Kouwenhoven et al., 2015a;Lin-Shiao et al., 2018Qu et al., 2018;Somerville et al., 2018) and the difficulty to associate such enhancers with target gene regulation adds another level of complexity to the quest of describing the GRN. To overcome these limitations, we utilize a recently developed metaanalysis approach (Fischer et al., 2016a), which helped us to dissect the GRNs of the mouse and human orthologue of p53 (Fischer, 2019). The analysis rests upon a ranking of potential p63 target genes based on the number of datasets supporting a p63-dependent regulation. In addition, we utilize the wealth of recent p63 and p53 ChIP-seq studies to establish a more precise global distinction between p53 and p63 binding sites and their underlying REs. This approach could serve as a blueprint to distinguish binding site specificities of TF siblings.
Further integration of gene expression studies with the binding data and enhancer:gene associations enables us to predict high-probability direct p63 target genes.
To illustrate the utility of our approach, we selected 30 genes from various p63 Expression Score groups reflecting commonly up-and down-regulated ones ( Figure 1C). We noted lower consistency across the data on p63-dependent gene regulation as compared to previous meta-analyses on human and mouse p53 Fischer et al., 2016a). In contrast to the recent investigations, data integrated here are based on a higher number of experiments in primary cells and a comparably lower number of replicates. Thus, the reduced consistency may also reflect the higher variance as opposed to data from more homogenous cell lines. Furthermore, p63-depleted cells are less viable, and the global decrease in mRNA levels may confound effects. Despite this, our approach identified genes that are commonly altered by p63.
We next performed gene set enrichment analysis (GSEA) for p63-dependently regulated genes using MSigDB gene sets (Subramanian et al., 2005). In agreement with the function of p63 as an essential proliferation factor (McDade et al., 2011;Senoo et al., 2007;Truong et al., 2006), epidermal development regulator (Mills et al., 1999;Yang et al., 1999), and MYC network activator (Wu et al., 2012), we find that genes commonly up-regulated by p63 significantly enrich gene sets associated with cell cycle, epidermis development, and MYC targets ( Figure 2A). In line with previous reports (Mehta et al., 2018), genes down-regulated by p63 enrich gene sets connected with interferon response ( Figure 2B). Corroborating the role of p63 in mammary stem cell activity (Chakrabarti et al., 2014) and SCC growth (Ramsey et al., 2013), we find that p63 up-and down-regulated genes enrich respective gene sets upand down-regulated in mammary stem cells ( Figure 2C) and across SCCs ( Figure 2D). In addition to pathways that have been linked to p63 earlier, we find that p63 up-regulated genes enrich for mTORC1 signaling genes and p63 down-regulated genes enrich for gene sets associated with oxidative phosphorylation and aerobic respiration ( Figure 2E).
Consistent with previous reports, our analysis also identifies KLF4  and SMAD4 (Calleja et al., 2016) as potential mediators of p63-dependent gene regulation. In addition, our analysis reveals that androgen receptor (AR), its co-factor ZMIZ1, as well as SP1, FLI1, and NANOG are novel candidates for mediating the p63-dependent up-regulation of multiple genes. Surprisingly, our analysis identified only SOX2 as a frequent binder of genes down-regulated by p63 ( Figure 3A). Consistent with the strong association of p63 up-regulated genes with the cell cycle ( Figure 2A) and with cell cycle regulators ( Figure 3A), we find that p63 up-regulated genes enrich DREAM (dimerization partner, RB-like, E2F, and multi-vulval class B) and E2F target genes ( Figure 3B), and DREAM target genes appear to be modestly but consistently down-regulated when p63 is lost ( Figure 3C). Notably, when DREAM targets were already down-regulated through Nutlin treatment in MCF10A cells, loss of p63 decreased their expression even further ( Figure 3C) (Karsli Uzunbas et al., 2019), which indicates a possible cumulative effect that is independent of p53 regulatory functions as has been suggested previously (McDade et al., 2014).
Together, the meta-analysis approach overcomes the limitations of individual studies and identifies target genes supported by multiple datasets. The extensive and integrated resource on p63-regulated genes enables researchers to compare their results quickly and to identify the most promising targets.

p63 and p53 regulate largely distinct gene sets
Given that p63 and p53 share a significant number of binding sites and thus potential target genes, we next compared the p63 Expression Score to the previously established p53 Expression Score (Fischer et al., 2016a). In agreement with the up-regulation of cell cycle genes and DREAM targets through p63 (Figure 2A  appear to be largely unaffected by p63. Consistently, expression data for 346 target genes with strong evidence for direct up-regulation by p53, do not show consistent expression changes upon knockdown or induction of p63 ( Figure 4B). Together, these results indicate that basal expression of the majority of p53 target genes is not affected by p63.

Common and distinct properties of p63 and p53 DNA binding
To identify shared p63 and p53 bound sites, we compared the 20 p63 ChIP-seq datasets (Supplementary Table 1) to 28 p53 ChIP-seq datasets we compiled recently . While the majority of all p53 ChIP-seq peaks occurs in only one of the experiments, more than half of the p63 peaks are present in two or more datasets ( Figure 5A and B). Even though we were able to integrate substantially more p53 datasets, the number of identified p63 binding sites was still higher ( Figure 5C). This indicates that p63 occupies many more binding sites as compared to p53. Importantly, when more datasets agree on p53 and p63 binding sites, these sequences are more likely to harbor a canonical p53 and p63RE, facilitating the motif discovery by tools such as HOMER (Heinz et al., 2010) and enriching bona fide binding sites ( Figure 5D). Earlier meta-analyses employed a similar strategy (Fischer et al., 2016a;Nguyen et al., 2018;Verfaillie et al., 2016). To dissect the binding preferences of p63 and p53, we generated three distinct peak sets ( Figure 5E). The 'p53+p63' set contained all binding sites with evidence in at least five p63 and five p53 ChIP-seq datasets. The 'p53 unique' (hereafter 'p53') set contained all binding sites that were supported by at least five p53 ChIPseq datasets but not a single p63 dataset. We also generated a 'p63 unique' (hereafter 'p63') set vice versa.
We employed an iterative de novo motif search using HOMER to identify frequent binding site motifs. After each round, we removed all peaks harboring the best motif and repeated the search. We identified similar yet distinct binding motifs for the three groups ( Figure 5F). Comparison of the primary 'p53+p63', 'p53' and 'p63' suggests that p63 binding sites display a highly conserved C, G, C and G at positions 4, 7, 14 and 17, respectively. The second round reveals a p53RE containing a 1bp spacer (p53 secondary motif), supporting a previously published model that p53 can bind to spacer-containing p53REs (Vyas et al., 2017).
The results further indicate that p53 can bind to a single half-site (p53 tertiary motif) and that this single half-site is more constrained at positions 5 and 6 as well as the flanking regions than half-sites in the canonical p53RE (e.g. primary p53+p63 and p53 motifs). Of note, these single half-sites may also include p53REs with spacers longer than 1 bp that are not detected separately because of their very low abundance. Sole half-sites together with spacercontaining p53REs underlie only ~5 % of p53 bound sites ( Figure 6). Furthermore, p53 and p63 appear to be able to bind to three-quarter sites (secondary and quaternary p53+p63 and p63 motifs) while p63 can generally bind to a broader spectrum of sequences as compared to p53 ( Figure 5F).
It is important to note that the vast majority (~70 %) of p53 and p63 binding sites harbor full-length p53 and p63REs ( Figure 6 and 7, Supplementary Table S3). There is a good correlation between p53 and p63 binding site occupation, and most sites commonly bound by p53 are also frequently bound by p63 (Supplementary Figure 1A). However, p63 binds many sites that are not bound by p53 (Supplementary Figure 1B). More importantly, p53 binding is strongly constrained to canonical p53RE . In contrast, p63 binding appears not to benefit from a more canonical p63RE . These data suggest that sequence-specific binding is particularly important to recruit p53, while p63 only requires minimal sequence identity and could require additional cofactors to bind and ultimately regulate its target genes. Therefore, we also searched for potential cooperating TFs that may be co-enriched at p53 and p63 binding sites. Consistent with earlier analyses (Verfaillie et al., 2016), no additional motif was substantially enriched in the vicinity of 'p53' or 'p53+p63' binding sites.

Identification of direct p63 target genes
Given that p63 regulates many target genes through enhancers (Kouwenhoven et al., 2015a;Lin-Shiao et al., 2018Qu et al., 2018;Somerville et al., 2018), straight forward integration of differential gene regulation data and p63 binding data based on proximity binding to a gene's TSS is unlikely to capture all direct p63 target genes. To resolve this issue, we integrated the p63 binding data and the p63 Expression Score based on enhancer:gene association information (Fishilevich et al., 2017) in addition to proximity binding to TSSs to predict direct p63 target genes. Given the large number of p63 binding sites identified ( Figure   5C and 5E) and the high variance in p63-dependent gene regulation ( Figure 1B), we employed conservative thresholds to identify high-probability target genes of p63. We only used p63 binding sites supported by at least half of the datasets (≥10) that are linked through TSS proximity (within 5 kb) or double-elite enhancer:gene associations (Fishilevich et al., 2017)  Therefore, we asked whether the set of 28 up-regulated direct p63 targets correlates with patient survival. To this end, we employed data of head and neck squamous cell carcinoma (HNSC) patients from The Cancer Genome Atlas (TCGA). Notably, it is known that this cancer type frequently harbors amplified TP63 (Lawrence et al., 2015). We find that expression levels of our gene set indeed correlate significantly negatively with HNSC patient survival (COX likelihood ratio test p=0.032). To determine whether expression levels of the set have an influence on the survival of HNSC patients, we subdivided the samples according to the average expression levels into four equally sized groups (low, low-med, med-high, high). While the sample group with low expression had the most favorable prognosis, the null hypothesis could not be rejected in the direct comparison with patients with high average expression levels (p=0.090; Figure 9a). However, upon contrasting the low-expression group with all remaining samples, a significant improvement of survival was detected (p=0.024; Figure 9b). Together, these findings indicate that the genes commonly up-regulated by p63 and in SCC influence the prognosis of HNSC patients. Thus, this finding calls for a more detailed assessment of ubiquitous p63/SCC genes as biomarkers in the future.

Discussion
Although p63 (ΔNp63) is known as master regulator in epidermis development and more recently emerged as a key oncogenic factor in SCC, a comprehensive assessment of the GRN commonly controlled by p63 and its comparison to the GRN commonly controlled by the closely related tumor suppressor p53 has been missing. An increasing number of available high-throughput datasets enabled us to generate ranked lists of p63-regulated genes and p63 bound DNA sites that together reveal high-probability direct p63 target genes regulated by p63 across cells of multiple origins. Because p63 target genes, very much like p53 target genes , differ substantially between mouse and human (Sethi et al., 2017), many p63 target genes initially described in mouse could not be confirmed to be p63-regulated in this study using human data. Given that p63 binding sites are frequently associated with enhancer regions and enhancer identity, we have integrated enhancer:gene associations to identify target genes that are regulated by p63 through direct binding to associated enhancers. This approach enabled the identification of novel direct target genes that are missed by standard analyses that employ only TSS proximity (Figure 8).
Given the similarity between their DBDs, it has been a long-standing question how p53 and p63 bind to distinct sites in the genome and how these sites differ from another. Several studies found differences in the biochemical properties of p53 and p63 that could affect their DNA binding specificity (Enthart et al., 2016;Lokshin et al., 2007;Sauer et al., 2008;Tichý et al., 2013). Various studies aimed to identify the precise p63 recognition motif and its difference from the p53RE using either SELEX (Ortt and Sinha, 2006;Perez et al., 2007) or ChIP-seq data McDade et al., 2014;Yang et al., 2006), yet these studies reported different features as being unique for p63 compared to p53 DNA recognition. By combining multiple ChIP-seq datasets we have contributed here to better distinguish between sites commonly bound by p53 and p63 across cell types and sites that are unique to p53 or p63 ( Figure 5E). Most importantly, our results could explain why a substantial fraction of DNA sites is occupied exclusively by p53 or p63. While most sites bound by p53 are also commonly occupied by p63 ( Figure 5E and Supplementary Figure 1A), single half-sites and half-sites separated by spacers underlie many sites that are only bound by p53 ( Figure 5F and 6), supporting earlier findings whereby p53 can be recruited through spacer-containing motifs (Vyas et al., 2017). However, while spacers reportedly have been identified in fifty percent of 200 analyzed p53REs (Vyas et al., 2017), our genome-wide quantification of motifs underlying 7705 high confidence p53 peaks based on an unbiased motif search using HOMER revealed that only 1.1 to 5.1% of the p53 peaks contain p53REs with 1 bp spacers or half sites that are possibly separated by longer spacers (Figure 6). Mechanistically, our results imply that relying on the CWWG core motif and the flanking regions may enable p53 to bind to those sites. In contrast, the two CNNG core motifs that underlie p63, but not p53REs, offer an explanation why a substantial fraction of DNA sites is bound exclusively by p63 ( Figure 5F and 7), supporting one of the models established earlier (McDade et al., 2014). In addition, our motif search indicates that factors bound to AP-1 (bZIP) and bHLH motifs may specifically support p63 binding (Supplementary Figure 3), and transcription factor enrichment analysis identified the bZIP TF MAF, the TF GRHL2, the chromatin remodeler BANF1, the histone methyltransferase PRMT1, and the ZNF750/KDM1A/KLF4 complex, which was previously shown to operate downstream of p63 (Boxer et al., 2014), as potential co-binders that could help to facilitate p63 binding to certain genomic loci (Supplementary Figure 4). Considering its pioneer role, p63 could vice versa enable the binding of these TFs to the respective loci. Given that p63 and p73 form stable heterotetramers (Gebel et al., 2016), p73 may possess binding specificities that are highly similar to those identified for p63. Our results indicate that our approach could serve as a blueprint to distinguish DNA recognition motifs, binding sites, cofactors, and target genes of TF siblings more precisely. Our iterative de novo search algorithm enabled the identification of spacer-containing p53REs, indicating that our approach uncovers second-tier TF binding motifs invisible to standard approaches. Moreover, the results provide insights to the p63 DNA binding repertoire in unprecedented depth ( Figure 5F). Consistent with results from an earlier genome-wide study (Yang et al., 2006), our findings imply that p63 is more frequently involved in a direct up-regulation as opposed to a direct down-regulation of target genes ( Figure 3A and 8). Mechanistically, p63 has been shown to up-regulate target genes through its alternative TAD located at the N-terminus while the Cterminus is important for down-regulation (Helton et al., 2006). Exogenous expression of different isoforms of p53 family members and their antagonistic effects on target gene promoters in luciferase reporter assays suggested a model whereby p63 exhibits a dominant negative effect on other p53 family members (Mundt et al., 2010;Westfall et al., 2003;Yang et al., 1998). Inconsistent with its reputation as dominant negative regulator of p53, however, genome-wide studies showed that the groups of p63-regulated genes and p53-regulated genes show only very little overlap (Gallant-Behm et al., 2012). A recent analysis of DNA sites bound and of genes regulated by p53 and p63 revealed that p63 is more likely to support than to inhibit p53 activity (Karsli Uzunbas et al., 2019). Our analysis further supports the notion that p63 does not commonly interfere with target gene up-regulation by p53 but that except for cell cycle genes they regulate largely distinct gene sets (Figure 4).
We identify several candidate TFs that may operate downstream of p63 and that may serve as transitional nodes in the p63 GRN. In addition to known mediators of p63-dependent gene regulation, such as MYC and KLF4, we identify AR and its co-factor ZMIZ1, SP1, FLI1, and NANOG as novel candidate nodes in the p63 GRN ( Figure 3A). In agreement with the tumor suppressor role of p53 and the oncogenic role of p63, we find that cell cycle genes are antagonistically regulated by p53 and p63 (Figure 2A Truong et al., 2006). In addition to indirect effects, we also identified multiple cell cycle genes as high-probability direct p63 targets (Figure 8). Thus, p63 appears to directly drive the expression of some cell cycle genes and to indirectly up-regulate the whole set of cell cycle genes through affecting the function of cell cycle regulators. While the up-regulation of cell cycle genes occurs in most cancers (Whitfield et al., 2006), we find that p63 additionally regulates genes that are specifically altered across SCCs ( Figure 2D). These results underscore the critical role of p63 and its target genes in determining the transcriptional profile of SCC. An example of a p63 target in SCC is NRG1, which can be inhibited to block SCC proliferation and tumor growth (Hegde et al., 2019). The resource of genes commonly regulated by p63 provided here may help to identify targets that can be exploited therapeutically. Expression levels of the 28 p63 target genes that are commonly up-regulated by p63 and in SCC (Figure 8) correlate significantly with poorer survival of HNSC patients ( Figure 9). This 28-gene set may contain particularly promising candidates for therapeutic interventions and for the use as biomarkers.

Re-analysis and integration of publicly available gene expression profiling datasets
We re-analyzed publicly available p63-dependent gene expression profiling datasets.
As a first quality requirement, we only included datasets for re-analysis that contained at least two biological replicates for the treatment as well as for the control condition. All microarray datasets were available at a pre-processed stage at the Gene Expression Omnibus (GEO) and we re-analyzed these datasets with GEO2R to obtain fold expression changes and Benjamini . For read alignment, we used the splice-aware mapping software segemehl (Hoffmann et al., 2009(Hoffmann et al., , 2014 v0.3.4 with adjusted accuracy (95%). Mappings were filtered by Samtools v1.9 (Li et al., 2009) for uniqueness and properly aligned mate pairs. Read quantification was performed on exon level using featureCounts v1.6.5 (Liao et al., 2014), parametrized according to the strand specificity inferred through RSeQC v3.0.0 (Wang et al., 2012). Differential gene expression and its statistical significance was identified using DESeq2 v1.20.0 (Love et al., 2014). Information on the samples that were compared for each dataset is included in and 'TEC' in our analysis. Common thresholds for adj. p-value ≤ 0.05 were applied.

Generation of the p63 Expression Score
For 19,156 genes covered by at least three datasets including a minimum of one RNAseq dataset, a p63 Expression Score was calculated as the number of datasets that find the gene to be significantly up-regulated minus the number of datasets that find the gene to be significantly down-regulated in dependence on p63. This meta-analysis resulted in 27 p63 Expression Score gene groups because no gene was identified as up-regulated in all 16 or 15 datasets or down-regulated in all 16, 15, 14 or 13 datasets.

Integration of publicly available p63 and p53 binding data
Peak datasets from p63 ChIP-seq experiments were retrieved from CistromeDB (Zheng et al., 2019b) (Supplementary Table 1). When replicate experiments were available, all peaks were used that have been identified in at least two replicates. A similar collection of p53 peak datasets has been described previously (Fischer, 2019). To intersect multiple peak files Bedtools 'multiinter' was used and to identify overlapping and non-overlapping peaks Bedtools 'intersect' was employed (Quinlan and Hall, 2010).

Survival analysis
For the 28-gene set, single-sample enrichment scores were derived from FPKM  (Lawrence et al., 2015), survival analyses were right-censored at 60 months (1800 days) to avoid non-cancer-related events. The Cox proportional hazards model was used to investigate the association of patient survival time and the combined expression levels of the 28-gene set. Subsequently, we subdivided the expression scores into four equally sized categorical groups (high, med-high, med-low, low). The rates of occurrence of events over time were compared between these groups using the fitted COX PH model.

Author Contributions
MF conceived and supervised the study. MF, SH, KR, HK, and AS designed the analyses. KR, MF, and HK performed the analyses. MF, SH, and SSM interpreted the data.
MF wrote the initial manuscript draft. All authors edited and approved the manuscript.

Conflicts of interest
The authors declare no conflict of interest.

Reference List
Abe      peak sets. The first round of motif search identified the 'primary' motif in each peak set. Using an iterative approach, all peaks that contained the 'primary' motif were removed and the de novo motif search was repeated. This iterative approach was followed until no more p53/p63like motif was identified.

Figure 6
The DNA binding landscape of p53. DNA sites occupied by p53 in at least 5 datasets were searched iterative with the motifs identified by our iterative de novo search ( Figure 5F).
We searched first for the primary 'p53+p63' motif and among all remaining sites for the primary 'p53' motif. All other 'p53+p63' and 'p53' motifs were searched subsequently.

Figure 7
The DNA binding landscape of p63. DNA sites occupied by p63 in at least 5 datasets were searched iterative with the motifs identified by our iterative de novo search ( Figure 5F).
We searched first for the primary 'p53+p63' motif and among all remaining sites for the primary 'p63' motif. All other 'p53+p63' and 'p63' motifs were searched subsequently (Supplementary   Table S3).

Figure 8
High-probability direct p63 target genes. Genes identified as significantly up-or down-regulated in at least the half of all datasets (|p63 Expression Score| ≥ 8) that are linked to p63 binding sites supported by at least half of all datasets (≥ 10) through binding within 5 kb from their TSS or through double-elite enhancer:gene associations (Fishilevich et al., 2017).
Using these thresholds we identified 138 and 42 high-probability candidates as directly up-and down-regulated by p63, respectively. Gene names marked in red are also up-or downregulated across SCCs (Campbell et al., 2018). Correlation between HOMER motif score for primary and secondary 'p53+p63' motifs and dataset support for (C and D) p53 binding or (E and F) p63 binding.

Supplementary Figure 2
Correlation between HOMER motif score for primary, secondary and tertiary (A to C) 'p53' motifs or (D to F) 'p63' motifs and dataset support for (A to C) p53 binding or (D to F) p63 binding.

Supplementary Figure 4
Top 20 TFs with ChIP-seq peak sets similar to (A) the common p53+p63 sites, (B) the unique p53 sites, and (C) the unique p63 sites ( Figure 5E) as identified using CistromeDB toolkit. Of note, some TP53 ChIP-seq datasets are wrongly labeled "T" in the database.