Genome-wide analysis of p63 binding sites identifies AP-2 factors as co-regulators of epidermal differentiation.

The p63 transcription factor (TP63) is critical in development, growth and differentiation of stratifying epithelia. This is highlighted by the severity of congenital abnormalities caused by TP63 mutations in humans, the dramatic phenotypes in knockout mice and de-regulation of TP63 expression in neoplasia altering the tumour suppressive roles of the TP53 family. In order to define the normal role played by TP63 and provide the basis for better understanding how this network is perturbed in disease, we used chromatin immunoprecipitation combined with massively parallel sequencing (ChIP-seq) to identify >7500 high-confidence TP63-binding regions across the entire genome, in primary human neonatal foreskin keratinocytes (HFKs). Using integrative strategies, we demonstrate that only a subset of these sites are bound by TP53 in response to DNA damage. We identify a role for TP63 in transcriptional regulation of multiple genes genetically linked to cleft palate and identify AP-2alpha (TFAP2A) as a co-regulator of a subset of these genes. We further demonstrate that AP-2gamma (TFAP2C) can bind a subset of these regions and that acute depletion of either TFAP2A or TFAP2C alone is sufficient to reduce terminal differentiation of organotypic epidermal skin equivalents, indicating overlapping physiological functions with TP63.

INTRODUCTION

The p53 (TP53) family member p63 (TP63) was discovered in the quest for novel TP53 like tumour suppressor genes (1,2). Initial studies in knockout mice revealed that TP63 did not represent a canonical tumour suppressor, rather it played a critical role in development and differentiation of stratifying epithelia (2,3). TP63 knockout mice die perinatally as a result of dehydration, due to an almost complete absence of mature epidermis (2,3). They also exhibit severe limb and craniofacial defects, phenotypes that are mirrored in human diseases caused by TP63 mutation (4,5). These conditions follow an autosomal dominant pattern of inheritance and have been identified in five syndromic and two non-syndromic conditions (4), characterized by ectodermal dysplasia, split-hand/foot malformation (SHFM), orofacial-clefting (OFC) including cleft lip/cleft palate (CL/P) (5).

The precise role played by TP63 in normal growth and development of stratifying epithelia and how this is perturbed in these congenital abnormalities remains unclear, likely the result of different mutations affecting different subsets of TP63 isoforms (5). The TP63 locus encodes at least six isoforms with differing transcriptional activities, as the result of alternative promoter usage (TA and deltaN) and C-terminal splicing (alpha, beta and gamma) (1,6). The deltaN isoforms were initially shown to act in a dominant negative manner leading to inhibition of transcriptional activities of TP53, TP73 and TA-TP63 isoforms (7) but more recently have been shown to be capable of inducing a subset of target genes via a second transactivation domain (8). The alpha variants of the TA and deltaNTP63 isoforms possess a C-terminal extension containing a sterile alpha motif (SAM), a transcriptional inhibitory domain (TID) and the deltaNTP63alpha isoform is predominant isoform expressed in basal keratinocytes and stratifying epithelia (1).

The presence of this TID allows deltaNTP63alpha to repress TP53-mediated activation of growth suppressive and pro-apoptotic genes such as p21 (CDKN1A), FAS and NOXA in response to DNA damage (9). Recent data indicate that TP63 plays dual TP53-dependent and -independent roles in the epidermis, maintaining proliferative capacity by opposing TP53-mediated activation of anti-proliferative signals such as CDKN1A, while promoting transcriptional induction of genes required for proliferation and terminal differentiation (10–16). However, complementation experiments in TP63 knockout mice indicates that restoration of individual or multiple TP63 isoforms is insufficient to restore epidermal differentiation (13), likely indicating that the full complement of TP63 isoforms may be required to orchestrate this complex process.

This complex interplay between TP63 isoforms and TP53 family members is reflected in a dual role in tumorigenesis, where TP63 is suggested to play both oncogenic and tumour suppressive roles [reviewed in (6,17)] dependent on the cell type and isoforms expressed. TP63 and particularly the deltaNTP63alpha isoform are frequently overexpressed in squamous cell carcinomas [reviewed in (6,17)], where elevated levels have been shown to repress activation of growth repressive and apoptotic TP53 family target genes to maintain tumour growth. Furthermore, there is increasing evidence to suggest that this TP63 function, in particular, the tumour suppressive TAP63alpha isoform, may be altered by mutant gain of function TP53 during oncogenic progression, promoting invasion and metastasis (18,19). Definition of the normal TP63-regulated network is therefore critical in better understanding its role in development and disease.

Like the TP53 family, the AP-2 family of transcription factors has been shown to play important roles in development, differentiation and tumorigenesis (20). AP-2 alpha (TFAP2A) is one of five AP-2 family members (alpha, beta, gamma, delta and epsilon) all of which have the ability to bind the same consensus sequence (21). Like TP63 mutations, heterozygous mutation of AP-2 alpha (TFAP2A) is associated with syndromic cleft palate in branchiooculofacial syndrome, which is additionally characterized by cutaneous/ocular anomalies, facial abnormalities (22).

The advent of next-generation sequencing and its application to chromatin immunoprecipitation (ChIP-seq) [reviewed in (23)] have enabled us to measure TP63 interactions on a genome-wide basis. Here, we use ChIP-seq to catalogue genome-wide TP63 binding regions in primary human neonatal foreskin keratinocytes and identified 7574 highly conserved TP63-binding sites. Furthermore, we delineate the differences between TP63- and p53-regulated genes and pathways and identify a role for TP63 in transcriptional regulation of multiple genes involved in palatal fusion. In addition, using in silico strategies, we identify AP-2alpha (TFAP2A), a transcription factor also associated with facial clefting (22), as a co-regulator of TP63 target genes. Also we demonstrate that TFAP2A and its family member AP-2gamma (TFAP2C) can co-associate with overlapping subsets of TP63 binding sites resulting in modulation of TP63 transcriptional activity. Finally, we demonstrate that TFAP2A and TFAP2C depletion as with TP63 knockdown results in alterations in terminal differentiation of human foreskin keratinocytes in organotypic raft cultures, underlining the complex interplay between TP63 and AP-2 factors.

MATERIALS AND METHODS

Cell culture

Primary human Foreskin keratinocytes (HFKs) were isolated from human neonatal foreskins as described previously (24) and cultured in Epilife supplemented with human keratinocyte growth supplement (Gibco). HFKs were transfected using FuGeneHD (Roche) with siRNA-targeting TFAP2A (20), TP63 or a scrambled control (10) and harvested 48 h post-transfection for protein or RNA analysis. siRNA-transfected cells were used to generate organotypic raft cultures as previously described (10) and allowed to differentiate at the air–liquid interface for 5 days before fixation in 4% paraformaldehyde. Fixed raft cultures were paraffin embedded and sections stained for haemotoxylin and eosin in the Centre for Cancer Research and Cell Biology (CCRCB) core tissue facility.

Chromatin immunoprecipitation

