Dissecting the DNA binding landscape and gene regulatory network of p63 and p53
Abstract
The transcription factor p53 is the best-known tumor suppressor, but its sibling p63 is a master regulator of epidermis development and a key oncogenic driver in squamous cell carcinomas (SCC). Despite multiple gene expression studies becoming available, the limited overlap of reported p63-dependent genes has made it difficult to decipher the p63 gene regulatory network. Particularly, 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. We use a novel iterative de novo motif search approach in conjunction with extensive ChIP-seq data to achieve a precise global distinction between p53-and p63-binding sites, recognition motifs, and potential co-factors. We integrate these data with enhancer:gene associations to predict p63 target genes and identify those that are commonly de-regulated in SCC representing candidates for prognosis and 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 (ΔNp63) governs epidermis development (Mills et al., 1999; Yang et al., 1999) and is an oncogenic driver of squamous cell carcinoma (SCC) (Cancer Genome Atlas Research Network 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, 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, for example 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) (Kouwenhoven et al., 2010; 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., 2019; Lin-Shiao et al., 2018; Qu 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 meta-analysis approach (Fischer et al., 2016a), which helped us to dissect the GRNs of the mouse and human orthologue of p53 (Fischer, 2020; 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.
Results
The p63 gene regulatory network
To identify genes commonly regulated by p63 across cell types and tissues, we employed a previously established meta-analysis approach, that has been helpful to infer core GRNs for human and mouse p53, the viral oncoprotein E7 and the cell cycle GRN (Fischer, 2019; Fischer et al., 2017; Fischer et al., 2016a; Fischer et al., 2014). From 11 genome-wide studies (Abraham et al., 2018; Bao et al., 2015; Carroll et al., 2006; Gallant-Behm et al., 2012; Karsli Uzunbas et al., 2019; Lin-Shiao et al., 2019; Saladi et al., 2017; Somerville et al., 2018; Watanabe et al., 2014; Wu et al., 2012; Zarnegar et al., 2012; Supplementary file 1), 16 publically available gene expression datasets were integrated to generate a specific p63 Expression Score (Supplementary file 2). The datasets have been obtained from knockdown (n = 12) or overexpression experiments (n = 4) of p63 in primary keratinocytes (n = 3), the keratinocyte cell line HaCaT (n = 2), the foreskin fibroblast cell line BJ (n = 1), the breast epithelial cell line MCF10A (n = 4), the squamous carcinoma cell lines H226 (n = 2), KYSE70 (n = 1), and FaDu (n = 1), as well as the pancreatic ductal adenocarcinoma cell lines BxPC3 (n = 1) and SUIT2 (n = 1) (Figure 1A and B and Supplementary file 1).
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, 2019; 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 up- and 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).
Further, we performed TF binding enrichment analysis for p63-dependently regulated genes using Enrichr (Kuleshov et al., 2016). In agreement with its established roles, we identify cell cycle gene regulators (E2F4, E2F6, SIN3A, E2F1, FOXM1, NFYA, and NFYB Fischer and Müller, 2017) and the MYC/MAX TFs as being enriched among p63-upregulated genes. Consistent with previous reports, our analysis also identifies KLF4 (Sen et al., 2012) and SMAD4 (Rodriguez 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, most datasets on p63-dependent gene expression were derived from cells in which p63 was overexpressed or depleted, without additional treatments. However, one dataset was derived from Nutlin-treated MCF10A cells (Karsli Uzunbas et al., 2019). MCF10A cells harbor wild-type p53 and DREAM targets are down-regulated in response to Nutlin treatment. Strikingly, depletion of p63 decreased their expression even further without affecting CDKN1A (p21) levels (Figure 3D), which indicates a possible cumulative effect that is independent of p53 regulatory functions.
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 (Figures 2A and 3) and their down-regulation through p53 (Fischer et al., 2016a; Fischer et al., 2016b; Schade et al., 2019; Uxa et al., 2019), we noted that negative p53 Expression Scores tend to correlate with positive p63 Expression Scores (Figure 4A). Furthermore, the results indicate that p53-induced genes (positive p53 Expression Scores) appear to be largely unaffected by p63. Consistently, expression data for 343 target genes with strong evidence for direct up-regulation by p53 (Fischer, 2017), 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 file 1) to 28 p53 ChIP-seq datasets we compiled recently (Fischer, 2019). Notably, p63 and p53 data was collected from cells with strong basal p63 expression and stimulated p53 expression, respectively. 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 ChIP-seq 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’ motifs 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 revealed a p53RE containing a 1 bp spacer (p53 secondary motif), supporting the 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 spacer-containing 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). This broader binding repertoire likely underlies p63’s capacity to engage with substantially more binding sites than does p53.
It is important to note that the vast majority (~70%) of p53- and p63-binding sites harbor full-length p53 and p63REs (Figures 6 and 7, Supplementary file 3). 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 (Figure 5—figure supplement 1). However, p63 binds many sites that are not bound by p53 (Figure 5E and Figure 5—figure supplement 2). More importantly, p53 binding is strongly constrained to canonical p53RE (Figure 5—figure supplements 1C–D and 2A–C). In contrast, p63 binding appears not to benefit from a more canonical p63RE (Figure 5—figure supplements 1E–F and 2D–F). These data suggest that sequence-specific binding is particularly important to recruit p53, while p63 only requires minimal sequence identity and could require additional co-factors 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. Consistent with the co-enrichment of AP-1 and p63 at enhancers (Lin-Shiao et al., 2018), we found that unique p63-binding sites were consistently enriched for AP-1 (bZIP) in addition to bHLH motifs (Figure 5—figure supplement 3). Using the CistromeDB toolkit (Zheng et al., 2019b), we identified TFs that significantly enrich for binding to the ‘p53+p63’, ‘p53’, and ‘p63’ sites. As expected, the analysis identified the p53 family members p53, p63, and p73 as best hits for the common sites, but only p53 and p73 for the unique p53 and only p63 and p73 for the unique p63 peak sets (Figure 5—figure supplement 4). In agreement with earlier studies, the analysis identified p300 (Katoh et al., 2019), MAF (Lopez-Pajares et al., 2015), SOX2 (Watanabe et al., 2014), BANF1 (also known as BAF) (Bao et al., 2015), and KMT2D (Lin-Shiao et al., 2018) as potential co-binders of p63; as well as TRIM28 (Doyle et al., 2010), BRD4 (Stewart et al., 2013), p300 (Lill et al., 1997), ZBTB33 (KAISO) (Koh et al., 2014), CDK9 (Claudio et al., 2006), and HEXIM1 (Lew et al., 2012) as potential co-binders of p53. Moreover, our analysis identified potential co-binders that to our knowledge have not been identified before, such as KDM1A, PRMT1, and GRHL2 for p63 and BRD9, ZNF131, and C17orf49 for p53. Importantly, these new potential co-binders appear to be unique to either p63 or p53, suggesting that they may contribute to shaping the DNA-binding landscapes that are specific to p63 and p53 (Figure 5—figure supplement 4).
Identification of direct p63 target genes
Given that p63 regulates many target genes through enhancers (Kouwenhoven et al., 2015a; Lin-Shiao et al., 2019; Lin-Shiao et al., 2018; Qu 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 E) 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) to genes with a |p63 Expression Score|| ≥ 8 (Table 1 and Figure 7—figure supplement 1). Of note, many genes are associated with proximal and enhancer p63 binding, because many proximal promoters are also identified as double-elite enhancers in the database. The 180 (138 up-regulated and 42 down-regulated) genes that passed our conservative filtering contain many genes that are known direct p63 targets, such as RAB38 (Barton et al., 2010), S100A2 (Kirschner et al., 2008; Lapi et al., 2006), HAS3 (Compagnone et al., 2017), IRF6 (Thomason et al., 2010), PTHLH (Somerville et al., 2018), GPX2 (Yan and Chen, 2006), JAG1 (Sasaki et al., 2002), MMP14 (Lodillinsky et al., 2016), NRG1 (Forster et al., 2014), and PLAC8 (Gallant-Behm et al., 2012). The identification of these well-established p63 target genes indicates the ability of our approach to identify bona fide candidates. Importantly, the integration of enhancer:gene associations enabled the identification of genes that are likely regulated by p63 through enhancers, such as IL1B, MREG, MYO5A, RRP12, SNCA, AK4, and EHD4 (Table 1 and Figure 7—figure supplement 1).
A p63/SCC 28-gene set correlates with HNSC patient survival
Out of the 180 high-probability p63 target genes 32 (28 up- and four down-regulated) are also identified as being commonly up- or down-regulated in SCCs compared to non-SCC cancers (Cancer Genome Atlas Research Network et al., 2018; Table 1). Importantly, several of the genes commonly up-regulated by p63 as well as in SCC have been identified to promote SCC growth or invasion, such as LAD1 (Abe et al., 2019), TMEM40 (Zhang et al., 2019), FGFBP1 (Czubayko et al., 1997), IL1B (Lee et al., 2015), FAT2 (Dang et al., 2016), FOSL1 (Usui et al., 2012), LPAR3 (Brusevold et al., 2014), MMP14 (Pang et al., 2016), and RASSF6 (Zheng et al., 2019a). 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 SCC (HNSC) patients from The Cancer Genome Atlas (TCGA). Notably, it is known that this cancer type frequently harbors amplified TP63 (Cancer Genome Atlas Network 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 8A). However, upon contrasting the low-expression group with all remaining samples, a significant improvement of survival was detected (p=0.024; Figure 8B). Expression of the 28-gene set correlated positively with p63 expression when p63 expression was rather low (FPKM <20), but showed a saturation and no further correlation when p63 expression was high (FPKM >20; Figure 8C and D, and Figure 8—figure supplement 1). This indicates that p63 levels influence the 28-gene set in a switch-like manner where a saturation of p63-dependent activation is quickly reached in HNSC cells. Together, these findings indicate that the genes commonly up-regulated by p63 and in SCC influence the prognosis of HNSC patients. Taken together, 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 (Fischer, 2020; Fischer, 2019), 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 (Table 1 and Figure 7—figure supplement 1).
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 (Kouwenhoven et al., 2010; 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 Figure 5—figure supplement 1A), single half-sites and half-sites separated by spacers underlie many sites that are only bound by p53 (Figures 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 (Figures 5F and 7), supporting one of the models established earlier (McDade et al., 2014). Notably, p63’s ability to bind to a greater variety of recognition motifs likely underlies the markedly greater number of p63 compared to p53-binding sites in the genome. In addition, our motif search indicates that factors bound to AP-1 (bZIP) and bHLH motifs may specifically support p63 binding (Figure 5—figure supplement 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 (Figure 5—figure supplement 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, co-factors, 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 Figure 7—figure supplement 1). Mechanistically, p63 has been shown to up-regulate target genes through its alternative TAD located at the N-terminus while the C-terminus 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 (Figures 2A and 4A). On the one hand, cell cycle genes are well-known to be down-regulated by p53 indirectly through the cyclin-dependent kinase inhibitor p21 and the cell cycle repressor complexes DREAM and RB-E2F (Fischer et al., 2016a; Fischer et al., 2016b; Schade et al., 2019; Uxa et al., 2019). On the other hand, cell cycle genes are down-regulated upon loss of p63 and this p63-dependent regulation reportedly occurs through regulating p21 signaling and the DREAM component p130 (McDade et al., 2011; Truong et al., 2006). In addition to indirect effects, we also predicted multiple cell cycle genes as direct p63 targets (Table 1). Consequently, a loss of p63 may substantially contribute to the effect of p53 in reducing cell cycle gene expression (Figure 3D). In addition of p63’s role in driving the expression of some cell cycle genes, the entire set of cell cycle genes may be subsequently up-regulated indirectly through p63’s pro-proliferative targets. 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. We provided a showcase example, where expression levels of the 28 p63 target genes that are commonly up-regulated by p63 and in SCC (Table 1) correlate significantly with poorer survival of HNSC patients (Figure 8). Thus, this 28-gene set may contain particularly promising candidates for therapeutic interventions and for the use as biomarkers.
Materials and methods
Re-analysis and integration of publicly available gene expression profiling datasets
Request a detailed protocolWe 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 Hochberg-corrected p-values (Clough and Barrett, 2016). Gene identifiers were mapped to Ensembl Gene IDs using the Ensembl annotation data (Cunningham et al., 2019). All RNA-seq datasets have been retrieved through GEO from the Sequence Read Archive (SRA) (Leinonen et al., 2011). We employed our RNA-seq analysis pipeline to obtain fold expression changes and p-values adjusted for multiple testing. Briefly, we utilized Trimmomatic (Bolger et al., 2014) v0.39 (5nt sliding window approach, mean quality cutoff 22) for read quality trimming according to inspections made from FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) v0.11.8 reports. Clipping was performed using Cutadapt v2.3 (Martin, 2011). Potential sequencing errors were detected and corrected using Rcorrector v1.0.3.1 (Song and Florea, 2015). Ribosomal RNA (rRNA) transcripts were artificially depleted by read alignment against rRNA databases through SortMeRNA v2.1 (Kopylova et al., 2012). The preprocessed data was aligned to the reference genome hg38, retrieved along with its gene annotation from Ensembl v.92 (Cunningham et al., 2019). For read alignment, we used the splice-aware mapping software segemehl (Hoffmann et al., 2014; Hoffmann et al., 2009) 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 Supplementary file 1. Given that all RNA-seq data was derived from PolyA-enriched samples, we only included Ensembl transcript types ‘protein_coding’, ‘antisense’, ‘lincRNA’ and ‘TEC’ in our analysis. Common thresholds for adj. p-value≤0.05 were applied.
Generation of the p63 expression score
Request a detailed protocolFor 19,156 genes covered by at least three datasets including a minimum of one RNA-seq 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.
Enrichment analyses
Request a detailed protocolGene set enrichment analysis (GSEA) was performed using GSEA (http://software.broadinstitute.org/gsea/) with ‘H’, ‘C2’, and ‘C6’ gene sets from MSigDB v7.0 (Subramanian et al., 2005) and custom panSCC gene sets derived from Table S1C in Cancer Genome Atlas Research Network et al., 2018. GSEA was performed on a pre-ranked list of genes that were ranked primarily by p63 Expression Score and secondarily by median log2(fold-change) to obtain unique ranks.
Enrichment of transcription factor binding to genes with high (≥8) or low (≤ −8) p63 Expression Score was identified using the results section ‘ENCODE and ChEA Consensus TFs from ChIP-X’ from Enrichr (Kuleshov et al., 2016).
Integration of publicly available p63 and p53 binding data
Request a detailed protocolPeak datasets from p63 ChIP-seq experiments were retrieved from CistromeDB (Zheng et al., 2019b; Supplementary file 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).
Motif search
Request a detailed protocolKnown p53 and p63REs were identified using the ‘known motifs’ in HOMER v4.10 with default options and -size given (Heinz et al., 2010). De novo motif discovery was performed with options -size given -len 10,15,20,25 -mis 5 S 10.
Identification of potential co-factors
Request a detailed protocolWe used the CistromeDB toolkit (Zheng et al., 2019b) to identify TFs that display ChIP-seq peaksets that are significantly similar to our ‘unique p53’, ‘unique p63’, and ‘p53+p63’ peaksets.
Survival and expression correlation analysis
Request a detailed protocolFor the 28-gene set, single-sample enrichment scores were derived from FPKM normalized gene expression values of 546 HNSC patient samples. To this end, we utilized the official GenePattern single sample gene set enrichment analysis (ssGSEA) codebase v10.0.3 (Barbie et al., 2009; https://github.com/GSEA-MSigDB/ssGSEA-gpmodule). A sample score represents the coordinately up- or down-regulated expression of all genes within one set as its signature (Barbie et al., 2009). Kaplan-Meier plots and correlation analyses were performed on TCGA time to event and event occurrence information using the R survival package v3.2–3. Following the TCGA standard for HNSC (Cancer Genome Atlas Network 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.
We retrieved read quantification data 'HTSeq - Counts' from 546 samples of the TCGA project HNSC utilizing the R package TCGAbiolinks v2.18.0 (Colaprico et al., 2016). Per sample, all read counts of the 28-gene set were merged into an artificially created metagene. Subsequently, we calculated normalized expression values per gene as fragments per kilobase million, where the length of a gene corresponds to the lengths of its exons assigned to either the canonical transcript (CCDS) or the longest transcript according to hg38 Ensembl annotation v92. TP63 FPKM values were plotted against the meta-gene FPKM value or the ssGSEA derived gene set scores (see above).
Data availability
All data generated or analysed during this study are included in the manuscript and supporting files.
-
NCBI Gene Expression OmnibusID GSE67382. A novel ATAC-seq approach reveals lineage-specific reinforcement of the open chromatin landscape via cooperation between BAF and p63.
-
NCBI Gene Expression OmnibusID GSE17611. Genome-Wide Profiling of p63 DNA-Binding Sites Identifies an Element that Regulates Gene Expression during Limb Development in the 7q21 SHFM1 Locus.
-
NCBI Gene Expression OmnibusID GSE59824. Transcription factor p63 bookmarks and regulates dynamic enhancers during epidermal differentiation.
-
NCBI Gene Expression OmnibusID GSE32061. Genome-wide analysis of p63 binding sites identifies AP-2 factors as co-regulators of epidermal differentiation.
-
NCBI Gene Expression OmnibusID GSE56674. Genome-wide characterization reveals complex interplay between TP53 and TP63 in response to genotoxic stress.
-
NCBI Gene Expression OmnibusID GSE43111. p63 attenuates epithelial to mesenchymal potential in an experimental prostate cell model.
-
NCBI Gene Expression OmnibusID GSE88859. ACTL6A Is Co-Amplified with p63 in Squamous Cell Carcinoma to Drive YAP Activation, Regenerative Proliferation, and Poor Prognosis.
-
NCBI Gene Expression OmnibusID GSE60814. Ras and TGF-β signaling enhance cancer progression by promoting the ΔNp63 transcriptional program.
-
NCBI Gene Expression OmnibusID GSE46837. SOX2 and p63 colocalize at genetic loci in squamous cell carcinomas.
-
NCBI Gene Expression OmnibusID GSE33571. Genomic profiling of a human organotypic model of AEC syndrome reveals ZNF750 as an essential downstream target of mutant TP63.
-
NCBI Gene Expression OmnibusID GSE111619. ΔNp63α Suppresses TGFB2 Expression and RHOA Activity to Drive Cell Proliferation in Squamous Cell Carcinomas.
-
NCBI Gene Expression OmnibusID GSE20286. p63 regulates an adhesion programme and cell survival in epithelial cells.
-
NCBI Gene Expression OmnibusID GSE40462. ΔNp63α represses anti-proliferative genes via H2A.Z deposition.
-
NCBI Gene Expression OmnibusID GSE111009. Control of p53-dependent transcription and enhancer activity by the p53 family member p63.
-
NCBI Gene Expression OmnibusID GSE126397. p63 establishes epithelial enhancers at critical craniofacial development genes.
-
NCBI Gene Expression OmnibusID GSE88833. ACTL6A Is Co-Amplified with p63 in Squamous Cell Carcinoma to Drive YAP Activation, Regenerative Proliferation, and Poor Prognosis.
-
NCBI Gene Expression OmnibusID GSE88832. ACTL6A Is Co-Amplified with p63 in Squamous Cell Carcinoma to Drive YAP Activation, Regenerative Proliferation, and Poor Prognosis.
-
NCBI Gene Expression OmnibusID GSE115462. TP63-Mediated Enhancer Reprogramming Drives the Squamous Subtype of Pancreatic Ductal Adenocarcinoma.
-
NCBI Gene Expression OmnibusID GSE47058. SOX2 and p63 colocalize at genetic loci in squamous cell carcinomas.
-
NCBI Gene Expression OmnibusID GSE17394. p63 regulates human keratinocyte proliferation via MYC-regulated gene network and differentiation commitment through cell adhesion-related gene network.
-
NCBI Gene Expression OmnibusID GSE33495. Genomic profiling of a human organotypic model of AEC syndrome reveals ZNF750 as an essential downstream target of mutant TP63.
References
-
Ladinin-1 is involved in cell motility and proliferation of oral squamous cell carcinoma cellsOral Surgery, Oral Medicine, Oral Pathology and Oral Radiology 128:e81–e82.https://doi.org/10.1016/j.oooo.2019.02.205
-
Trimmomatic: a flexible trimmer for illumina sequence dataBioinformatics 30:2114–2120.https://doi.org/10.1093/bioinformatics/btu170
-
p63 regulates an adhesion programme and cell survival in epithelial cellsNature Cell Biology 8:551–561.https://doi.org/10.1038/ncb1420
-
Cdk9 phosphorylates p53 on serine 392 independently of CKIIJournal of Cellular Physiology 208:602–612.https://doi.org/10.1002/jcp.20698
-
The gene expression omnibus databaseMethods in Molecular Biology 1418:93–110.https://doi.org/10.1007/978-1-4939-3578-9_5
-
TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA dataNucleic Acids Research 44:e71.https://doi.org/10.1093/nar/gkv1507
-
Solution structure and binding specificity of the p63 DNA binding domainScientific Reports 6:26707.https://doi.org/10.1038/srep26707
-
The p53-p21-DREAM-CDE/CHR pathway regulates G2/M cell cycle genesNucleic Acids Research 44:164–174.https://doi.org/10.1093/nar/gkv927
-
Census and evaluation of p53 target genesOncogene 36:3943–3956.https://doi.org/10.1038/onc.2016.502
-
Mice are not humans: the case of p53Trends in Cancer 17:30239.https://doi.org/10.1016/j.trecan.2020.08.007
-
Cell cycle transcription control: dream/MuvB and RB-E2F complexesCritical Reviews in Biochemistry and Molecular Biology 52:638–662.https://doi.org/10.1080/10409238.2017.1360836
-
Np63 represses anti-proliferative genes via H2A.Z depositionGenes & Development 26:2325–2336.https://doi.org/10.1101/gad.198069.112
-
Mechanism of TAp73 inhibition by δnp63 and structural basis of p63/p73 hetero-tetramerizationCell Death & Differentiation 23:1930–1940.https://doi.org/10.1038/cdd.2016.83
-
Fast mapping of short sequences with mismatches, insertions and deletions using index structuresPLOS Computational Biology 5:e1000502.https://doi.org/10.1371/journal.pcbi.1000502
-
Control of p53-dependent transcription and enhancer activity by the p53 family member p63Journal of Biological Chemistry 294:10720–10736.https://doi.org/10.1074/jbc.RA119.007965
-
Diversity in DNA recognition by p53 revealed by crystal structures with Hoogsteen base pairsNature Structural & Molecular Biology 17:423–429.https://doi.org/10.1038/nsmb.1800
-
Gene regulatory mechanisms orchestrated by p63 in epithelial development and related disordersBiochimica Et Biophysica Acta (BBA) - Gene Regulatory Mechanisms 1849:590–600.https://doi.org/10.1016/j.bbagrm.2015.03.003
-
Enrichr: a comprehensive gene set enrichment analysis web server 2016 updateNucleic Acids Research 44:W90–W97.https://doi.org/10.1093/nar/gkw377
-
Il‐1β promotes malignant transformation and tumor aggressiveness in oral CancerJournal of Cellular Physiology 230:875–884.https://doi.org/10.1002/jcp.24816
-
Identification of HEXIM1 as a positive regulator of p53Journal of Biological Chemistry 287:36443–36454.https://doi.org/10.1074/jbc.M112.374157
-
The Sequence Alignment/Map format and SAMtoolsBioinformatics 25:2078–2079.https://doi.org/10.1093/bioinformatics/btp352
-
KMT2D regulates p63 target enhancers to coordinate epithelial homeostasisGenes & Development 32:181–193.https://doi.org/10.1101/gad.306241.117
-
p53 and p73 display common and distinct requirements for sequence specific binding to DNANucleic Acids Research 35:340–352.https://doi.org/10.1093/nar/gkl1047
-
p63 maintains keratinocyte proliferative capacity through regulation of Skp2-p130 levelsJournal of Cell Science 124:1635–1643.https://doi.org/10.1242/jcs.084723
-
Dominant negative (ΔN) p63α induces drug resistance in hepatocellular carcinoma by interference with apoptosis signaling pathwaysBiochemical and Biophysical Research Communications 396:335–341.https://doi.org/10.1016/j.bbrc.2010.04.093
-
p53/p63/p73 isoforms: an orchestra of isoforms to harmonise cell differentiation and response to stressCell Death & Differentiation 13:962–972.https://doi.org/10.1038/sj.cdd.4401914
-
Revealing a human p53 universeNucleic Acids Research 46:8153–8167.https://doi.org/10.1093/nar/gky720
-
FGFR2 signaling underlies p63 oncogenic function in squamous cell carcinomaJournal of Clinical Investigation 123:3525–3538.https://doi.org/10.1172/JCI68899
-
Tumor suppressor p53: from engaging DNA to target gene regulationNucleic Acids Research 48:8848–8869.https://doi.org/10.1093/nar/gkaa666
-
The p53 family member genes are involved in the notch signal pathwayJournal of Biological Chemistry 277:719–724.https://doi.org/10.1074/jbc.M108080200
-
RB, p130 and p107 differentially repress G1/S and G2/M genes after p53 activationNucleic Acids Research 47:11197–11208.https://doi.org/10.1093/nar/gkz961
-
Cooperation between the transcription factors p63 and IRF6 is essential to prevent cleft palate in miceJournal of Clinical Investigation 120:1561–1569.https://doi.org/10.1172/JCI40266
-
Redox state of p63 and p73 core domains regulates sequence-specific DNA bindingBiochemical and Biophysical Research Communications 433:445–449.https://doi.org/10.1016/j.bbrc.2013.02.097
-
p63 regulates proliferation and differentiation of developmentally mature keratinocytesGenes & Development 20:3185–3197.https://doi.org/10.1101/gad.1463206
-
DREAM and RB cooperate to induce gene repression and cell-cycle arrest in response to p53 activationNucleic Acids Research 47:9087–9103.https://doi.org/10.1093/nar/gkz635
-
RSeQC: quality control of RNA-seq experimentsBioinformatics 28:2184–2185.https://doi.org/10.1093/bioinformatics/bts356
-
SOX2 and p63 colocalize at genetic loci in squamous cell carcinomasJournal of Clinical Investigation 124:1636–1645.https://doi.org/10.1172/JCI71545
-
GPX2, a direct target of p63, inhibits oxidative Stress-induced apoptosis in a p53-dependent mannerJournal of Biological Chemistry 281:7856–7862.https://doi.org/10.1074/jbc.M512655200
-
Genomic profiling of a human organotypic model of AEC syndrome reveals ZNF750 as an essential downstream target of mutant TP63The American Journal of Human Genetics 91:435–443.https://doi.org/10.1016/j.ajhg.2012.07.007
-
RASSF6-TRIM16 Axis promotes cell proliferation, migration and invasion in esophageal squamous cell carcinomaJournal of Genetics and Genomics 46:477–488.https://doi.org/10.1016/j.jgg.2019.10.004
-
Cistrome data browser: expanded datasets and new tools for gene regulatory analysisNucleic Acids Research 47:D729–D735.https://doi.org/10.1093/nar/gky1094
Article and author information
Author details
Funding
Deutsche Forschungsgemeinschaft (FI 1993/2-1)
- Martin Fischer
Bundesministerium für Bildung und Forschung (031L016D)
- Steve Hoffmann
Leibniz Association (Open Access Fund)
- Martin Fischer
- Steve Hoffmann
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Copyright
© 2020, Riege et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 2,551
- views
-
- 393
- downloads
-
- 34
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Chromosomes and Gene Expression
- Neuroscience
Pathogenic variants in subunits of RNA polymerase (Pol) III cause a spectrum of Polr3-related neurodegenerative diseases including 4H leukodystrophy. Disease onset occurs from infancy to early adulthood and is associated with a variable range and severity of neurological and non-neurological features. The molecular basis of Polr3-related disease pathogenesis is unknown. We developed a postnatal whole-body mouse model expressing pathogenic Polr3a mutations to examine the molecular mechanisms by which reduced Pol III transcription results primarily in central nervous system phenotypes. Polr3a mutant mice exhibit behavioral deficits, cerebral pathology and exocrine pancreatic atrophy. Transcriptome and immunohistochemistry analyses of cerebra during disease progression show a reduction in most Pol III transcripts, induction of innate immune and integrated stress responses and cell-type-specific gene expression changes reflecting neuron and oligodendrocyte loss and microglial activation. Earlier in the disease when integrated stress and innate immune responses are minimally induced, mature tRNA sequencing revealed a global reduction in tRNA levels and an altered tRNA profile but no changes in other Pol III transcripts. Thus, changes in the size and/or composition of the tRNA pool have a causal role in disease initiation. Our findings reveal different tissue- and brain region-specific sensitivities to a defect in Pol III transcription.
-
- Biochemistry and Chemical Biology
- Chromosomes and Gene Expression
The mRNA 5'-cap structure removal by the decapping enzyme DCP2 is a critical step in gene regulation. While DCP2 is the catalytic subunit in the decapping complex, its activity is strongly enhanced by multiple factors, particularly DCP1, which is the major activator in yeast. However, the precise role of DCP1 in metazoans has yet to be fully elucidated. Moreover, in humans, the specific biological functions of the two DCP1 paralogs, DCP1a and DCP1b, remain largely unknown. To investigate the role of human DCP1, we generated cell lines that were deficient in DCP1a, DCP1b, or both to evaluate the importance of DCP1 in the decapping machinery. Our results highlight the importance of human DCP1 in decapping process and show that the EVH1 domain of DCP1 enhances the mRNA-binding affinity of DCP2. Transcriptome and metabolome analyses outline the distinct functions of DCP1a and DCP1b in human cells, regulating specific endogenous mRNA targets and biological processes. Overall, our findings provide insights into the molecular mechanism of human DCP1 in mRNA decapping and shed light on the distinct functions of its paralogs.