Cells were cross-linked by addition of formaldehyde to a final concentration of 1% and incubated for 10 min at room temperature. Cross-linking was stopped by addition of Glycine to a final concentration of 0.125 M washed and washed twice with ice-cold PBS. Cross-linked cells were re-suspended in cell collection buffer (100 mM Tris–HCl, pH 9.4, 10 mM DTT), incubated on ice for 10 min and sequentially re-suspended in PBS, NCP1 (10 mM EDTA, 0.5 mM EGTA, 10 mM HEPES, pH 6.5, 0.25% Triton X100), NCP2 (1 mM EDTA, 0.5 mM EGTA, 10 mM HEPES, pH 6.5, 0.25%, 200 mM NaCl) and lysis buffer (10 mM EDTA, 10 mM Tris–HCl, pH 8.1, 0.5% Empigen BB, 1% SDS) supplemented with complete protease inhibitor cocktail (Roche). Chromatin was fragmented to an average size of 200–400 bp using a bio-ruptor (Diagenode) sonicating water bath (10 min ×3, 15 s on and 15 s off). Debris was removed by centrifugation for 10 min at 15 000g. For chromatin immunoprecipitations, chromatin was supplemented with IP buffer (to final concentration of 0.1% sodium deoxycholate, 0.1% Triton X-100, 10 mM Tris–HCl, pH 8.0, 1 mM EDTA) and incubated overnight at 4°C with antibodies or matched IgG controls, pre-bound to protein-G-coated Dynabeads (Invitrogen). ChIPs were washed five times with RIPA buffer (50 mM HEPES, pH 8, 1 mM EDTA, pH 8, 1% NP-40, 0.7% sodium deoxycholate, 0.5 M LiCl), once with TE (10 mM Tris–HCl, 1 mM EDTA) and cross-links reversed in elution buffer (10 mM Tris–HCl, pH 8.0, 1 mM EDTA, 1% SDS) at 65°C for 16–18 h. Eluted samples were diluted in an equal volume of TE, treated with DNAse-free RNAse (Ambion) for 30 min at 37°C, Proteinase K (Invitrogen) for 1–2 h at 55°C and purified using QIAQUICK columns (Qiagen) according to the manufacturer’s instructions.

For re-ChIP assays, a similar protocol was used with the following modifications. After the first ChIP, beads were washed twice with wash buffer 1 (2 mM EDTA, 20 mM Tris–HCl, pH 8.1, 0.1% SDS, 150 mM NaCl) and two washes with wash buffer 2 (2 mM EDTA, 20 mM Tris–HCl, pH 8.1, 0.1% SDS, 300 mM NaCl) and eluted twice with 10 mM DTT before second round of ChIP. Antibodies used for ChIP in this study were monoclonal antibodies TP63 (Santa Cruz: 4A4), TP53 (DO1), TFAP2A (3B5), TFAP2C (H77) or polyclonal TFAP2A (Santa Cruz, C-18). The specificity of the 4A4 antibody and suitability for ChIP-seq was confirmed by western blot of keratinocytes depleted of TP63. The 4A4 recognizes three bands correlating with the predominant deltaNTP63alpha isoform and the less abundant deltaNTP63beta and deltaNp63gamma isoforms (10), which are all effectively depleted by the TP63 siRNA targeting the DNA-binding domain and therefore all p63 isoforms (Supplementary Figure S1A).

Quantitative ChIP-PCR

Quantitative PCR was carried out on ChIP-enriched DNA using FastStart SYBR Green Master (Roche) according to the manufacturer’s instructions. Fluorescence was monitored on a DNA Engine® Peltier Thermal Cycler (Bio-Rad) equipped with a Chromo4 Real-Time PCR Detection System (Bio-Rad) with the following cycling conditions: initial denaturation 95°C for 10 min 40 cycles of 95°C, 15 s; 58°C, 15 s; 60°C, 40 s and melting curve analysis also performed. ChIP enrichment was determined for IP and IgG control for each region in triplicate from at least two independent biological replicates and expressed as percentage input using the pfaffl method. Results were further normalized to any enrichment observed in a previously described negative control region of chromosome 5 (10). All PCR products were further confirmed by separation on a 2% agarose gel. For siRNA recruitment ChIP-qPCR, fold enrichment over input was calculated compared to a negative control region (CCND1) and expressed relative to scrambled control. Primers were batch designed by extracting fasta files of sequences within the TP63 peaks using the fetch genomic sequence tool within the galaxy environment (25,26) and using this as input for batchprimer3 (27) with optimal product size 100 bp, Tm 60°C (for sequences see Supplementary Methods).

Sequencing and peak detection

End-repaired sequencing libraries were prepared from two independent TP63 ChIPs and matched input control using NE Next sample preparation kit (New England Biolabs) according to the manufacturer’s instructions. To remove residual adapters, the resulting adapter ligated DNA was purified with Ampere SPRI beads (Beckman Coulter) and subjected to 16 cycles of adapter-mediated PCR using Hercules II Fusion DNA Polymerase (Agilent). The resulting libraries were size fractionated in a 2% agars gel and purified to obtain 2–400 bp fragments. Each library was validated using an Agilent 2100 bioanalyser and sequenced on the Illumina GAII sequencer according to the manufacturer’s instructions. Fastq files were generated with Illumina pipeline software (OLB v1.8and CASAVA 1.7.0 using the default chastity base call thresholds) and subsequently filtered to remove suboptimal reads (>5% N bases and/or >95% same base) and PCR duplicates. High-quality reads were then mapped against the NCBI37/Hg19 (Ensembl 56) genome with bwa (28) allowing for gapped alignment and maximum five alignment locations for each read. Following removal of inferior alignments (reads partially aligned or with more than two indels or more than three mismatches to the reference sequence).

The MACS algorithm (version 1.3.7.1) (29) was used to analyse the resulting alignments and identify TP63-binding regions, distribution of aligned reads of TP63 ChIP libraries with matched input libraries. The MACS algorithm uses a dynamic Poisson distribution and the input control to account for local biases in the genome and sample. TP63-binding peaks were called using the MACS algorithm, using bam files as inputs, MFOLD = 20 and P = 4 × 10−5. The resulting bed files identifying 9174 and 16 282 were overlapped and 7547 high-confidence TP63-binding regions called when present in both biological replicates (Supplementary Table S1). The individual resulting bed and genome-wide histograms of read pile-ups (wig files) from MACS were visualized using the integrated genome viewer (IGV) (30) and are available from GEO under accession number GSE32061.

Data analysis

As an initial analysis of the genomic location of the 7574 TP63 peaks was carried out using the ChIP-seq toolset (31), this describes the location of each peak as compared to the nearest transcription start site in addition to giving information as to whether the region is within an Intron, exon, etc. (Figure 2A).

Figure 2.

Genome-wide comparison of TP63- and TP53-binding sites. (A) Venn diagram showing overlap of TP63-bindings sites identified by ChIP-seq in HFKs compared with published TP53 ChIP-seq following DNA damage [endogenous TP53 in IMR90 fibroblasts treated with 5FU (54) and exogenous TP53 in SAOS2 cells treated with Etoposide or Actinomycin D (53)]. (B) Examples of TP63 only bound sites (IRF6, HRAS and EGFR) and TP63/TP53 bound sites (MDM2, CDKN1A and PTEN/KLLN). (C) Quantitative PCR ChIP validation of TP63 (top panel) and TP53 (bottom panel) binding to the above sites in HFKs ± adrimaycin treatment and western blot showing effect of adriamycin treatment on TP63 levels and TP53 levels and activation (inset panel). (D) Venn diagram comparing genes within 25 kb of TP63 only, TP63/TP53 and TP53 only binding sites. (E) Gene ontology analysis results (David—GOBPFat), comparing TP63 only (1420) and pooled TP63/TP53 (2620) gene lists.

A TP63-binding motif was identified from the 7574 TP63 high confidence binding sites using the de-motif motif identification module of the Partek® genomic suite software (Partek® software. Copyright, Partek Inc., St Louis, MO) using default parameters (Figure 2B). The advantage of using this software was that it readily accommodated the inclusion of all of the 7574 TP63-binding regions. Motifs were generated from the resulting position-weighted matrices using the SeqLogo package (32) from Bioconductor within the R-environment (Figures 2B and Supplementary Figure S2). The distance of TP63-binding motifs from the centre of peaks was determined and plotted as a frequency distribution with bin size 10 bp (Figure 2C).

A conservation plot for the 3 kb regions surrounding the 7574 TP63 peaks was generated by plotting PhastCons for vertebrates using the conservation plot tool from CEAS (33) within the Cistrome galaxy environment (25,26).

To correlate TP63 peaks with potential target genes, we chose to annotate each peak to any gene within 25 kb of the identified TP63 peak. This was achieved by overlapping peaks with Refseq genes extended by 25 kb downloaded from the UCSC browser (34) within the Galaxy environment (25,26).

Seqminer (Version 1.3) (35) was used to extract histone modification data and for K-means clustering (raw) of this data for respective sets of binding sites.

To identify potential TP63 co-factors, the 7574 TP63-binding regions were input to the Cis-Element Annotation System web-server CEAS (33). This counts the occurrences of ∼800 characterized motifs from the Transfac and Jaspar database within ChIP-regions and reports those, which are significantly enriched as compared to their occurrences in the genome. To identify the number of peak containing predicated binding regions, we mapped the predicted motifs to the ChIP. The distance of predicted transcription factor-binding motifs from the centre of the ChIP-regions was calculated and plotted as a frequency distribution with bin size 10 bp (Figure 5C).

Figure 5.

TFAP2A is functionally associated with TP63-binding sites. (A) Effect on the IRF6 enhancer activity of the addition of exogenous TP63 in H1299 cells as measured by luciferase activity. All luciferase graphs represent the mean fold increase of at least three independent biological ± SEM; Student’s t-test *P < 0.05. (B) Luciferase assay measuring the effect of the addition of increasing amounts of TFAP2A on TP63-mediated activation of the IRF6 enhancer region in H1299 cells. (C) Western blot analysis of exogenous TP63 and TFAP2A expression in H1299 cells. (D) Luciferase assay measuring the effect of the addition of increasing amounts of TFAP2A on TP63-mediated activation of the PVRL1 enhancer region in H1299 cells. (E) Quantitative ChIP-PCR comparing TP63 and TFAP2A binding to IRF6, FGFR2, JAG2, PVRL1 and PDGFC associated binding sites, expressed as fold enrichment compared to a negative control region within the same ChIP. Data represent the mean of two independent biological replicates. (F) Western blot validation of TP63 and TFAP2A depletion in HFKs used in recruitment ChIPs. (G) Quantitative ChIP-PCR comparing TP63 and TFAP2A binding to IRF6, FGFR2, JAG2 and PVRL1 associated binding sites in HFKs depleted for TP63 or TFAP2A compared to scrambled control. Results calculated as fold enrichment compared to a negative control region within the same ChIP and normalized to scrambled control to compare between antibodies and experiments. Data represent the mean of two independent biological replicates. (H) Quantitative ChIP-PCR of second negative control binding site (BCL2 upstream) in HFKs depleted for TP63 or TFAP2A expressed as fold change compared to scrambled control. Data represent the mean of two independent biological replicates. (I) Quantitative RT-PCR quantification of expression of CL/P associated genes in TFAP2A-depleted human foreskin keratinocytes compared to scrambled control. Graph represents the mean fold increase of at least three independent biological ± SEM; Student’s t-test *P < 0.05, **P < 0.01.

Microarray analyses

Raw data were downloaded from the Gene Expression Omnibus (GEO) (36) for MCF10A cells retrovirally depleted for all TP63 isoforms and scrambled control (GSE20286) (37) and murine keratinocytes siRNA depleted for all TP63 isoforms and scrambled control (GSE10564) (38). Raw.cel files were imported into Partek software using gcRMA normalization. Changes in the expression level were determined by comparing scrambled control and TP63-depleted arrays using one-way ANOVA and a cut-off of P < 0.05. This was further shortlisted to probesets exhibiting a >1.5fold change of Log2 transformed data.

Heatmaps were generated from extracted data using Matrix2png web application (39).

Hypergeometric distribution for gene list overlaps (40), which predicts the probability of the overlap of the two gene lists based on chance alone, given the length of the two lists and the number of genes that could have been present on both lists (R-code available on request).

Quantitative RT-PCR

RNA extraction was carried out using Trizol (Invitrogen) according to the manufacturer’s instructions. RNA was reverse transcribed to cDNA using the Transcriptor high-fidelity cDNA synthesis kit (Roche), with random hexamers according to the manufacturer’s instructions. Amplification of PCR products was monitored using Lightcycler® 480 SYBR Green I Master (Roche) according to the manufacturer’s instruction and fluorescence monitored on a Roche 480 Lightcycler and melting curve analysis also performed. In brief, cDNA samples were diluted 1:50 and quantified compared to a standard dilution series using the absolute relative quantitation method. The cycling conditions were as follows: initial denaturation 95°C for 10 min 45 cycles of 95°C, 15 s; 58°C, 15 s; 72°C, 20 s. Expression levels were assessed in triplicate, normalized to RPLPO levels and graphs represent the mean of three independent biological replicates, ± standard error of the mean (SEM) (for primer sequences see Supplementary Methods).

Luciferase assays and cloning

The regions contained within theTP63 peaks associated with IRF and PVRL1 (PVRL1 = 844 bp fragment from Intron-1, IRF6 = 880 bp fragment located 10 kb upstream of IRF6) were amplified and cloned upstream of the luciferase reporter in the pGL3 promoter vector (Promega). Plasmids were sequenced to confirm inserts (for primer sequences see Supplementary Materials). TP63 plasmids used were as previously described (41). The TFAP2A construct was obtained from Addgene [plasmid SP(RSV); 12100) (42)], and matched TFAP2A and TFAP2C image clones in pCMVsport6 were obtained from MRC Geneservice (Cambridge). For luciferase assays, H1299 cells were plated at 7 × 104 cells per well of a six-well plate and transfected the following day with a total of 1.6 µg DNA. Each well was transfected with 500 ng of luciferase reporter construct, 50 ng renilla luciferase, 300 ng of TP63 construct (pcDNA3.1 CMV) and 300/750 ng TFAP2A construct and made up to 1.6 µg with empty vector. Transfections were carried out with Genejuice (Novagen) in H1299 cells according to the manufacturer’s instructions. For HFK transfection, cells were plated at 7.5 × 104 per well in 12-well plate and transfected with 0.8 µg DNA using Fugene6 (Roche), according to the manufacturer’s instructions. Transfected cells were harvested 48 h post-transfection and assayed for firefly and renilla luciferase activity using the Dual Luciferase Reporter Assay (Promega). Each assay was carried out in triplicate, and luciferase activity is measured as the ratio of firefly/renilla luciferase activity and normalized to control. Results represent the mean of three independent biological replicates ± SEM.

Western blot analysis

Western blots were carried out as previously described (43). Primary antibodies used were Santa Cruz monoclonal anti-TP63 (4A4), TP53 (DO1), CDKN1A (C19) TFAP2A (3B5), TFAP2C (6E4/4) and polyclonal anti-TP63 (H129), TFAP2A (C-18), TFAP2C (H77), polyclonal TP53-pS15 (Cell Signalling), anti-monoclonal anti-beta-actin (Sigma). Secondary antibodies used in this study were goat anti-mouse and rabbit-HRP (Santa Cruz). Luminescence was revealed by incubation with Western Lightning ECL (Perkin-Elmer) and signal detected on an Alpha Innotech FluorChem™ SP imaging system.

Organotypic raft section staining

Tissue embedding, haematoxylin and eosin staining, and indirect immunofluorescence were carried as previously described (44). Primary antibodies used were polyclonal anti-keratin 1, Loricrin (Covance), monoclonal anti-Keratin 10 (Biodesign), transglutaminase-1 (Santa Cruz), filaggrin (Neomarkers), BrDU (Santa Cruz). The number of proliferating cells in organotypic raft cultures was determined by counting the BrDU positive cells for a minimum of 10 fields of view (1000 µM).

RESULTS

Genome-wide identification of TP63-binding sites in primary human foreskin keratinocytes

Chromatin immunoprecipitation (ChIP) was carried out from two independent pooled batches of actively proliferating HFKs with a ChIP-validated TP63-specific monoclonal antibody (4A4) (Supplementary Figure S1A), which detects all TP63 isoforms. Prior to library preparation and sequencing, successful target enrichment, was confirmed by semi-quantitative (Supplementary Figure S1B) and quantitative PCR (ChIP-qPCR) (Supplementary Figure S1C) of two well-characterized constitutive TP63-binding sites, one 2.2 kb upstream of the CDKN1A transcription start site (TSS) and a second in intron 1 of MDM2 (9,45). Subsequent to sequencing, visual inspection of these control regions of the CDKN1A and MDM2 gene using the Integrated Genome Viewer (IGV) (30) indicated significant enrichment of reads in these positive control regions in both replicates (Figure 1A).

Figure 1.

Genome-wide characterization and validation TP63-binding sites. (A) ChIP-seq results for the CDKN1A and MDM2 loci for two independent biological replicates. (B) Comparison of tag number between peaks identified by MACS common to the two independent TP63 ChIP-seq experiments. (C) Quantitative PCR validation of range of 10 TP63-binding sites from a range of peaks (Intra/intergenic, proximal promoter) and positive and negative control (data represents the mean of three biological replicates expressed as % input ± SEM). (D) Correlation of qPCR validation and average tag number from each peak from ChIP-seq for 42 validated target regions. (E) Biological Processes GO result from Database for Annotation, Visualization and Integrated Discovery (DAVID).

To identify genome-wide TP63-binding sites, we utilized the model-based analysis of ChIP-seq (MACS) peak detection algorithm (29) for the two independent TP63 ChIPs and their matched inputs as control. To generate a list of high confidence peaks, we compared the peaks called by MACS from the two biological replicates and identified 7574 overlapping regions present in both replicates (Supplementary Table S1). Significant correlation of peak heights was observed between the two experiments (r = 0.7724, P < 0.0001) (Figure 1B).

We confirmed TP63 binding to 10 regions from a range of peak heights and locations (inter/intragenic, proximal promoter) by quantitative PCR (ChIP-qPCR) of three independent TP63 ChIPs from cycling HFKs, all of which, exhibited significant enrichment compared to an IgG control and a negative control region (Figure 1C). To date, we have validated 43 binding sites, of which only one did not show significant enrichment compared to control (data not shown). Comparison of peak height and ChIP-qPCR validation of these 42 binding sites reveals positive correlation between average peak height and qPCR validation (r = 0.6364, P < 0.0001) (Figure 1D).

The majority of TP63 peaks are located within genes or in regions upstream of core promoters (Supplementary Figure S2A and B). De novo motif analysis of all 7574 TP63-binding regions reveals a TP63-binding motif, similar to those previously described (Supplementary Figure S2C and Table S2) (46). The TP63-binding motifs are focused in the middle of the predicted peaks (Supplementary Figure S2D) and are highly conserved across mammals (Supplementary Figure S1E). Of note, the TP63-binding motif is more degenerate than that described for p53 (47). Of the 7574 TP63 peaks, 6424 (85%) contain at least one occurrence of the full TP63-binding motif, whereas 1150 (15%) peaks do not contain the identified TP63 consensus-binding motif (Supplementary Table S2).

Comparison of the TP63-binding sites with the Broad histone modification datasets available for normal human epidermal keratinocytes (NHEKs) from the ENCODE project (48) reveals that TP63-binding sites are globally enriched for H3K4 mono and di-methylation and low tri-methylation in both flanking regions, which is characteristic of enhancer regions (Supplementary Figure S3A, S3B, S3C, Clusters i–iii, 6921 total sites) (49) [K-Means clustering using SeqMiner (35) (Supplementary Figure S4C)]. However, a small proportion are associated with histone profiles characteristic of promoter regions indicated by high levels of H3K27 acetylation and H3K4 trimethylation (49) (Clusters v–vi, 653 total sites) (Supplementary Figure S3C), the binding profiles of which can be attributed to gene preferentially transcribed from the positive strand (Cluster vi), negative strand (Cluster v) or in both directions (Cluster iv) (Supplementary Figure S3C).

Annotation of TP63-binding sites

Histone modification data and preliminary inspection of TP63-binding regions indicated that multiple binding sites lie within potential bidirectional promoters, regions of gene overlap or close proximity to multiple genes (data not shown). Mapping each peak to a single gene may consequently lead to omission of regulatory interactions. In addition, many of the TP63 binding sites were in regions distal to the nearest gene (Supplementary Figure S2B); therefore, we annotated peaks to any Refseq gene within 25 kb of the identified peaks (Supplementary Table S3). Using this annotation strategy, we identified 6172 potential target genes within 25 kb of a TP63-binding site. This likely omits longer range interactions; however, for the purpose of this study represents a reasonable compromise.

We utilized the Database for Annotation, Visualization and Integrated Discovery website (DAVID) (50,51) to carry out gene ontology (GO) analysis to gain a broad overview of the gene pathways potentially regulated by TP63 and to support our annotation strategy. Iterative analysis of the 6172 associated gene list (within 25 kb) or alternative gene lists within 10 kb (4420 genes) and 5 kb (3.721 genes) (Supplementary Table S3) of TP63-binding sites indicated enrichment of genes involved in biological processes including adhesion, epidermis/ectoderm development, cytoskeletal organization, apoptosis, intracellular signalling, cell motion and proliferation (Figure 1E). In support of our annotation of genes within 25 kb, a pathway known to be associated with p63, epidermal development, is ranked third in GO analysis of biological processes and includes 95 genes, whereas it is ranked 35th (63 genes) and 42nd (38 genes) in the 10 and 5 kb analyses, respectively (Supplementary Table S4), indicating that TP63 likely plays an important role in long-range interactions. In addition, analysis of TP63 binding sites using the Genomic Regions Enrichment of Annotations Tool (GREAT) (52) with iterative decreasing maximum distances for annotation of peaks (1 Mb, 500 kb, 250 kb, 100 kb and 25 kb) indicates increasing enrichment for epidermal-related biological processes up to 500 kb (Supplementary Table S5). Inspection of these gene lists indicates that TP63 potentially regulates genes with key roles in all aspects of epidermal development, growth control and differentiation and truly represents a master regulator of stratifying epithelia.

TP63 and TP53 bind to distinct and overlapping genomic regions

The canonical role of the TP63 in maintaining proliferation by repressing TP53-activated target genes in cycling cells, and the switch between TP63 and TP53 binding at these sites upon DNA damage has been studied for a small number of anti-proliferative and pro-apoptotic genes (9,45). To estimate the proportion of TP63-binding sites identified by ChIP-seq that are potentially co-ordinately regulated by TP53 in response to DNA damage, we compared our HFK TP63-binding sites with those identified in two recent TP53 ChIP-seq studies [pooled binding sites from etoposide/actinomycin D-treated SAOS2 osteosarcoma cells expressing exogenous wild-type TP53 (53) and 5FU-treated IMR90 normal lung fibroblasts (54)]. We used DNA damage response as it is a well-recognized function of TP53. Employing this strategy, we identify 1316 of our 7574 TP63-binding sites that are bound by TP53 in response to DNA damage and 6258 TP63 sites not bound by TP53 (Figure 2A) (Supplementary Table S1). To confirm the validity of this approach we carried out quantitative ChIP-PCR for TP63 and TP53 on three TP63 only sites (IRF6, HRAS, EGFR) and three TP63/TP53 sites (MDM2, CDKN1A, PTEN) (Figure 2B) from HFKs in the presence and absence of Adriamycin (Figure 2C). TP63 is significantly enriched on all six sites and is decreased upon Adriamycin treatment when TP63 is degraded (Figure 2C). As expected TP53 is bound to the three predicted TP53/TP63 sites and increased upon Adriamycin treatment, whereupon TP53 levels and activation are increased as measured by p53 protein levels and serine-15 phosphorylation, respectively (Figure 2C). In contrast, no enrichment is observed for TP63 only sites for TP53 (Figure 2C).

To get an overview of the genes with TP63 and/or TP53 binding sites, we utilized Gene Ontology analysis in DAVID (50,51). We compared the genes within 25 kb of these sites and derived three gene lists; Genes associated with TP63 only (4229), TP63/TP53 common genes (1946) and TP53 only genes (1348) (Figure 2D). The validity of the approach is supported in that epidermis development is the highest ranked term for the TP63 only genes, whereas within the TP63/TP53 common genes, this is ranked 260th (Figure 2E and Supplementary Table S4). Furthermore, terms such as negative regulation of cell proliferation and regulation of programmed cell death are highly ranked in the TP63/TP53 list, indicating that TP63 binds constitutively to a large proportion of DNA damage induced TP53-binding sites in normal cycling HFKs. Furthermore, in support of the GO analysis of KEGG pathways, analysis with DAVID (50,51) on the same gene lists indicates that the ‘TP53 signalling pathways’ was the top-ranked KEGG pathway for the TP63/TP53 common gene list, but this pathway was not enriched when the same analysis was carried out on either the TP63 or TP53 only lists (Supplementary Table S4). In summary, a subset of TP63-binding sites are predicted to be co-ordinately regulated by TP53 in response to DNA damage, while the majority of the TP63-binding sites identified in cycling HFKs are not bound by TP53 in response to DNA damage.

Integrative analysis identifies TP63 target genes associated with cleft palate

Analysis of >6000 potential TP63 targets is prohibitive, and the binding of TP63 to these sites may only affect transcription in specific conditions or physiological processes such as DNA damage or differentiation. For the purpose of this study, we adopted an integrative strategy utilizing publicly available microarray data for epithelial cells depleted for TP63 for which the raw data were available from GEO (36). We obtained raw data files for murine keratinocytes and MCF10A basal breast epithelial cells depleted for all isoforms of TP63 and scrambled control (37,38). These datasets were chosen since TP63 has been shown to be important for the differentiation of these cells (2,3,10,16,37,41) and like HFKs they express high levels of TP63 and have wild-type p53. Both datasets were analysed as described in ‘Materials and Methods’ section and differentially expressed genes determined using one-way ANOVA (Supplementary Table S6).

To identify direct TP63 target genes, we overlapped the list of genes which exhibited a >1.5-fold change compared to scrambled control (Log2 transformed) and a P < 0.05 with the list of 6171 genes within 25 kb of a p63-binding site (Figure 3A) and compared the significance of overlaps based on the sampling of a hypergeometrically distributed random variable (Supplementary Figure S4A). Examining the murine keratinocyte dataset revealed 1484 genes were altered >1.5-fold, of which 611 were within 25 kb of TP63-binding sites (hypergeometric test P = 0) (Supplementary Figure S4B). In the MCF10A dataset, there were 2053 genes significantly altered, 776 of which were within 25 kb of TP63-binding sites (hypergeometric test P = 1.13 × 10−11) (Supplementary Figure S4C). Combining these list identifies 1213 TP63-bound genes significantly changed in either murine keratinocytes or MCF10A experiments, 174 of which are altered in both cell types (Figure 3A) (Supplementary Table S7) (Hypergeometric test P = 0). The differences in the genes regulated is likely reflective of differences in the cell types and may be reflective of the TP63 isoforms and co-factors expressed. We have previously shown that the MCF10A cell line express a similar complement of TP63 isoforms to that expressed in HFKs (41).

Figure 3.

Integrative analysis of ChIP-seq and microarray data identifies multiple TP63 targets involved in palatal fusion. (A) Venn diagram of overlap of TP63 ChIP-seq genes with genes whose expression is significantly altered in TP63-depleted MCF10A and murine keratinocyte microarray experiments. (B) Genetic Association Database enrichment results from DAVID. (C) GeneGO pathway analysis using MetacoreTM software. (D) Venn diagram of overlap of TP63 ChIP-seq genes, microarray overlap gene list and cleft palate (CL/P) manually curated gene list. (E) Heatmap of microarray data for a 41-cleft palate gene list from TP63-depleted MCF10A and murine keratinocyte microarray experiments. (F) Quantitative PCR ChIP comparison of TP63 binding sites with Pan TP63 targeting 4A4 monoclonal and TP63-alpha targeting polyclonal antibody (H129). Data represent the mean of three biological replicates expressed as % input ± SEM. (G) Quantitative RT-PCR for CL/P-associated genes in TP63-depleted human foreskin keratinocytes compared to scrambled control (Graphs represents the mean fold increase of at least three independent biological ± SEM; Student’s t-test *P < 0.05, **P < 0.01).

We concentrated on genes implicated in palatal fusion as this is one of the phenotypes observed in patients with TP63 mutations and our GO analysis of genes associated with TP63 binding (6172 genes) identified enrichment for genes whose mutation is associated with cleft palate in the genetic association database (GAD) (55) (Figure 3B). Similar to epidermal regeneration, palatal fusion requires intricate control of multiple biological processes including growth, adhesion, apoptosis, motility and epithelial/mesenchymal interactions, all of which were identified in the GO analyses. Regulation of palatal fusion has additionally been shown to require TGF-beta signalling (56), which we have recently shown to be regulated by TP63 in the MCF10A breast cell line (41) and the TGF-beta/Wnt signalling pathway was the most enriched GeneGO pathway analysis using the MetacoreTM software (Figure 3C) (Supplementary Table S4).

Since gene ontology is limited by the annotations within the databases interrogated, we chose to generate a more comprehensive manually curated list of cleft palate associated genes. We identified Refseq genes whose mutation was associated with cleft palate in the Online Mendelian Inheritance in Man database (OMIM) (57) and the Genetic Association Database (GAD) (55). We also incorporated genes whose perturbation had been shown to result in cleft palate in the Mouse Genome Informatics (MGI) phenotype database (58) and those described in a review of the molecular biology of palate development (59). This resulted in a list of 304 genes associated with cleft palate (Supplementary Table S8), which we compared to our results from our p63 ChIP-seq and the microarray analysis and identified 41 genes that were significantly altered in at least one of the two array experiments and associated with TP63-binding sites in our TP63 ChIP-seq (Figure 3D). The changes in expression of these 41 genes are graphically represented in a heatmap of fold change in the two TP63 depleted cell types relative to scrambled control (Figure 3E) (for genes that did not reach significance in one of the datasets, the probeset with the lowest P-value is included). Of these 41 genes, 25 are similarly up/down-regulated in response to TP63 depletion, 14 are differentially regulated (oppositely up/down upon TP63 depletion in two cell types, likely reflective of tissue-specific regulation) and 2 are only present on one of the two arrays (PDSS2, GSDMC). These include genes associated with both syndromic and non-syndromic CL/P in humans such as fibroblast growth factor receptor 2 (FGFR2) (Apert OMIM 101200; Crouzon OMIM 123500 and Pfeiffer OMIM 101600 syndromes) (60,61), Poliovirus receptor-like 1 (PVRL1) (CLPED1; OMIM 225060) (62,63) and Interferon regulatory factor 6 (IRF6), the most commonly mutated gene in humans with cleft palate (Van Der Woude VWS: OMIM 119300 and Popliteal Pterygium syndromes PPS: OMIM 119500) (64). Interestingly, two groups independently identified IRF6 as a TP63 target during preparation of this manuscript (65,66).

We validated TP63 binding to a subset of the 41 cleft palate associated genes by ChIP-qPCR and if more than one peak occurred within 25 kb of a gene, the highest ranked peak was chosen for validation. Binding was assessed by ChIP-qPCR of three independent TP63 ChIPs from cycling HFKs using the pan-TP63 4A4 antibody and a second polyclonal antibody (H129), which recognizes the alpha isoforms. Significant enrichment compared to negative control regions was observed for eight regions examined with both antibodies on CL/P gene predicted binding regions and positive controls (MDM2, CDKN1A) (Figure 3F). To determine if TP63 affected constitutive expression of these targets, we examined mRNA levels of ten of these genes in HFKs depleted for TP63 with an siRNA targeting the TP63 DNA-binding domain and therefore all TP63 isoforms (10). Quantitative RT-PCR indicated that expression of the 10 genes was altered (eight significantly) (Figure 3G), suggesting that TP63 binding contributes to transcription of these genes in HFKs.

Identification of TFAP2A as a co-factor for TP63

We utilized transcription factor binding site enrichment analysis utilizing the Cis-Element Annotation System (CEAS) (33), to identify potential TP63 co-factors. This revealed that the most significantly enriched transcription factor-binding motif was AP2alpha (Figure 4A), along with a number of other AP-2 motifs represented in the top 100 significantly enriched motifs. Examination of predicted TFAP2A sites within the TP63-binding sites identified using CEAS (33) (Figure 4A) indicates that 37.8% (2867) of the 7574 TP63-binding sites contain at least one of the pooled predicted AP-2 sites (Supplementary Table S9). These predicted AP-2-binding motifs are enriched towards the centre of the TP63-binding sites (Figure 4B) but are not coincident with TP63 motifs as there is an average of 110 bp between the TP63 and AP-2 motifs (Supplementary Figure S5A). Independent analysis of TP63-peak regions using the CisFinder (tool for finding over represented transcription factor binding motifs) (67) returned similar results for predicted TFAP2A (Supplementary Figure S5B). Therefore, given the overlapping phenotypes observed in patients with TP63 and AP2alpha mutations and observations from in silico analyses, we hypothesized that TFAP2A may represent a potential co-factor for TP63.

Figure 4.

Identification of TFAP2A as a potential TP63 co-factor. (A) Table of results from transcription factor motif enrichment analysis of 7574 TP63-binding sites. (B) Cumulative frequency distribution plots comparing distance of predicted TFAP2A sites and TP63-binding motifs from the centre of TP63-binding sites/peaks. (C) Comparison of localization of TP63-binding sites with or without predicted AP-2 sites. (D and E) Quantitative PCR ChIP validation of TP63 (D) and TFAP2A (E) binding to TP63 sites associated with CL/P genes (data represent the mean of three biological replicates expressed as % input ± SEM). (F) Semi-quantitative PCR results for sequential re-ChIP assay for TP63 and TFAP2A for the IRF6, FGFR2, TGFB1 and PVRL1 associated TP63-binding sites, showing both factors co-ChIP.

Comparing the subset of TP63-binding regions associated with predicted AP-2 sites indicates that these regions are frequently closer to TSSs than those devoid of predicted AP-2 sites (Figure 4C, 0–1 kb) and this correlates with higher H3K4 trimethylation, characteristic of promoter regions (Supplementary Figure S5C) (49). Furthermore, comparison of TP63/AP-2 associated genes shows significant enrichment within the 1213 TP63 bound genes found to be altered in at least one of the TP63-depleted microarray studies (Figure 3A) (hypergeometric test P = 3.17 × 10−05) indicating that TP63-controlled genes are enriched for AP-2 sites. In addition, 393 of the 1354 TP63 sites, predicted to be bound by TP53 in response to DNA damage, contain predicted AP-2 sites. Analysis of the regulated genes indicates a similar but less significant enrichment is observed for the subset of these genes associated with predicted TP53-binding sites (hypergeometric test P = 0.00023).

To determine if TFAP2A is associated with TP63-binding regions in human keratinocytes, we carried out ChIP for TP63 and TFAP2A from cycling HFKs using a polyclonal TFAP2A antibody (C-18), followed by ChIP-qPCR. We examined TP63 binding to the TP63 binding regions of 11 of the 41 cleft palate-associated genes which contained predicted AP-2 sites (IRF6, FGFR2, PVRL1, JAG2, TGFA, MMP14, EGFR, CD44, TGFB1, PTPRF and IL1A), two cleft palate-associated genes without predicted AP-2 sites (PDGFC and ETS1) and a negative control region on chromosome 5 (Figure 4D and E). Significant enrichment was observed in TFAP2A ChIPs for TP63-bound regions with predicted AP-2 sites compared to IgG and negative control regions. No such enrichment was observed for TP63bound regions without predicted TFAP2A sites (Figure 4D and E). To examine if TP63 and TFAP2A are simultaneously associated on TP63-binding regions genes, we carried out sequential ChIP for TP63 followed by TP63 or TFAP2A and show enrichment of both TP63 and TFAP2A compared to IgG controls for the IRF6-, PVRL1-, FGFR2-binding regions (Figure 4F), indicating that TFAP2A and TP63 can bind to these regions simultaneously, suggesting that they may cooperatively regulate target genes.

TFAP2A co-operates with TP63

To test this hypothesis, we generated a construct containing the IRF6 enhancer spanning the TP63-binding region (encompassing multiple predicted TFAP2A sites and SNP) (Figure 2B), cloned upstream of the luciferase gene under the control of the SV40 promoter (PGL3 promoter). Initially, H1299 cells, which are p53 null and express low levels of TP63, were co-transfected with this luciferase construct and one of the six TP63 isoforms along with renilla luciferase as a transfection control. Co-transfection of all six TP63 isoforms resulted in increased transactivation of the luciferase reporters to varying degrees, with the predominant deltaNTP63alpha isoform eliciting the smallest increase and the deltaNTP63beta isoform the greatest (Figure 5A). We then sought to determine the effect of the addition of increasing amounts of TFAP2A on deltaNTP63beta-mediated activation of this enhancer region (Figure 5B). DeltaNTP63alpha was also included since it is the predominant isoform expressed in keratinocytes and stratifying epithelia. Our results indicate that the addition of increasing amounts of TFAP2A enhances deltaNTP63beta-mediated activation of the IRF6 enhancer in H1299 cells but has no such effect on deltaNTP63alpha-mediated activation (Figure 5B). We also carried out western blot analysis of exogenous proteins and show that the enhancement of transactivation of the IRF6 enhancer is not caused by TFAP2A-mediated changes in TP63 transgene expression (Figure 5C). To determine if this effect occurs in other TP63/TFAP2A-bound regulatory regions, we cloned the TP63-binding region from Intron 1 of the PVRL1 gene, which is predicted to contain multiple TFAP2A-binding sites and obtained similar results (Figure 5D).

To determine if the recruitment of TFAP2A requires TP63 or vice versa, we carried out ChIP-qPCR from cycling keratinocytes depleted for TP63 or TFAP2A using siRNA. For these studies, we used an alternative TFAP2A monoclonal antibody (Santa Cruz 3B5), with better specificity for TFAP2A (Supplementary Figure S6). In these assays, we observed significant fold-enrichment over the negative control region for the four TP63/TFAP2A-bound regions examined (IRF6, FGFR2, JAG2 and PVRL1) with both the TP63 (4A4) and TFAP2A (3B5) antibodies (Figure 5E), but we do not see any enrichment by TFAP2A in a TP63-binding site devoid of predicted TFAP2A-binding sites (PDGFC). As expected depletion of TP63 or TFAP2A (Figure 5F) results in their decreased relative binding to the IRF6, PVRL1, FGFR2 and JAG associated regions, indicating that these interactions are specific (Figure 5G). TP63 depletion results in an almost equivalent decrease in TFAP2A recruitment to all four of these regions (Figure 5G). In contrast, recruitment of TP63 to the JAG2-associated region is decreased by TFAP2A depletion, whereas recruitment of TP63 to the IRF6, PVRL1 and FGFR2 is unaffected. To check for normalization effects, we examined levels of a second negative control region (2 kb upstream of BCL2) and found no difference (Figure 5H). Finally, depletion of TFAP2A resulted in a reduction in the expression of IRF6, FGFR2B and PVRL1 in HFKs (Figure 5I). In summary, TP63 and TFAP2A factors co-operate to increase expression of certain TP63-regulated genes, and the recruitment of TFAP2A is affected by the presence of TP63 on the promoter/enhancer regions containing both binding sites. Interestingly, it is the beta and gamma isoforms of TP63, which exhibit the greatest transcriptional activation of the IRF6 enhancer and cooperate with TFAP2A. This suggests that these isoforms, despite their relatively low abundance compared with the predominant deltaNTP63alpha isoform may play an important role in hetero-complexes the activity of which may be regulated by the transcriptional inhibitory domain of deltaNTP63alpha.

TFAP2A and TFAP2C can bind to overlapping but distinct subsets of TP63-binding regions

We had focused our studies on the TFAP2A AP-2 family member based on the biology described earlier; however, other AP-2 family members have the potential to bind to the same sites (21). TFAP2C has been shown to be expressed throughout the murine epidermis in contrast to TFAP2A, which is predominantly expressed in basal layers (68). We found that TFAP2C was indeed expressed in cycling HFKs, but this was at ∼30-fold lower levels than TFAP2A as assessed by qPCR estimation of copy number (Figure 6A). These low levels of expression were further confirmed by western blot analysis of cycling HFKs with TFAP2C-specific antibodies, which are validated for western blot (6E4/4) and ChIP (H77) (69) (Supplementary Figure S6).

Figure 6.

TFAP2C interacts with a subset of TP63-binding sites. (A) Quantitative PCR estimation of TFAP2A and TFAP2C mRNA copy number in cycling HFKs. Data represent the mean of three independent biological replicates ± SEM. (B) Quantitative PCR ChIP validation comparing TFAP2A (3B5) and TFAP2C (H77) interaction with a subset of TP63-binding regions associated with CL/P genes. Data represent the mean of three biological replicates expressed as % input ± SEM). (C and D) Luciferase assay measuring the effect of the addition of TFAP2A or TFAP2C to six TP63 isoforms on activation of the IRF6 (C) or PVRL1 (D) enhancer region in H1299 cells. (E) Luciferase assay measuring the effect of the addition of TFAP2A or TFAP2C on TP63-mediated activation of the IRF6 enhancer region in primary human foreskin keratinocytes (HFKs). Luciferase graphs represent the mean fold increase of at least three independent biological ± SEM; Student’s t-test *P < 0.05, **P < 0.01.

Despite the relatively low levels of TFAP2C compared to TFAP2A, this does not preclude it from playing an important role in transcription of the TP63 genes associated with TFAP2A binding. We therefore compared TFAP2A (3B5 antibody) and TFAP2C (H77) binding to a subset of these regions by ChIP-qPCR compared to a negative control region in cycling HFKs. This analysis reveals that TFAP2A and TFAP2C are both enriched on the IRF6- and PVRL1-associated binding sites in cycling HFKs; however, TFAP2C is not enriched on the FGFR2-, JAG2- and CDKN1A-associated regions, which are bound by TFAP2A (Figure 6B). This correlates with the ability of TFAP2A and TFAP2C to cooperate with TP63 to activate the IRF6 and PVRL1 luciferase enhancer constructs in H1299 cells with TFAP2C exhibiting the greater enhancement (Supplementary Figure S7A and D), and this was not due to changes in TP63 transgene expression (Supplementary Figure S7B and C). Furthermore, TFAP2A and TFAP2C factors exert differential effects on TP63-mediated transactivation of the IRF-6 and PVRL-1 enhancers in H1299 cells, dependent on the specific TP63 isoform expressed (Figure 6C and D). In general, the deltaNTP63 isoforms co-operate with the two AP-2 factors, whereas the TAP63 isoforms do not. This is consistent with the fact that in primary HFKs deltaNTP63 isoforms are expressed, whereas the TAP63 isoforms are at barely detectable levels (70). To examine the effects of TFAP2A and TFAP2C addition in the context of the endogenous complement of TP63 isoforms, we carried out similar luciferase assays as described above in HFKs. Interestingly, the addition of TFAP2A or TFAP2C alone consistently but not significantly increased activation of the IRF6 enhancer (Figure 6E) and both TFAP2A and TFAP2C significantly augmented deltaNTP63beta-mediated activation of the IRF6 enhancer (Figure 6E). Taken together, these results indicate that TFAP2A and TFAP2C can bind to different but overlapping subsets of TP63-binding regions and differentially modulate transcriptional activities of TP63 isoforms.

Depletion of TFAP2A and TFAP2C inhibits differentiation of organotypic raft cultures

Finally, we sought to determine the effects of acute depletion of TFAP2A and TFAP2C on growth and differentiation in organotypic raft cultures. Organotypic raft cultures of HFKs transiently depleted with siRNAs specific for TFAP2A, or TFAP2C (69) (20) or TP63 (10) were generated (Figure 7A). TP63 depletion as described previously (10) resulted in a complete absence of early, intermediate and late markers of differentiation (KRT1, TGM1, FLG, respectively) (Figure 7B). Depletion of TFAP2A or TFAP2A also resulted in inhibition of differentiation as measured by morphology and expression of the same markers of differentiation. However, the expression of markers of differentiation such as TGM1 and KRT1 is diminished to a lesser extent in TFAP2C compared to TFAP2A-depleted cells (Figure 7B). Furthermore, there was no significant effect on proliferation in either TFAP2A or TFAP2C depleted rafts as assessed by number of BrDU-incorporating cells (Figure 7C), in contrast to TP63 depletion, which results in hypoproliferation as previously described (10).

Figure 7.

TFAP2A and TFAP2C are required for efficient differentiation of organotypic raft cultures. (A) Western blot analysis of TFAP2A, TFAP2C, TP63 and actin loading control of the HFKs transfected with TFAP2A, TFAP2C, TP63 targeting and scrambled control siRNA. (B) Sections of organotypic raft cultures generated from TFAP2A-, TFAP2C- or TP63-depleted HFKs stained for haematoxylin & Eosin (H&E) and indirect immunofluorescent staining of early [keratins 1 (KRT1), intermediate (transglutaminase-1 (TGM1)] and late [Filaggrin (FLG)] markers of differentiation (scale bar = 100 µM). (C) Number of bromodeoxyuridine (BrDU) incorporating cells in organotypic raft culture assessed by immunofluorescent staining. Graph represents the average number of BrDU-positive cells in the basal epithelial layer from organotypic raft cultures per 1000 µM, expressed as percentage scrambled control (mean ± SE at least 10 counts each of two independent biological replicates).

DISCUSSION

TP63 is a master transcriptional regulator of genes involved in development, growth and differentiation of stratifying epithelia, although the precise pathways controlled by this factor are unclear. Here, we report genome-wide mapping and identification of 7574 TP63-binding sites in human neonatal foreskin keratinocytes. We use integrative in silico and in vitro approaches to identify direct TP63 target genes, regulated processes and co-factors.

The genomic distribution of our 7574 TP63 peaks is in agreement with results from the burgeoning ChIP-seq literature identifying an unexpected number of transcription factor binding sites in regions distal to core promoter regions. Indeed, the majority of TP63 sites are associated with the characteristic histone modifications associated with enhancer regions, namely H3K4 mono and di-methylation and low tri-methylation, while a small proportion of TP63-binding sites are associated with transcription start sites and the cognate histone modifications, high levels of H3K27 acetylation and H3K4 trimethylation (49). This correlates with TP63 effecting long-range effects on transcription, and there is increasing evidence that many transcription factor-binding events elicit effects via long-range 3D interactions (71). This model for TP63 actions is supported by recent TP63 ChIP-seq data published while this article was in preparation, showing that in normal adult epidermal keratinocytes (NHEKs) TP63 regulated a long-range enhancer responsible for control of DLX4/5 as the causative genetic abnormality in a non-syndromic split-hand foot malformation (SHFM1) (72). Remarkably, comparing our dataset to the one in the Kouwenhoven paper (72), 5840 of the 7574 peaks identified in neonatal foreskin (77%) are present in their 11 368 reported binding regions underlining the robustness and reproducibility of the ChIP-seq strategies employed. The differences observed likely represent a combination of the variation in cut-offs applied, differences in antibodies, cell types used and depth of sequencing.

Given the potential for transcription factors such as TP63 to mediate long-range effects, the assignation of transcription factor binding sites to their cognate regulated genes has become a complex problem [reviewed in (73)]. For the purposes of this report, we defined a potential TP63 regulated gene as any gene, which falls within 25 kb of a TP63-binding site. This strategy allowed us to include genes that, for example, may be regulated by bi-directional promoters or overlapping regulatory regions, which could potentially be co-ordinately regulated by TP63. This revealed 6171 potential TP63 target genes from the Refseq database and included a large proportion of previously confirmed TP63 target genes. This annotation strategy is supported by Gene Ontology analysis using DAVID (50,51), which indicated greater enrichment for ‘epidermis development’ genes when we compared annotation distances of 5, 10 and 25 kb. In fact utilizing the recently published Genomic Enrichment of Annotations Tool (GREAT) (52), we observed increasing enrichment of the same ‘epidermis development’ term up to a maximum extensions of 500 kb.

The biological significance of this role in normal epidermal homeostasis is highlighted by observations in organotypic skin culture models, wherein depletion of all TP63 isoforms results in dramatic hypoproliferation and failure of stratification and differentiation (10,16). This hypoproliferative phenotype can be rescued by concomitant TP53 (16), CDKN1A, or p130 depletion (10); however, in both cases this is insufficient to restore differentiation. Furthermore, recent work in our laboratory indicates that TP63 plays an additional role in maintaining proliferative capacity through the positive regulation of the S-Phase Kinase-associated protein 2 (SKP2) (10). Thus, TP63 has been suggested to play a dual role in the epidermis maintaining proliferative capacity while also being required to regulate expression of genes required for differentiation, the latter in a TP53-independent manner. Comparison of our TP63-binding sites with recently published ChIP-seq studies of DNA damage-induced TP53-binding sites indicates that the majority of TP63 binding sites are not bound by TP53 in response to DNA damage. This stratification is further supported by GO analysis of annotated genes. Despite the majority of TP63 sites being devoid of TP53 binding in response to DNA damage, there are greater than 1300 sites predicted to be co-ordinately regulated by TP63 and TP53, which include many described TP63/TP53 targets including CDKN1A, MDM2, BAX, BBC3 (Puma), PMAIP1 (Noxa).

Intriguingly, initial analyses indicated enrichment of genes whose mutation is associated with cleft palate. This was of particular interest as we have recently shown an important role for TP63 in the TGF-beta-induced epithelial to mesenchymal transition in MCF10A basal breast myoepithelial cells (41), and TGF-beta and related signalling pathways have been shown to be critical for palatal fusion [reviewed in (74)]. In addition, GeneGO pathway analysis revealed TGF, Wnt and cytoskeletal remodelling as the most significantly enriched canonical pathways in our TP63 gene set, and these signal transduction pathways play well documented roles in palatal development [reviewed in (59,74)]. Finally, GO analysis of biological processes revealed significant enrichment for cell adhesion, development of epidermis/ectoderm and regulation of cell death and proliferation. These processes are necessary events in palatal fusion, which requires co-ordinate regulation of epithelial and mesenchymal growth, adhesion of the palatal shelves and removal of the midline cells by a combination of apoptosis and EMT for successful fusion [reviewed in (75)].

The integrative strategy employed allowed us to focus on a smaller number of genes in an unbiased manner for further validation, wherein we identify a number of the genes associated with both human syndromic and non-syndromic CL/P as direct targets of TP63. Interestingly, the AP-2 family of transcription factors that have been shown to be involved in development [reviewed in (76)], in particular AP-2alpha (TFAP2A), were also the most enriched predicted transcription factor binding sites within the region of DNA precipitated by TP63 antibodies in our in silico analyses. Potential co-operation with TP63 is suggested by overlapping phenotypes observed of patients with TP63 and TFAP2A which is the causative mutation in Branchiooculofacial syndrome (BOFS OMIM 113620) (22), a syndrome characterized by cutaneous/ocular anomalies, facial abnormalities cleft palate and like TP63 mutations, autosomal dominant inheritance. Furthermore, the TP63 regulated enhancer we and others have identified (65,66) in the IRF6 enhancer contains a SNP in the TFAP2A binding site, the absence of which in the minor allele predisposes to non-syndromic CL/P (77).

Given that TFAP2A has been shown to exhibit a similar epidermal pattern of expression to TP63, being predominantly expressed in the basal layers of the skin (68), we therefore characterized this potential co-regulation of genes by TFAP2A and TP63. We demonstrate that these factors can bind simultaneously to the same regions in vivo and that TP63 is required for TFAP2A recruitment to these regions. These interactions transcriptionally enhance TP63-mediated activation of the IRF6 and PVRL1 promoter/enhancer regions. We found that the deltaNTP63beta isoform exhibited greatest activation of the IRF6 region, whereas the deltaNTP63gamma isoform activated the PVRL1 region to the greatest extent, which correlates with observations in microarray studies where the deltaNTP63 isoforms have been individually overexpressed (78). In concordance with these studies, we show that the beta and gamma isoforms exhibit the greatest transcriptional activation capability, with the beta eliciting the largest response (79). The deltaNTP63 isoforms most frequently co-operated with AP-2 in H1299 cells and addition of TFAP2A and TFAP2C alone was sufficient to activate the IRF6 enhancer in HFKs, which express all three deltaNTP63 isoforms. This suggests that the deltaNTP63beta and gamma isoforms contribute to TP63 target gene regulation, but whether this is elicited by homo- or hetero- oligomers involving the predominant deltaNTP63alpha isoform is as yet undetermined.

Analysis of the role of TFAP2A and TFAP2C is complicated by embryonic lethality of knockouts (80); however, epithelial-specific conditional knockout mice reveal functional redundancy between the two factors, ablation of both revealing an essential role for AP-2 in induction of epidermal differentiation (81). Epithelial-specific ablation of TFAP2A alone was only sufficient to result in hyperproliferation and delayed differentiation (82), whereas epiblastic deletion of TFAP2C resulted only in a delay of skin development and barrier formation (83). Therefore, our results extend these studies to show that acute TFAP2A or TFAP2C depletion can inhibit epidermal differentiation and that this is likely due to cooperative regulation of TP63 target genes, adding additional levels of control of gene expression based on the temporal expression of the different AP-2 isoform.

In summary, the genome-wide identification of direct TP63 targets and their regulatory elements establishes a framework for the more detailed analysis of the TP63 function. In particular, for elucidation of the roles of individual TP63 isoforms, their interplay with other TP53 family members and co-factors, and how these are perturbed in disease. Furthermore, our observations in relation to TP63 and AP-2-mediated regulation of a network of genes implicated in cleft palate provides a valuable resource in the search for genetic determinants of CL/P and underlines the role of TP63 as a master regulator of development, growth and differentiation.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online: Supplementary Tables 1–10 and Supplementary Figures 1–7.

DATA ACCESSION

The ChIP-seq data from this study have been submitted to NCBI Gene Expression Omnibus (http://www.ncbi.nih.gov/geo/) under accession GSE32061 and the raw data passed to the sequence read archive (http://www.ncbi.nlm.nih.gov/sra).

FUNDING

The Medical Research Council project [G0700754]; Northern Ireland Cancer Translational Cancer Research Group [RRG/3284/05]; Breakthrough Breast Cancer. Funding for open access charge: MRC project grant.

Conflict of interest statement. None declared.